AtCoder ARC #021 : D - だいたい最小全域木

本番には解けなかったので、後日他人の解説をみて解いた。
http://arc021.contest.atcoder.jp/tasks/arc021_4

問題

M=200次元の座標系において、N=5000個の点が(原点以外で)おおむねランダムで配置されていたとする。
ここでN点から最小全域木を作りたい。

２点間の距離は、原点から各点に伸びるベクトルを正規化したのち、２ベクトルの内積cを求めた場合、1-cとする。
最小全域木の1.01倍以下の総辺長になるグラフを返せ。

解法

O(N^3)かかる一般的な最小全域木は構築できない。
よって何らかの方法で探索空間を削減する必要がある。

本番自分が試したダメ解法は、まずN点の重心を求め、そこに近い順に処理済みの中の最寄の点に連結することだった。
これはO(MN^2)かかって若干TLEするうえ、最小全域木の1.02倍程度の総辺長となりダメだった。

よって、もっと高速で最適な探索法が必要。
まずは点同士の比較を軽量化する。
事前にN点の各頂点の正負から200bitのbitmaskを作る。
2点間のおおよその類似度を両者のxorより求めることができる。
これは各点ペアにつきO(M)かかるが、200個のベクトル要素の内積を取るのに比べると、200bitのxorは64bit値のxorを4回取るだけで高速である。
よって同じO(MN^2)でも定数項がかなり小さい。

以下の例では、符号の正負が110個以上等しい点同士を類似であるとみなすことができる。
この110を調整することで、速度と精度のトレードオフを調整できる。
後は類似な点同士間での辺のみを対象としてクラスカル法を動かせばよい。

なお、正負が等しい点が大量にあったらどうなるの…と心配になるかもしれないが、各点はランダムに生成されるため、極端な偏りはないと仮定できる。

ll P[5001][201];
double S[5001][201];
double R[5001][5001];
double BB[5001];
ll MB[5001][4];

class UF {
	public:
	static const int ufmax=5002;
	int ufpar[ufmax],ufrank[ufmax];
	UF() { init();}
	void init(){int i; FOR(i,ufmax) { ufpar[i]=i; ufrank[i]=0; } }
	int find(int x) {	return (ufpar[x]==x)?(x):(ufpar[x] = find(ufpar[x]));}
	int operator[](int x) {return find(x);}
	void unite(int x,int y) {
		x = find(x); y = find(y);
		if(x==y) return;
		if(ufrank[x]<ufrank[y]) ufpar[x]=y;
		else {ufpar[y]=x; if(ufrank[x]==ufrank[y]) ufrank[x]++;}
	}
};

void solve() {
	int f,i,j,k,l;
	unsigned int x,y,t;
	unsigned int w,z;
	cin>>w;
	
	x = 123456789;
	y = 362436069;
	z = 521288629;
	FOR(i,5000) {
		FOR(j,200) {
			t=x^(x<<11);
			x=y;
			y=z;
			z=w;
			w=(w ^ (w >> 19)) ^ (t ^ (t >> 8));
			ll v = w % 100000;
			v -= 50000;
			if(v>=0) v++;
			P[i][j]=v;
			if(P[i][j]>0) MB[i][j/50] |= 1LL<<(j%50);
			BB[i]+=v*v;
		}
		BB[i]=sqrt(BB[i]);
		FOR(j,200) S[i][j]=P[i][j]/BB[i];
	}
	
	set<pair<double,pair<int,int> > > V;
	FOR(x,5000) for(y=x+1;y<5000;y++) {
		int diff=0;
		diff += __builtin_popcountll(MB[x][0]^MB[y][0]);
		diff += __builtin_popcountll(MB[x][1]^MB[y][1]);
		diff += __builtin_popcountll(MB[x][2]^MB[y][2]);
		diff += __builtin_popcountll(MB[x][3]^MB[y][3]);
		
		if(diff<90) {
			double co=0;
			FOR(i,200) co+=S[x][i]*S[y][i];
			V.insert(make_pair(1-co,make_pair(x,y)));
		}
	}
	
	UF uf;
	ITR(it,V) {
		x=it->second.first;
		y=it->second.second;
		
		if(uf[x]!=uf[y]) {
			_P("%d %d\n",x+1,y+1);
			uf.unite(x,y);
		}
	}
	
}

まとめ

本番も何かしら探索空間を狭めようとしたが、正負による分類までは思いついてもbitmask化することは思いつかなかった…。