mayoko’s diary

プロコンとかいろいろ。

何かを解きました

何か楽しそうなキャンペーンをやっていたので応募しました。プレゼントがほしかったのでみんなには内緒でこっそりやってました(といっても知ってる人多そう)。

問題

recruit.elysium.co.jp

(もしかしたら このページが消えるかもしれないので問題文も書いておきます)
平面上に n 個の赤い点と n 個の青い点が与えられ(これら 2n 個の点はどの 3 点を取っても同一直線上にないことが保証される), 赤い点と青い点を結んで n 個の線分を作ることを考える。このとき,
(1) うまく線分を作ると, n 個の線分がすべて交差しないようにできることを示せ。
(2) n 個の線分がすべて交差しないようにするアルゴリズムを作れ。

10 <= n <= 10^7

解法

(2) の話だけします(といっても (1) の証明を使って (2) を解いたので少しだけその話もします)。僕の解法は計算量の期待値が O(N log^2 N) です。n = 10^7 に突っ込むと 5 分くらい結果が返ってきません。

(1) の証明の概要は, 数学的帰納法です。
(a) n = 1 の時答えがあるのは明らか

(b) n = 1, 2, ..., k で正しいとすると, n = k+1 でも正しいことを示す。平面上にある直線 l が存在して, 直線 l の上側の点は 赤い点, 青い点が x 個ずつ, 下側では k+1-x 個ずつあるようにすることができる(0 < x < k+1)。問題はこの直線を介して別の部分問題に分かれるので, n = k+1 でも OK

というような感じです(直線が存在することの証明が若干難しいです)。これと同じようにして, 問題をどんどん分割していくことを考えましょう。

上記の性質を満たす直線は, 少し動かしてもその性質を変えません。性質が変わるのは, この直線がある点の上に乗ったときです。よって, 考える直線は, ある点を通るものとしても良いです。この点 r から, ほかの点を偏角ソートします。また, 最初に r を通り x 軸に平行な直線を考え, これの上側にある点と下側にある頂点を求めておきます。

後は直線を少しずつ回しながら, 直線の上にある頂点の数を更新していきます。で, 赤い点と青い点が等しい数になったら終わりにして, 部分問題を解かせます。

点 r をランダムに選べば, 毎回の計算で頂点の数が半分になることが期待できる(ホントか) ので, 深さは O(logN) になりそうです。また, 各 dfs では持っている点のサイズ sz に対して, 計算量 O(sz*log(sz)) です。なので, 全体の計算量はだいたい O(N log^2 N) です。

この問題では座標がかなり大きくなりうる(|x| <= 10^9, |y| <= 10^9) ので, 偏角ソートを atan2 の double で持つと怪しいです。なので, 下では偏角ソートを外積でやるやつをやっています。

typedef pair<int, int> pii;
typedef vector<int> vi;
typedef long long ll;
std::random_device rnd;
std::mt19937 mt(rnd());

// geometry
struct Pint {
    int x;
    int y;
    Pint() {}
    Pint(int x, int y) : x(x), y(y) {}
    Pint operator+(Pint p) const {return Pint(x+p.x, y+p.y);}
    Pint operator-(Pint p) const {return Pint(x-p.x, y-p.y);}
    ll det(Pint p) const {return (ll)x*p.y - (ll)y*p.x;} // 外積
    bool operator<(const Pint& rhs) const {
        return det(rhs) > 0;
    }
};

inline int sign(ll x) {
    if (x > 0) return 1;
    else if (x == 0) return 0;
    else return -1;
}

// 線分p1-p2 と線分q1-q2 が交差しないか判定
bool intersect(Pint p1, Pint p2, Pint q1, Pint q2) {
    if (sign((q1-p1).det(q1-p2)) * sign((q2-p1).det(q2-p2)) >= 0) return false;
    if (sign((p1-q1).det(p1-q2)) * sign((p2-q1).det(p2-q2)) >= 0) return false;
    return true;
}

const int MAXN = 10001000;
Pint red[MAXN], blue[MAXN];
inline Pint trans(const pii& p) {
    Pint ret;
    if (p.first == 0) ret = red[p.second];
    else ret = blue[p.second];
    return ret;
};

// 頂点数
int n;

// 分割統治で解を求める
// 計算量の期待値は O(Nlog^2N)
void dfs(vi& r, vi& b) {
    int sz = r.size();
    assert(sz == (int)(b.size()));
    // サイズが 2 以下の場合は直接求める
    if (sz == 2) {
        if (!intersect(red[r[0]], blue[b[0]], red[r[1]], blue[b[1]])) {
            cout << r[0] << " " << b[0] << endl;
            cout << r[1] << " " << b[1] << endl;
        } else {
            assert(!intersect(red[r[0]], blue[b[1]], red[r[1]], blue[b[0]]));
            cout << r[0] << " " << b[1] << endl;
            cout << r[1] << " " << b[0] << endl;
        }
    }
    if (sz == 1) {
        cout << r[0] << " " << b[0] << endl;
    }
    if (sz <= 2) return;
    // 部分問題に分ける
    int rIndex = mt()%sz;
    const Pint rndPnt = red[r[rIndex]];
    // 上側にある頂点と下側にある頂点
    // first: color second: index
    vector<pii> upVertex, downVertex;
    // 上側にある頂点の集合
    // first: color second: index
    set<pii> S;
    // red で上側にある頂点数, blue で上側にある頂点数
    int rup = 1, bup = 0;
    S.insert(pii(0, r[rIndex]));
    // red
    for (int i = 0; i < sz; i++) if (i != rIndex) {
        Pint diff = red[r[i]] - rndPnt;
        if (diff.y < 0) {
            downVertex.emplace_back(0, r[i]);
        } else {
            rup++;
            S.insert(pii(0, r[i]));
            upVertex.emplace_back(0, r[i]);
        }
    }
    // blue
    for (int i = 0; i < sz; i++) {
        Pint diff = blue[b[i]] - rndPnt;
        if (diff.y < 0) {
            downVertex.emplace_back(1, b[i]);
        } else {
            bup++;
            S.insert(pii(1, b[i]));
            upVertex.emplace_back(1, b[i]);
        }
    }
    // 偏角ソートの代わりに外積でソート
    auto cmp = [&](const pii& lhs, const pii& rhs) {
        Pint lp = trans(lhs), rp = trans(rhs);
        lp = lp-rndPnt;
        rp = rp-rndPnt;
        return lp < rp;
    };
    sort(upVertex.begin(), upVertex.end(), cmp);
    sort(downVertex.begin(), downVertex.end(), cmp);
    int upSize = upVertex.size();
    int downSize = downVertex.size();
    int upCnt = 0, downCnt = 0;
    while (upCnt < upSize && downCnt < downSize) {
        if (rup == bup && rup < sz) break;
        // 上と下どちらを先に処理すべきなのかを決める
        Pint upDiff = trans(upVertex[upCnt]), downDiff = trans(downVertex[downCnt]);
        upDiff = upDiff - rndPnt, downDiff = downDiff - rndPnt;
        downDiff.x *= -1; downDiff.y *= -1;
        if (downDiff < upDiff) {
            // down のほうを先に調べる
            pii tmp = downVertex[downCnt++];
            if (tmp.first == 0) {
                rup++;
            } else {
                bup++;
            }
            S.insert(tmp);
        } else {
            // up のほうを先に調べる
            pii tmp = upVertex[upCnt++];
            if (tmp.first == 0) {
                rup--;
            } else {
                bup--;
            }
            S.erase(tmp);
        }
    }
    // 上側で残ってるのがあれば処理
    while (upCnt < upSize) {
        if (rup == bup && rup < sz) break;
        pii tmp = upVertex[upCnt++];
        if (tmp.first == 0) {
            rup--;
        } else {
            bup--;
        }
        S.erase(tmp);
    }
    // 下側も同様
    while (downCnt < downSize) {
        if (rup == bup && rup < sz) break;
        pii tmp = downVertex[downCnt++];
        if (tmp.first == 0) {
            rup++;
        } else {
            bup++;
        }
        S.insert(tmp);
    }
    // 失敗してたらやり直し
    if (!(rup == bup && rup < sz)) {
        upVertex.clear(); downVertex.clear(); S.clear();
        return dfs(r, b);
    }
    // 成功してたら問題を分割
    vi nr1, nr2, nb1, nb2;
    for (int i = 0; i < sz; i++) {
        if (S.find(pii(0, r[i])) == S.end()) {
            nr1.push_back(r[i]);
        } else {
            nr2.push_back(r[i]);
        }
        if (S.find(pii(1, b[i])) == S.end()) {
            nb1.push_back(b[i]);
        } else {
            nb2.push_back(b[i]);
        }
    }
    r.clear(); b.clear(); S.clear(); upVertex.clear(); downVertex.clear();
    dfs(nr1, nb1);
    dfs(nr2, nb2);
}

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    cin >> n;
    for (int i = 0; i < n; i++) {
        cin >> red[i].x >> red[i].y;
    }
    for (int i = 0; i < n; i++) {
        cin >> blue[i].x >> blue[i].y;
    }
    vi r(n), b(n);
    for (int i = 0; i < n; i++) {
        r[i] = b[i] = i;
    }
    dfs(r, b);
    return 0;
}