かっつのメモ帳

主に競プロ 時々日記

AOJ 1403 - Twin Trees Bros.

これは 帰ってきた AOJ-ICPC Advent Calendar 2022 の2日目です。

問題概要

問題のリンク

3次元平面上に $N$ 頂点の木が2つ与えられる。2つの木の頂点の対応付けであって、次の幾何学的変換を自由に行うことで、2つの木が完全に一致するようなものは何通り存在するか求めよ。

  • 各座標値に何らかの定数を足し合わせる平行移動

  • 全ての座標値に正の定数を掛け合わせる拡大縮小

  • x, y, z軸を中心とした回転移動

制約

  • $3 \leq N \leq 200$

  • $-1000 \leq x _ {i} , y _ {i} , z _ {i} \leq 1000$

楽しそう

解法

基本方針は、座標変換の方法を固定する→「2つの木に含まれる座標が等しい」かつ「各頂点の隣接リストが等しい」ことを調べ、一致判定を行うという流れです。

色々やりようはあると思っていて、多項式時間で解けるところまで落とし込んで、正しくそれを実装すれば大体通る気がしてます(多分)。

次に示すのは、木の葉ノードの1つに着目した変換方法の一例。

(1) 葉ノードを並行移動によって原点と一致させる

(2) 葉ノードとその親ノード間の距離を縮小により1にする

(3) 親ノードがxz平面上に来るように、z軸中心の回転移動を行う

(4) 親ノードがx軸上に来るように、y軸中心の回転移動を行う

(5) 親ノードに接続する葉ノード以外の点を(存在するならば)1つ選び、y軸上に来るように、x軸中心の回転移動を行う

(3)の変換のイメージ図

この変換は、葉ノードと(5)で登場する頂点を固定すれば一意に定まります。

木の片方は適当に変換すれば良く、もう片方について考えられる上記の変換を全て試して、得られた木が一致するか調べれば良いです。候補は $O(N^{2})$ で、木の同型判定に $O(N^{2})$ かけてもAOJでは余裕持って通すことが出来ました(0.09ms)。

自分の1WAの原因が重複を取り除けていない原因に依るもので、確かにちゃんと重複取り除こうとすると面倒になるなと分かったので、素朴にsetに全部入れて数えました。

追記(2022/12/02修正)

この変換は、葉ノードと(5)で登場する頂点を固定すれば一意に定まります。

↑これ大嘘でした。次の画像のように、固定した葉の親と、更にその唯一の隣接ノードが一直線上に並んだ時、変換方法が一意でなく、(自分の実装では)正しく判定できなくなっていました。

hackケースの1例*1(正しい答えは1)

6

0 0 0

1 0 0

2 0 0

3 0 0

3 1 0

3 2 0

1 2

2 3

3 4

4 5

5 6

0 0 0

1 0 0

2 0 0

3 0 0

3 -1 0

3 -2 0

7 8

8 9

9 10

10 11

11 12

パスグラフのようになっている部分をスキップすることで、同一の方針で正しく解くことが出来ます。以下の実装は修正済みの提出です。

実装

#pragma GCC optimize("Ofast")
#include <bits/stdc++.h>
using namespace std;
typedef long long int ll;

const double eps = 1e-8;
int sgn(double a){return a<-eps?-1:a>eps?1:0;}
int sgn(double a, double b){return sgn(a-b);}
struct point{
    double x,y,z;
    point() {}
    point(double x,double y,double z): x(x), y(y), z(z) {}
    point operator+ (const point& p){
        return point(x+p.x, y+p.y, z+p.z);
    }
    point operator- (const point& p){
        return point(x-p.x, y-p.y, z-p.z);
    }
    point operator/ (const double d){
        return point(x/d, y/d, z/d);
    }
    double dist(const point& p){
        return sqrt((x-p.x)*(x-p.x) + (y-p.y)*(y-p.y) + (z-p.z)*(z-p.z));
    }
};

int main(){
    cin.tie(nullptr);
    ios::sync_with_stdio(false);
    // input
    int n; cin >> n;
    vector<point> A(n);
    for(int i=0;i<n;i++){
        int x,y,z; cin >> x >> y >> z;
        A[i] = point(x,y,z);
    }
    vector<vector<int>> g1(n);
    for(int i=1;i<n;i++){
        int x,y; cin >> x >> y;
        x--; y--;
        g1[x].push_back(y);
        g1[y].push_back(x);
    }
    vector<point> B(n);
    for(int i=0;i<n;i++){
        int x,y,z; cin >> x >> y >> z;
        B[i] = point(x,y,z);
    }
    vector<vector<int>> g2(n);
    for(int i=1;i<n;i++){
        int x,y; cin >> x >> y;
        x -= n+1; y -= n+1;
        g2[x].push_back(y);
        g2[y].push_back(x);
    }

    // 木の葉ノードを基準に変形
    auto getTree = [&](vector<point> A, vector<vector<int>> &g, int leaf, int idx = -1)->vector<point>{
        // step1: 葉ノードを原点に平行移動させる
        {
            auto p = A[leaf];
            for(int i=0;i<n;i++){
                A[i] = A[i]-p;
            }
        }
        // step2: 葉ノードと親の距離を1にする
        int par = g[leaf][0];
        {
            double d = A[leaf].dist(A[par]);
            vector<point> nA(n);
            nA[leaf] = A[leaf];
            auto dfs = [&](auto self, int s, int p)->void{
                for(int t:g1[s]){
                    if(t == p) continue;
                    point v = (A[t]-A[s])/d;
                    nA[t] = nA[s]+v;
                    self(self, t, s);
                }
            };
            dfs(dfs, leaf, -1);
            swap(A, nA);
        }
        // step3: z軸中心に回転させ、親ノードをxz平面上へ
        {
            double theta = -atan2(A[par].y, A[par].x);
            for(int i=0;i<n;i++){
                double nx = A[i].x*cos(theta) - A[i].y*sin(theta);
                double ny = A[i].x*sin(theta) + A[i].y*cos(theta);
                A[i].x = nx, A[i].y = ny;
            }
        }
        // step4: y軸中心に回転させ、親ノードをx軸上へ
        {
            double theta = -atan2(A[par].z, A[par].x);
            for(int i=0;i<n;i++){
                double nx = A[i].x*cos(theta) - A[i].z*sin(theta);
                double nz = A[i].x*sin(theta) + A[i].z*cos(theta);
                A[i].x = nx, A[i].z = nz;
            }
        }
        // step5: x軸中心に回転させ、親ノードに接続する元の葉以外をy軸上へ
        // idx = -1 ならば任意、それ以外はA[idx]を合わせる
        if(idx == -1){
            while(1) {
                int nxtpar = -1;
                for (int t: g[par]) {
                    if (t == leaf) continue;
                    nxtpar = t;
                    if (abs(A[t].y) < eps and abs(A[t].z) < eps) continue;
                    idx = t;
                    break;
                }
                if (idx != -1 or nxtpar == -1) break;
                leaf = par;
                par = nxtpar;
            }
        }
        if(idx != -1){
            double theta = -atan2(A[idx].z, A[idx].y);
            for(int i=0;i<n;i++){
                double ny = A[i].y*cos(theta) - A[i].z*sin(theta);
                double nz = A[i].y*sin(theta) + A[i].z*cos(theta);
                A[i].y = ny, A[i].z = nz;
            }
        }
        return A;
    };

    int lA = -1;
    for(int i=0;i<n;i++){
        sort(g1[i].begin(), g1[i].end());
        if(g1[i].size() == 1) lA = i;
    }
    auto gA = getTree(A, g1, lA);

    // 比較用の関数
    set<vector<int>> st;
    auto compare = [&](vector<point> &A, vector<point> &B)->void{
        vector<int> Bidx(n,-1),Binv(n);
        // A[i] = B[Bidx[i]]
        for(int i=0;i<n;i++){
            for(int j=0;j<n;j++){
                if(sgn(A[i].x, B[j].x) or sgn(A[i].y, B[j].y) or sgn(A[i].z, B[j].z)) continue;
                Bidx[i] = j; Binv[j] = i; break;
            }
            if(Bidx[i] == -1) return;
        }
        // 隣接リストがそれぞれ一致しているか調べる
        for(int i=0;i<n;i++){
            auto v = g1[i], vv = g2[Bidx[i]];
            for(int &j:vv){
                j = Binv[j];
            }
            sort(vv.begin(), vv.end());
            if(v != vv) return;
        }
        st.insert(Bidx);
    };

    // AとBの比較
    for(int i=0;i<n;i++){
        if(g2[i].size() != 1) continue;
        // iは葉ノード
        int leaf = i;
        int par = g2[i][0];
        // step5でy軸に合わせるidxを全て試す
        int idx = -1;
        while(1) {
            int nxtpar = -1;
            for (int t: g2[par]) {
                if (t == leaf) continue;
                nxtpar = t;
                auto p = B[par]-B[leaf];
                auto pp = B[t]-B[leaf];
                // 同一直線上にある場合
                if(sgn(p.x*pp.y, pp.x*p.y) == 0 and sgn(p.y*pp.z, pp.y*p.z) == 0 and sgn(p.z*pp.x, pp.z*p.x) == 0) continue;
                auto gB = getTree(B, g2, i, t);
                compare(gA, gB);
                idx = t;
            }
            if (idx != -1 or nxtpar == -1) break;
            leaf = par;
            par = nxtpar;
        }
        if(idx == -1){
            auto gB = getTree(B, g2, i, -1);
            compare(gA, gB);
        }
    }

    // output
    cout << st.size() << endl;
}

*1:ウッキウキでAOJの全提出を試したら1つだけ落ちた。ほぼ自分しか引っかかってない嘘を加筆修正してる時間ちょっと恥ずかしい。