TopCoder SRM 611, Division 1, Level 2 : Egalitarianism2

問題文

http://community.topcoder.com/stat?c=problem_statement&pm=13008&rd=15844

問題概要

　平面上に N 個の点がある。これらの点を頂点とし、二点間にそのユークリッド距離に等しい重みの辺を張ったグラフを考える。このグラフの全域木であって、辺重みの標準偏差が最小となるものの標準偏差を返せ。
　標準偏差は、利用された辺の内 i 番目のものの重みを 、辺重みの平均値を x とすると で計算される。

解法

　最適解な全域木に於ける辺重みの平均値 x を予め知っているとすると、辺重みと x との差を辺の優先度決定に使うよう変形した Kruskal 法で答えを計算することができます。実際には x の値は未知なので、x になりそうな平均値をいくつか試して最適解を得られるものを探すことになります。このとき、平均値は実数なので、可能な値を全て試すというようなことはできません。しかし、 個の値を試せばよいということが示せます。
　平均値を仮定するとき、同じ辺の優先度を得られるような値が複数あったととします。このとき、それらの値からは全く同じ全域木が得られるので、一つだけ試せば十分です。辺の優先度が変化するのは、ある二つの辺 i, j について、仮定した平均値を m として と の大小関係が入れ替わる点の周辺です。これは と の中点の周辺なので、 の周辺で着目した二辺の優先度が入れ替わります。従って、有効な i, j 全てについてその辺重みの中点を求めれば、辺重みが入れ替わる点を列挙できます。こうして求まった中点たちを昇順に列挙したとすると、隣り合う点の間の（開）区間内では辺の優先度は変化しません。辺の数は 個なので、区間の数は 個です。こうして求まったそれぞれの区間内から、一つずつ適当な値を取り出して仮の平均値とし、上記のような Kruskal 法を適用して、最小の標準偏差を達成するものを探すことで、答えを得ることができます。

コード

#define REP( i, m, n ) for ( int i = (int)( m ); i < (int)( n ); ++i )
#define FOR( v, c ) for ( auto &v : c )
#define ALL( c ) (c).begin(), (c).end()

#define PB( n ) push_back( n )

const int EPS = 1e-7;

struct Edge
{
	int a, b;
	double c;

	Edge( int a, int b, double c ) : a( a ), b( b ), c( c ) {};

	bool operator < ( const Edge &lhs ) const
	{
		return c < lhs.c;
	}
};

class UnionFind;
// UnionFind( N )
// find( x )
// same( x, y )
// unite( x, y )

class Egalitarianism2
{
public:
	int N;
	vector<Edge> es;

	double minStdev( vector <int> x, vector <int> y )
	{
		{ // initilize
			N = x.size();
			es.clear();	
		}

		REP( i, 0, N )
		{
			REP( j, 0, i )
			{
				es.emplace_back( i, j, sqrt( pow( 1. * x[i] - x[j], 2. ) + pow( 1. * y[i] - y[j], 2. ) ) );
			}
		}
		
		const int M = es.size();

		vector<double> ws;
		REP( i, 0, M )
		{
			REP( j, 0, i )
			{
				ws.PB( ( es[i].c + es[j].c ) / 2 );
			}
		}

		double res = DBL_MAX;
		FOR( w, ws )
		{
			res = min( res, solve( w + EPS ) );
		}
		return res;
	}

	double solve( const double mean )
	{
		sort( ALL( es ), [=]( const Edge &e1, const Edge &e2 ){ return abs( mean - e1.c ) < abs( mean - e2.c ); } );
		vector<double> as;
		UnionFind uf( N );

		FOR( e, es )
		{
			if ( !uf.same( e.a, e.b ) )
			{
				as.PB( e.c );
				uf.unite( e.a, e.b );
			}
		}

		const double b = accumulate( ALL( as ), 0. ) / ( N - 1 );
		double c = 0;
		FOR( a, as )
		{
			c += pow( a - b, 2. );
		}
		return sqrt( c / ( N - 1 ) );
	}
};