AtCoder Regular Contest 021 D - だいたい最小全域木

問題

arc021.contest.atcoder.jp

解法

5000*5000/2 のすべての辺を貼ってから最小全域木のアルゴリズムを適用しても間に合いませんが, 最小全域木に使われる可能性の高い辺のみに注目して辺を貼れば, 最小全域木に近いものを作ることが出来ます。

今回の場合, コストの作り方から２つの頂点 a, b について, a[i] と b[i] の符号がなるべく一致しているのが有利です(符号がバラバラだと類似度が自然に小さくなるので, コストは高くなる)。また, 問題文に書いてあるように, 頂点の座標は, すべての値が一様に出ると言われているので, 符号はプラスマイナスがそれぞれ 1/2 の確率で現れるとみなすことが出来ます。

ということで, bitset で各頂点の符号の分布を管理して, 例えば「符号が 110 個以上一致していない頂点については辺がないものとする」というようにすれば, ほとんど最小全域木を作ることが出来ます。

a と b の符号のそれぞれが一致する確率は p = 1/2 の二項分布で, その分散は試行回数 n に対して p * (1-p) * n = n/4 = 50 と書けます。で, なんかの定理で確率分布は正規分布に収束するみたいなことを言っていたような気がするので, 1σ区間を取った, 107 個以上一致してない頂点は (100-68)/2 = 16% ぐらいであることが言えます。かなりせめて 2σ区間をとっても AC しました。

struct UnionFind {
    vector<int> par;
    int n, cnt;
    UnionFind(const int& x = 0) {init(x);}
    void init(const int& x) {par.assign(cnt=n=x, -1);}
    inline int find(const int& x) {return par[x] < 0 ? x : par[x] = find(par[x]);}
    inline bool same(const int& x, const int& y) {return find(x) == find(y);}
    inline bool unite(int x, int y) {
        if ((x = find(x)) == (y = find(y))) return false;
        --cnt;
        if (par[x] > par[y]) swap(x, y);
        par[x] += par[y];
        par[y] = x;
        return true;
    }
    inline int count() const {return cnt;}
    inline int count(int x) {return -par[find(x)];}
};

ll pos[5000][200];
double N[5000];
bitset<200> bits[5000];
struct edge {
    int u, v;
    double cost;
    edge() {}
    edge(int u, int v, double cost) : u(u), v(v), cost(cost) {}
    bool operator<(const edge& rhs) const {return cost < rhs.cost;}
};

vector<edge> G;

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    uint x = 123456789, y = 362436069, z = 521288629, w;
    cin >> w;
    for (int i = 0; i < 5000; i++) for (int j = 0; j < 200; j++) {
        uint t = x^(x<<11);
        x = y;
        y = z;
        z = w;
        w = (w^(w>>19))^(t^(t>>8));
        int v = w%100000 - 50000;
        if (v >= 0) v = v+1;
        pos[i][j] = v;
    }
    for (int i = 0; i < 5000; i++) {
        for (int j = 0; j < 200; j++) {
            N[i] += pos[i][j]*pos[i][j];
            if (pos[i][j] < 0) bits[i].set(j);
        }
        N[i] = sqrt(N[i]);
    }
    double sum = 0;
    for (int i = 0; i < 5000; i++) for (int j = i+1; j < 5000; j++) {
        if ((bits[i]^bits[j]).count() > 90) continue;
        double cost = 0;
        for (int k = 0; k < 200; k++) {
            cost += pos[i][k] * pos[j][k];
        }
        cost /= N[i]*N[j];
        cost = 1-cost;
        G.emplace_back(i, j, cost);
    }
    sort(G.begin(), G.end());
    UnionFind uf(5000);
    vector<pii> ans;
    for (edge e : G) {
        if (uf.same(e.u, e.v)) continue;
        uf.unite(e.u, e.v);
        ans.emplace_back(e.u, e.v);
        sum += e.cost;
    }
    //cout << ans.size() << endl;
    //cout << sum << endl;
    for (pii p : ans) cout << p.first+1 << " " << p.second+1 << endl;
    return 0;
}