WARush

SRMの結果とか、解けた問題のコードを書いていきます

SRM611 Div1 Medium "Egalitarianism2"

問題

TopCoder Statistics - Problem Statement

平面上に都市がN個ある。便宜上これらの都市に0~N-1の番号を付ける。i番の都市は、座標(x[i], y[i])にある。

王様はちょうどN-1の道路を作り、都市を繋げたいと考えていた。それぞれの道路は2つの都市を結ぶものでなくてはならない。全ての道路は直線でなければならない。そのため、道路の長さは2点のユークリッド距離となる。道路は交差したり、重なったりしてもよい。(交差している道路は繋がっているわけではない。そのため、N-1の道路で全ての都市をつないだ時、その構造はツリーとなる)

王様は道路を短くすることには熱心ではない。しかし、道路が短かったり長かったりすると、国民は文句をいうだろう。そのため、王様はN-1の道路の距離の標準偏差をなるべく小さくしようと考えていた。

正確には、与えられた実数のシーケンス(a1, ..., aS)があり、これの標準偏差を求めるには次のようにする。まずb = (a1+...+aS) / Sとする。(すなわち平均である)次に、c = (b-a1)^2 + ... +(b-aS)^2とする。最後に、sqrt(c/S)とすれば、標準偏差が求まる。道路はN-1個作るため、SはN-1となる。

あなたはN個の都市の座標を表すint[] x, yが与えられる。標準偏差が最小になるようにN-1個の道路を作ったときの、その標準偏差を返せ。

制約

3 <= N <= 20
-10^6 <= x[i], y[i] <= 10^6


Editorialsを見て

平均をmと決め打ちしちゃおう。それぞれの道路の長さをwiとし、N-1個の道路を作ったときの標準偏差

sqrt( sum( m - wi ) / (N - 1) )

となる。
ここで、最小化したいのは

sum( m - wi )

の部分となる。

最小化するには、それぞれの道路のコストをm - wiとして、MST(最小全域木)のアルゴリズムを使用すればよい。

だが、いったいどれだけの仮平均を試せばよいのだろうか?

MSTは、小さいコストの辺を貪欲的に選んでいくアルゴリズムである。
例えば、道路の長さを小さい順に並べたら下図のような感じだったとする。
f:id:Ekaing:20140315225619j:plain

仮平均mを決めると・・
f:id:Ekaing:20140315225700j:plain

m-wiは次のようになり、MSTにとっての優先度が決まる
f:id:Ekaing:20140315225806j:plain

仮平均mを変えると、優先度もいろいろ変わってくる。
f:id:Ekaing:20140315225836j:plain
f:id:Ekaing:20140315225838j:plain

ここで、mを左右に動かして優先度の全てのセットを試すことができれば、その中のひとつに標準偏差が最小のものとなる優先度のものが必ず含まれることになる。

優先度の状態が変わるのはどのような瞬間かというと、2つの要素の平均点をmが超えたときに優先度の変動が起こる。
f:id:Ekaing:20140315230319j:plainf:id:Ekaing:20140315230323j:plainf:id:Ekaing:20140315230326j:plain

よって、2つの道路の平均の全ての場所を仮平均として試せばよい。ただし、平均そのままをmとすると、2つの道路の優先度が全く一緒になって何か嫌なので、平均+めっちゃ小さい数をmとすればよい。

計算量は、平均点がO(N^4)、MSTのアルゴリズムがO(N^2)となり、O(N^6)となる。

ソースコード

int n;
// i -> j の辺の距離
long long len2[25][25]; 
double len[25][25];

class Egalitarianism2 {

    public:

    double minStdev(vector <int> x, vector <int> y) {

        n = x.size();

        // 各辺の長さを出す
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                len2[i][j] = (long long)(x[i] - x[j]) * (x[i] - x[j])
                    + (long long)(y[i] - y[j]) * (y[i] - y[j]);
                len[i][j] = sqrt(len2[i][j]);
            }
        }

        // 仮平均の候補を出す
        vector <double> cand;
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                for (int k = 0; k < n; k++) {
                    for (int l = 0; l < n; l++) {
                        cand.push_back((len[i][j] + len[k][l]) / 2 + 1e-8);
                    }
                }
            }
        }

        double res = 1e18;
        
        for (int i = 0; i < cand.size(); i++) {
            double m = cand[i]; // 仮平均
            long long sum2;
            double sum;

            // プリム法
            prim(m, sum2, sum);

            // 標準偏差算出
            double ave = sum / (n - 1);
            double ans = sqrt((sum2 - 2 * ave * sum + (n - 1) * ave * ave) / (n - 1));
            res = min(res, ans);
        }

        return res;
    }

    // プリム法
    void prim(double m, long long &sum2, double& sum) {

        sum2 = 0LL;
        sum = 0.0;

        bool used[25];
        memset(used, false, sizeof(used));

        double mincost[25];
        int before[25];
        for (int i = 0; i < n; i++) {
            mincost[i] = abs(len[0][i] - m);
            before[i] = 0;
        }
        used[0] = true;

        for (int i = 1; i < n; i++) {
            int add = -1;
            double w = 1e18;
            for (int v = 0; v < n; v++) {
                if (!used[v] && mincost[v] < w) {
                    add = v;
                    w = mincost[v];
                }
            }

            used[add] = true;
            sum2 += len2[before[add]][add];
            sum += len[before[add]][add];

            for (int u = 0; u < n; u++) {
                if (abs(len[add][u] - m) < mincost[u]) {
                    mincost[u] = min(mincost[u], abs(len[add][u] - m));
                    before[u] = add;
                }
            }
        }
    }
};