すごくシンプルだけど詰まる部分もたくさんありそうな問題

問題へのリンク

問題概要

二次元平面上に、2 組の $N$ 個の点集合 $S = \{(a_{0}, b_{0}), \dots, (a_{N-1}, b_{N-1})\}$ 、 $T = \{(c_{0}, d_{0}), \dots, (c_{N-1}, d_{N-1})\}$ がある。

$S$ に含まれる $N$ 個の点に対して、一律に

原点を中心とした回転をする (角度は任意)
平行移動をする (移動量は任意)

を実施することによって、 $T$ とピッタリ一致させることが可能かどうかを判定せよ。

制約

$1 \le N \le 100$
座標値の絶対値は $10$ 以下

考えたこと (解法 1)

前提として、各点を複素数として考えることにする。 $S$ の点 $i$ を複素数 $s_{i}$ と表すことにして、 $T$ の点 $i$ を複素数 $t_{i}$ と表すことにする。

計算量的には結構無茶苦茶やっても大丈夫そうなので、まずは考えるべき回転量や平行移動量が離散値になるように頑張ることにする。まずは、

$S$ の点 $s_{0}$ に対応する $T$ の点を $t_{x}$ である

と固定して考えることにした。この時点で $O(N)$ 通りとなる。この固定のいいところは「平行移動量を fix できる」というところだ。つまり、 $s_{0} = 0$ となるように $S$ の各点を平行移動し、 $t_{x} = 0$ となるように $T$ の各点を平行移動すると、「残りの $N-1$ 点同士が原点を中心とした回転によって一致するかどうかを判定する問題」となると言える。ここで、さらに、

$S$ の点 $s_{1}$ に対応する $T$ の点を $t_{y}$ である

と固定して考えよう。 $|t_{y}| = |s_{1}|$ でない場合はスキップすることにする。

このように決めると、 $S$ の各点をどれだけ回転させるべきかが一意に決まる。全体をまとめると、具体的には、次のように解ける。

Step 1： $S$ の各点から $s_{0}$ を引き、 $T$ の各点から $t_{x}$ を引く (このとき、 $s_{0} = 0$ 、 $t_{x} = 0$ となる)
各について、以下の探索を行う ( でない場合はスキップ)
- Step 2： $S$ の各点に対して、複素数 $\frac{t_{y}}{s_{1}}$ をかける
- Step 3： $S$ の各点が $T$ に一致するかどうかを判定する

Step 3 が少し悩むが、各点を { $x$ 座標値, $y$ 座標値} の辞書順でソートして ( $x$ 座標が一致する場合の誤差に注意)、前から順に一致するかを試していけば OK。

計算量は全体で $O(N^{3} \log N)$ となる。

コード

ここでは、複素数平面上の点を表す構造体ライブラリを使った。

#include<bits/stdc++.h>
using namespace std;

// basic settings
using DD = long double;
const long double PI = acosl(-1.0L);
constexpr long double EPS = 1e-10;   // to be set appropriately

// 複素数平面上の点を表すクラス
struct Point {
    DD x, y;
    
    // constructor
    constexpr Point() : x(0), y(0) {}
    constexpr Point(DD x, DD y) : x(x), y(y) {}
    
    // various functions
    constexpr Point conj() const {return Point(x, -y);}
    constexpr DD dot(const Point &r) const {return x * r.x + y * r.y;}
    constexpr DD cross(const Point &r) const {return x * r.y - y * r.x;}
    constexpr DD norm() const {return dot(*this);}
    constexpr long double abs() const {return sqrt(norm());}
    constexpr long double amp() const {
        long double res = atan2(y, x);
        if (res < 0) res += PI*2;
        return res;
    }
    constexpr bool eq(const Point &r) const {return (*this - r).abs() <= EPS;}
    constexpr Point rot90() const {return Point(-y, x);}
    constexpr Point rot(long double ang) const {
        return Point(cos(ang) * x - sin(ang) * y, sin(ang) * x + cos(ang) * y);
    }
    
    // arithmetic operators
    constexpr Point operator - () const {return Point(-x, -y);}
    constexpr Point operator + (const Point &r) const {return Point(*this) += r;}
    constexpr Point operator - (const Point &r) const {return Point(*this) -= r;}
    constexpr Point operator * (const Point &r) const {return Point(*this) *= r;}
    constexpr Point operator / (const Point &r) const {return Point(*this) /= r;}
    constexpr Point operator * (DD r) const {return Point(*this) *= r;}
    constexpr Point operator / (DD r) const {return Point(*this) /= r;}
    constexpr Point& operator += (const Point &r) {
        x += r.x, y += r.y;
        return *this;
    }
    constexpr Point& operator -= (const Point &r) {
        x -= r.x, y -= r.y;
        return *this;
    }
    constexpr Point& operator *= (const Point &r) {
        DD tx = x, ty = y;
        x = tx * r.x - ty * r.y;
        y = tx * r.y + ty * r.x;
        return *this;
    }
    constexpr Point& operator /= (const Point &r) {
        return *this *= r.conj() / r.norm();
    }
    constexpr Point& operator *= (DD r) {
        x *= r, y *= r;
        return *this;
    }
    constexpr Point& operator /= (DD r) {
        x /= r, y /= r;
        return *this;
    }

    // friend functions
    friend ostream& operator << (ostream &s, const Point &p) {
        return s << '(' << p.x << ", " << p.y << ')';
    }
    friend constexpr Point conj(const Point &p) {return p.conj();}
    friend constexpr DD dot(const Point &p, const Point &q) {return p.dot(q);}
    friend constexpr DD cross(const Point &p, const Point &q) {return p.cross(q);}
    friend constexpr DD norm(const Point &p) {return p.norm();}
    friend constexpr long double abs(const Point &p) {return p.abs();}
    friend constexpr long double amp(const Point &p) {return p.amp();}
    friend constexpr bool eq(const Point &p, const Point &q) {return p.eq(q);}
    friend constexpr Point rot90(const Point &p) {return p.rot90();}
    friend constexpr Point rot(const Point &p, long long ang) {return p.rot(ang);}
};

int main () {
    int N;
    cin >> N;
    vector<Point> s(N), t(N);
    for (int i = 0; i < N; ++i) cin >> s[i].x >> s[i].y;
    for (int i = 0; i < N; ++i) cin >> t[i].x >> t[i].y;
    
    // 例外処理
    if (N == 1) {
        cout << "Yes" << endl;
        return 0;
    }

    // s[0] と t[x] が対応するとする
    bool res = false;
    for (int x = 0; x < N; ++x) {
        // 平行移動させて、s[0] と t[x] が原点に来るようにする
        vector<Point> s2(N), t2(N);
        for (int i = 0; i < N; ++i) s2[i] = s[i] - s[0], t2[i] = t[i] - t[x];
        
        // もう 1 点固定する
        for (int y = 0; y < N; ++y) {
            if (x == y) continue;
            if (abs(abs(t2[y]) - abs(s2[1])) > EPS) continue;
            
            // S の各点を回転させる
            vector<Point> s3(N), t3(N);
            Point kaiten = t2[y] / s2[1];
            for (int i = 0; i < N; ++i) {
                s3[i] = s2[i] * kaiten;
                t3[i] = t2[i];
            }
            
            // s3 と t3 が一致するかを判定する
            bool same = true;
            auto cmp = [&](Point a, Point b) -> bool {
                return (abs(a.x - b.x) > EPS ? a.x < b.x : a.y < b.y);
            };
            sort(s3.begin(), s3.end(), cmp);
            sort(t3.begin(), t3.end(), cmp);
            for (int i = 0; i < N; ++i) if (!eq(s3[i], t3[i])) same = false;
            if (same) res = true;
        }
    }
    cout << (res ? "Yes" : "No") << endl;
}

解法 (2)：点の相対的位置関係を符号化する

まず、

$S$ の点 $s_{0}$ に対応する $T$ の点を $t_{x}$

と固定して考えるところまでは、解法 (1) と同じ。さらに、 $S$ の各点から $s_{0}$ を引き、 $T$ の各点から $t_{x}$ を引くところまでも同じ。こうして、 $s_{0}$ と $t_{x}$ がともに原点になるように平行移動を完成させたあと、 $S$ と $T$ とが回転で移り合うかどうかを判定する方法を考えるところまでは同じ。

ここから先が解法 (1) と異なる。実際に回転しない方法がある。一般に、3 点 $0, s_{1}, s_{2}$ と、 $0, t_{y}, t_{z}$ が、回転によって同じ配置になるかどうかは次のように判定できる。

$|s_{1}| = |t_{y}|$ かつ $|s_{2}| = |t_{z}|$ である
「 $s_{1}, s_{2}$ の内積 = $t_{y}, t_{z}$ の内積」である
「 $s_{1}, s_{2}$ の外積 = $t_{y}, t_{z}$ の外積」である

最後の 3. が必要なのは、裏返しの位置にある関係を防ぐためだ。このことを踏まえると、次の解法が考えられる。

Step 1： $S$ の各点から $s_{0}$ を引き、 $T$ の各点から $t_{x}$ を引く (このとき、 $s_{0} = 0$ 、 $t_{x} = 0$ となる)
Step 2： $S, T$ ( $s_{0}$ と $t_{x}$ を除く) の各点を原点周りに偏角ソートする
Step 3： $S$ に対して、 $(|s_{i}|$ , $s_{i}, s_{i+1}$ の内積, $s_{i}, s_{i+1}$ の外積) の 3 つ組を $i$ に対して並べて得られる $3N-3$ 次元ベクトル $p$ と、 $T$ に対して同様の手続きで得られる $3N-3$ 次元ベクトル $q$ を作る
Step 4：ある $k$ が存在して、 $p$ を $3k$ だけ巡回シフトすることで $q$ と一致するかどうかを調べる

この解法のいいところは、整数値計算のみで完結できる ( $(|s_{i}|$ の部分は二乗すればよい)。Step 4 に $O(N^{2})$ の計算量を要することから、全体の計算量は $O(N^{3})$ となる。

なお、Step 4 の部分は Z-algorithm や Suffix Array (の LCP 配列) などを用いることで $O(N)$ や $O(N \log N)$ に改良することもできる。その場合、全体の計算量は、今度は Step 3 の偏角ソートがボトルネックとなって、 $O(N^{2} \log N)$ となる。

コード

ここでは、Step 4 において $O(N^{2})$ の計算量を要する単純な全探索した場合のコード (全体で $O(N^{3})$ の計算量) を示す。

#include<bits/stdc++.h>
using namespace std;

// basic settings
using DD = long long;
const long double PI = acosl(-1.0L);
constexpr long double EPS = 1e-10;    // to be set appropriately

// Point or Vector
struct Point {
    DD x, y;
    
    // constructor
    constexpr Point() : x(0), y(0) {}
    constexpr Point(DD x, DD y) : x(x), y(y) {}
    
    // various functions
    constexpr Point conj() const {return Point(x, -y);}
    constexpr DD dot(const Point &r) const {return x * r.x + y * r.y;}
    constexpr DD cross(const Point &r) const {return x * r.y - y * r.x;}
    constexpr DD norm() const {return dot(*this);}
    constexpr long double abs() const {return sqrt(norm());}
    constexpr long double amp() const {
        long double res = atan2(y, x);
        if (res < 0) res += PI*2;
        return res;
    }
    constexpr bool eq(const Point &r) const {return (*this - r).abs() <= EPS;}
    constexpr Point rot90() const {return Point(-y, x);}
    constexpr Point rot(long double ang) const {
        return Point(cos(ang) * x - sin(ang) * y, sin(ang) * x + cos(ang) * y);
    }
    
    // arithmetic operators
    constexpr Point operator - () const {return Point(-x, -y);}
    constexpr Point operator + (const Point &r) const {return Point(*this) += r;}
    constexpr Point operator - (const Point &r) const {return Point(*this) -= r;}
    constexpr Point operator * (const Point &r) const {return Point(*this) *= r;}
    constexpr Point operator / (const Point &r) const {return Point(*this) /= r;}
    constexpr Point operator * (DD r) const {return Point(*this) *= r;}
    constexpr Point operator / (DD r) const {return Point(*this) /= r;}
    constexpr Point& operator += (const Point &r) {
        x += r.x, y += r.y;
        return *this;
    }
    constexpr Point& operator -= (const Point &r) {
        x -= r.x, y -= r.y;
        return *this;
    }
    constexpr Point& operator *= (const Point &r) {
        DD tx = x, ty = y;
        x = tx * r.x - ty * r.y;
        y = tx * r.y + ty * r.x;
        return *this;
    }
    constexpr Point& operator /= (const Point &r) {
        return *this *= r.conj() / r.norm();
    }
    constexpr Point& operator *= (DD r) {
        x *= r, y *= r;
        return *this;
    }
    constexpr Point& operator /= (DD r) {
        x /= r, y /= r;
        return *this;
    }

    // friend functions
    friend ostream& operator << (ostream &s, const Point &p) {
        return s << '(' << p.x << ", " << p.y << ')';
    }
    friend constexpr Point conj(const Point &p) {return p.conj();}
    friend constexpr DD dot(const Point &p, const Point &q) {return p.dot(q);}
    friend constexpr DD cross(const Point &p, const Point &q) {return p.cross(q);}
    friend constexpr DD norm(const Point &p) {return p.norm();}
    friend constexpr long double abs(const Point &p) {return p.abs();}
    friend constexpr long double amp(const Point &p) {return p.amp();}
    friend constexpr bool eq(const Point &p, const Point &q) {return p.eq(q);}
    friend constexpr Point rot90(const Point &p) {return p.rot90();}
    friend constexpr Point rot(const Point &p, long long ang) {return p.rot(ang);}
};

// 偏角ソート
void amp_sort(vector<Point> &v) {
    sort(v.begin(), v.end(), [&](Point p, Point q) {
        return (abs(amp(p) - amp(q)) > EPS ? amp(p) < amp(q) : norm(p) < norm(q));
    });
}

int main () {
    int N;
    cin >> N;
    vector<Point> s(N), t(N);
    for (int i = 0; i < N; ++i) cin >> s[i].x >> s[i].y;
    for (int i = 0; i < N; ++i) cin >> t[i].x >> t[i].y;
    
    // 例外処理
    if (N == 1) {
        cout << "Yes" << endl;
        return 0;
    }

    // s[0] と t[x] が対応するとする
    bool res = false;
    for (int x = 0; x < N; ++x) {
        // 平行移動させて、s[0] と t[x] が原点に来るようにする
        vector<Point> s2, t2;
        for (int i = 0; i < N; ++i) {
            if (i != 0) s2.push_back(s[i] - s[0]);
            if (i != x) t2.push_back(t[i] - t[x]);
        }

        // 偏角ソート
        amp_sort(s2);
        amp_sort(t2);

        // ベクトルを作る
        vector<long long> p, q;
        for (int i = 0; i < s2.size(); ++i) {
            p.push_back(norm(s2[i]));
            p.push_back(dot(s2[i], s2[(i+1)%s2.size()]));
            p.push_back(cross(s2[i], s2[(i+1)%s2.size()]));
            q.push_back(norm(t2[i]));
            q.push_back(dot(t2[i], t2[(i+1)%t2.size()]));
            q.push_back(cross(t2[i], t2[(i+1)%t2.size()]));
        }

        // 3 の倍数シフトを試す
        for (int k = 0; k < p.size(); k += 3) {
            bool same = true;
            for (int i = 0; i < p.size(); ++i) {
                if (p[(i+k)%p.size()] != q[i]) same = false;
            }
            if (same) res = true;
        }
    }
    cout << (res ? "Yes" : "No") << endl;
}

解法 (3)：重心を合わせる

解法 (1)、解法 (2) はいずれも、 $s_{0}$ に対応する $t_{x}$ を決めた上で判定した。これによって平行移動量を fix することができて、残りは回転移動のみを考えれば良いという戦法だ。しかし、この部分ですでに $O(N)$ 通りの探索となっている。

実は、平行移動量を $O(1)$ で fix できる方法がある！　点集合 $S, T$ それぞれの重心を求め、それらを予め原点に来るように平行移動すればよいのだ。残りは「回転して一致するか」を判定する問題となる。

つまり重心を活用することによって、解法 (1)(2) それぞれについて、計算量のオーダーの次数を 1 ずつ下げることができる。解法 (1) は $O(N^{2} \log N)$ の計算量となるし、解法 (2) は $O(N^{2})$ または $O(N \log N)$ の計算量となる。なんと、 $N \le 2 \times 10^{5}$ のような制約でも解けることが分かる！

コード (解法 (2) の高速化)

ここでは、解法 (2) の Step 4 において、Suffix Array を用いた解法を示す。全体として $O(N \log N)$ の計算量となる。ただし、定数倍が重たいので、提出してみたら解法 (1) の方がかえって速かった。

#include<bits/stdc++.h>
using namespace std;

// Suffix Array
template<class MeetSemiLattice> struct SparseTable {
    vector<vector<MeetSemiLattice> > dat;
    vector<int> height;
    
    SparseTable() { }
    SparseTable(const vector<MeetSemiLattice> &vec) { init(vec); }
    void init(const vector<MeetSemiLattice> &vec) {
        int n = (int)vec.size(), h = 1;
        while ((1<<h) < n) ++h;
        dat.assign(h, vector<MeetSemiLattice>(1<<h));
        height.assign(n+1, 0);
        for (int i = 2; i <= n; i++) height[i] = height[i>>1]+1;
        for (int i = 0; i < n; ++i) dat[0][i] = vec[i];
        for (int i = 1; i < h; ++i)
            for (int j = 0; j < n; ++j)
                dat[i][j] = min(dat[i-1][j], dat[i-1][min(j+(1<<(i-1)),n-1)]);
    }
    
    MeetSemiLattice get(int a, int b) {
        return min(dat[height[b-a]][a], dat[height[b-a]][b-(1<<height[b-a])]);
    }
};

template<class Str> struct SuffixArray {
    // data
    Str str;
    vector<int> sa;    // sa[i] : the starting index of the i-th smallest suffix (i = 0, 1, ..., n)
    vector<int> rank;  // rank[sa[i]] = i
    vector<int> lcp;   // lcp[i]: the lcp of sa[i] and sa[i+1] (i = 0, 1, ..., n-1)
    SparseTable<int> st;  // use for calcultating lcp(i, j)

    // getter
    int& operator [] (int i) {
        return sa[i];
    }
    vector<int> get_sa() { return sa; }
    vector<int> get_rank() { return rank; }
    vector<int> get_lcp() { return lcp; }

    // constructor
    SuffixArray() {}
    SuffixArray(const Str& str_) : str(str_) {
        build_sa();
    }
    void init(const Str& str_) {
        str = str_;
        build_sa();
    }
    void build_sa() {
        vector<int> s;
        unordered_map<int,int> dict;
        for (int i = 0; i < (int)str.size(); ++i) {
            if (!dict.count(str[i])) dict[str[i]] = dict.size();
        }
        for (int i = 0; i < (int)str.size(); ++i) {
            s.push_back(dict[str[i]] + 1);
        }
        s.push_back(0);
        sa = sa_is(s, dict.size());
        calcLCP(s);
        buildSparseTable();
    }

    // SA-IS
    // upper: # of characters
    vector<int> sa_is(vector<int> &s, long long upper = 100000) {
        int N = (int)s.size();
        if (N == 0) return {};
        else if (N == 1) return {0};
        else if (N == 2) {
            if (s[0] < s[1]) return {0, 1};
            else return {1, 0};
        }

        vector<int> isa(N);
        vector<bool> ls(N, false);
        for (int i = N - 2; i >= 0; --i) {
            ls[i] = (s[i] == s[i + 1]) ? ls[i + 1] : (s[i] < s[i + 1]);
        }
        vector<int> sum_l(upper + 1, 0), sum_s(upper + 1, 0);
        for (int i = 0; i < N; ++i) {
            if (!ls[i]) ++sum_s[s[i]];
            else ++sum_l[s[i] + 1];
        }
        for (int i = 0; i <= upper; ++i) {
            sum_s[i] += sum_l[i];
            if (i < upper) sum_l[i + 1] += sum_s[i];
        }

        auto induce = [&](const vector<int> &lms) -> void {
            fill(isa.begin(), isa.end(), -1);
            vector<int> buf(upper + 1);
            copy(sum_s.begin(), sum_s.end(), buf.begin());
            for (auto d: lms) {
                if (d == N) continue;
                isa[buf[s[d]]++] = d;
            }
            copy(sum_l.begin(), sum_l.end(), buf.begin());
            isa[buf[s[N - 1]]++] = N - 1;
            for (int i = 0; i < N; ++i) {
                int v = isa[i];
                if (v >= 1 && !ls[v - 1]) {
                    isa[buf[s[v - 1]]++] = v - 1;
                }
            }
            copy(sum_l.begin(), sum_l.end(), buf.begin());
            for (int i = N - 1; i >= 0; --i) {
                int v = isa[i];
                if (v >= 1 && ls[v - 1]) {
                    isa[--buf[s[v - 1] + 1]] = v - 1;
                }
            }
        };
            
        vector<int> lms, lms_map(N + 1, -1);
        int M = 0;
        for (int i = 1; i < N; ++i) {
            if (!ls[i - 1] && ls[i]) {
                lms_map[i] = M++;
            }
        }
        lms.reserve(M);
        for (int i = 1; i < N; ++i) {
            if (!ls[i - 1] && ls[i]) {
                lms.push_back(i);
            }
        }
        induce(lms);

        if (M) {
            vector<int> lms2;
            lms2.reserve(isa.size());
            for (auto v: isa) {
                if (lms_map[v] != -1) lms2.push_back(v);
            }
            int rec_upper = 0;
            vector<int> rec_s(M);
            rec_s[lms_map[lms2[0]]] = 0;
            for (int i = 1; i < M; ++i) {
                int l = lms2[i - 1], r = lms2[i];
                int nl = (lms_map[l] + 1 < M) ? lms[lms_map[l] + 1] : N;
                int nr = (lms_map[r] + 1 < M) ? lms[lms_map[r] + 1] : N;
                bool same = true;
                if (nl - l != nr - r) same = false;
                else {
                    while (l < nl) {
                        if (s[l] != s[r]) break;
                        ++l, ++r;
                    }
                    if (l == N || s[l] != s[r]) same = false;
                }
                if (!same) ++rec_upper;
                rec_s[lms_map[lms2[i]]] = rec_upper;
            }
            auto rec_sa = sa_is(rec_s, rec_upper);

            vector<int> sorted_lms(M);
            for (int i = 0; i < M; ++i) {
                sorted_lms[i] = lms[rec_sa[i]];
            }
            induce(sorted_lms);
        }
        return isa;
    }

    // find min id that str.substr(sa[id]) >= T
    int lower_bound(const Str& T) {
        int left = -1, right = sa.size();
        while (right - left > 1) {
            int mid = (left + right) / 2;
            if (str.compare(sa[mid], string::npos, T) < 0)
                left = mid;
            else
                right = mid;
        }
        return right;
    }

    // find min id that str.substr(sa[id], T.size()) > T
    int upper_bound(const Str& T) {
        int left = -1, right = sa.size();
        while (right - left > 1) {
            int mid = (left + right) / 2;
            if (str.compare(sa[mid], T.size(), T) <= 0)
                left = mid;
            else
                right = mid;
        }
        return right;
    }

    // search
    bool is_contain(const Str& T) {
        int lb = lower_bound(T);
        if (lb >= sa.size()) return false;
        return str.compare(sa[lb], T.size(), T) == 0;
    }

    // find lcp
    void calcLCP(const vector<int> &s) {
        int N = (int)s.size();
        rank.assign(N, 0), lcp.assign(N - 1, 0);
        for (int i = 0; i < N; ++i) rank[sa[i]] = i;
        int h = 0;
        for (int i = 0; i < N - 1; ++i) {
            int pi = sa[rank[i] - 1];
            if (h > 0) --h;
            for (; pi + h < N && i + h < N; ++h) {
                if (s[pi + h] != s[i + h]) break;
            }
            lcp[rank[i] - 1] = h;
        }
    }
    
    // build sparse table for calculating lcp
    void buildSparseTable() {
        st.init(lcp);
    }

    // calc lcp of str.sutstr(a) and str.substr(b)
    int getLCP(int a, int b) {
        return st.get(min(rank[a], rank[b]), max(rank[a], rank[b]));
    }

    // debug
    void dump() {
        for (int i = 0; i < sa.size(); ++i) {
            cout << i << ": " << sa[i] << ", " << str.substr(sa[i]) << endl;
        }
    }
};


// 複素数平面上の点を管理するライブラリ
using DD = long long;
constexpr long double EPS = 1e-10;   // to be set appropriately
const long double PI = acosl(-1.0L);

// Point
struct Point {
    DD x, y;
    
    // constructor
    constexpr Point() : x(0), y(0) {}
    constexpr Point(DD x, DD y) : x(x), y(y) {}
    
    // various functions
    constexpr Point conj() const {return Point(x, -y);}
    constexpr DD dot(const Point &r) const {return x * r.x + y * r.y;}
    constexpr DD cross(const Point &r) const {return x * r.y - y * r.x;}
    constexpr DD norm() const {return dot(*this);}
    constexpr long double abs() const {return sqrt(norm());}
    constexpr long double amp() const {
        long double res = atan2(y, x);
        if (res < 0) res += PI*2;
        return res;
    }
    constexpr bool eq(const Point &r) const {return (*this - r).abs() <= EPS;}
    constexpr Point rot90() const {return Point(-y, x);}
    constexpr Point rot(long double ang) const {
        return Point(cos(ang) * x - sin(ang) * y, sin(ang) * x + cos(ang) * y);
    }
    
    // arithmetic operators
    constexpr Point operator - () const {return Point(-x, -y);}
    constexpr Point operator + (const Point &r) const {return Point(*this) += r;}
    constexpr Point operator - (const Point &r) const {return Point(*this) -= r;}
    constexpr Point operator * (const Point &r) const {return Point(*this) *= r;}
    constexpr Point operator / (const Point &r) const {return Point(*this) /= r;}
    constexpr Point operator * (DD r) const {return Point(*this) *= r;}
    constexpr Point operator / (DD r) const {return Point(*this) /= r;}
    constexpr Point& operator += (const Point &r) {
        x += r.x, y += r.y;
        return *this;
    }
    constexpr Point& operator -= (const Point &r) {
        x -= r.x, y -= r.y;
        return *this;
    }
    constexpr Point& operator *= (const Point &r) {
        DD tx = x, ty = y;
        x = tx * r.x - ty * r.y;
        y = tx * r.y + ty * r.x;
        return *this;
    }
    constexpr Point& operator /= (const Point &r) {
        return *this *= r.conj() / r.norm();
    }
    constexpr Point& operator *= (DD r) {
        x *= r, y *= r;
        return *this;
    }
    constexpr Point& operator /= (DD r) {
        x /= r, y /= r;
        return *this;
    }
    
    // comparison operators
    constexpr bool operator == (const Point &r) const {return eq(r);}
    constexpr bool operator != (const Point &r) const {return !eq(r);}
    constexpr bool operator < (const Point &r) const {
        return (::abs(x - r.x) > EPS ? x < r.x : y < r.y);
    }
    
    // friend functions
    friend ostream& operator << (ostream &s, const Point &p) {
        return s << '(' << p.x << ", " << p.y << ')';
    }
    friend constexpr Point conj(const Point &p) {return p.conj();}
    friend constexpr DD dot(const Point &p, const Point &q) {return p.dot(q);}
    friend constexpr DD cross(const Point &p, const Point &q) {return p.cross(q);}
    friend constexpr DD norm(const Point &p) {return p.norm();}
    friend constexpr long double abs(const Point &p) {return p.abs();}
    friend constexpr long double amp(const Point &p) {return p.amp();}
    friend constexpr bool eq(const Point &p, const Point &q) {return p.eq(q);}
    friend constexpr Point rot90(const Point &p) {return p.rot90();}
    friend constexpr Point rot(const Point &p, long long ang) {return p.rot(ang);}
};

// 偏角ソート
void amp_sort(vector<Point> &v) {
    sort(v.begin(), v.end(), [&](Point p, Point q) {
        return (abs(amp(p) - amp(q)) > EPS ? amp(p) < amp(q) : norm(p) < norm(q));
    });
}

int main () {
    // 入力受け取り
    int N;
    cin >> N;
    vector<Point> s(N), t(N);
    Point gs(0, 0), gt(0, 0);
    for (int i = 0; i < N; ++i) {
        cin >> s[i].x >> s[i].y;
        s[i] *= N;
        gs += s[i];
    }
    for (int i = 0; i < N; ++i) {
        cin >> t[i].x >> t[i].y;
        t[i] *= N;
        gt += t[i];
    }
    gs /= N, gt /= N;  // 重心
    
    // 例外処理
    if (N == 1) {
        cout << "Yes" << endl;
        return 0;
    }
    
    // 平行移動させて、S, T それぞれの重心が原点に来るようにする。
    // ただし、重心に一致する点がある場合は除去しておく。
    vector<Point> s2, t2;
    for (int i = 0; i < N; ++i) {
        if (s[i] != gs) s2.push_back(s[i] - gs);
        if (t[i] != gt) t2.push_back(t[i] - gt);
    }
    if (s2.size() != t2.size()) {
        cout << "No" << endl;
        return 0;
    }
    
    // 偏角ソート
    amp_sort(s2), amp_sort(t2);
    
    // ベクトルを作る (S 側と T 側を連結させてしまう)
    vector<long long> v;
    for (int i = 0; i < (int)s2.size(); ++i) {
        v.push_back(norm(s2[i]));
        v.push_back(dot(s2[i], s2[(i+1)%s2.size()]));
        v.push_back(cross(s2[i], s2[(i+1)%s2.size()]));
    }
    for (int i = 0; i < (int)s2.size(); ++i) {
        v.push_back(norm(t2[i]));
        v.push_back(dot(t2[i], t2[(i+1)%t2.size()]));
        v.push_back(cross(t2[i], t2[(i+1)%t2.size()]));
    }
    
    // Suffix Array を構築して、3 の倍数シフトを試す
    bool res = false;
    SuffixArray<vector<long long>> sa(v);
    int len = (int)s2.size() * 3;
    for (int k = 0; k < len; k += 3) {
        if (sa.getLCP(0, k+len) >= len-k && sa.getLCP(len-k, len) >= k) res = true;
    }
    
    cout << (res ? "Yes" : "No") << endl;
}

けんちょんの競プロ精進記録

競プロの精進記録や小ネタを書いていきます

AtCoder ABC 207 D - Congruence Points (2D, 黄色, 400 点)

問題概要

制約

考えたこと (解法 1)

コード

解法 (2)：点の相対的位置関係を符号化する

コード

解法 (3)：重心を合わせる

コード (解法 (2) の高速化)