torus711 のアレ

主に競技プログラミングの問題について書きます

TopCoder SRM 575, Division 1, Level 2 : TheSwapsDivOne

概要

数の列が与えられ、これに対し以下のような操作をする。

  1. 異なる二つの位置を選び、それらの数を交換する操作を k 回
  2. 空でない連続した部分列を一つ選び、その総和をとる

両操作に於いて位置がランダムに選ばれるとき、得られる部分列の総和の期待値を求めよ。

解法

k 回のスワップ操作のあと、各位置についてその位置にある数字が元のものと異なる( or 同一である)確率は等しくなります。
解法は、この確率を求めるパートと、答えとなる期待値を求めるパートに二分されます。

k 回のスワップ後に異なる( or 同一の)ものが入っている確率を求める

dp[ 操作回数 ][ 異なるかどうか ] := 確率
として動的計画法で求めます。
初期値は dp[0][0] = 1, dp[0][1] = 0 としてこれ以降の状態を逐次計算します。
ある時点で同一のものが入っている確率を p_0 、異なるものが入っている確率を p_1 、数列の項数を N とします。
この時、次の操作を行ったあとのそれぞれの確率 p_0'p_1' は次のようにして計算されます。
 p_0' = \frac{( N - 1 )( N - 2 )}{N( N - 1 )}p_0 + \frac{ 2 }{ N( N - 1 )}p_1 = \frac{ N - 2 }{ N }p_0 + \frac{ 2 }{ N( N - 1 )}p_1
 p_1' = \frac{ 2( N - 1 )}{ N( N - 1 )}p_0 + ( \frac{( N - 1 )( N - 2 )}{ N( N - 1 )} + \frac{ 2( N - 2 )}{ N( N - 1 )})p_1 = \frac{ 2 }{ N }p_0 + ( \frac{ N - 2 }{ N } + \frac{ 2( N - 2 )}{ N( N - 1 )})p_1
分数部分は選び方の条件を満たす選び方を、選び方の総数で割っています。
条件を満たす選び方は数式に出てくる順に、「着目している位置が選ばれない選び方」「着目している位置が選ばれ、他方は元あった数である選び方」「片方が着目している位置であるような選び方」です。

期待値を計算する

期待値は、その事象での値に発生する確率を掛けたものです。
従って、各位置毎にそこが部分列に含まれる確率を計算し、その位置にある数をそれに掛けて足し合わせます。
元の数と異なるものが入っている場合の数は、他の全ての値の平均であると考えます。
ある位置 i ( 0-indexed ) が部分列に含まれる確率 p は次のように計算されます。
 p = 2\frac{( i + 1 )( N - i )}{ N( N + 1 ) }
// なんで分母が N^2 ではなく N( N + 1 ) なのか分かっていません……
ともかく、あとはこの確率に値を掛けたものを足し合わせるだけです。

コード

typedef vector<int> VI;

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

class TheSwapsDivOne
{
public:
	double find( vector <string> sequence, int k )
	{
		VI s;
		{
			string tmp( accumulate( ALL( sequence ), string() ) );
			transform( ALL( tmp ), back_inserter( s ), bind2nd( minus<int>(), '0' ) );
		}

		const int N = s.size();

		vector< vector<double> > dp( k + 1, vector<double>( 2, 0. ) );
		// 0 : 同じものが残っている確率 / 1 : 違うものが入っている確率
		dp[0][0] = 1;
		REP( i, 0, k )
		{
			dp[ i + 1 ][0] = dp[i][0] * ( N - 2. ) / N + dp[i][1] * 2. / ( N * ( N - 1 ) );
			dp[ i + 1 ][1] = dp[i][0] * 2. / N + dp[i][1] * ( ( N - 2. ) / N + ( 2. * ( N - 2 ) ) / ( N * ( N - 1 ) ) );
		}

		int sum = accumulate( ALL( s ), 0 );
		double res = 0.;
		REP( i, 0, N )
		{
			double p = ( ( i + 1. ) * ( N - i ) ) / ( N * ( N + 1 ) ) * 2;
			res += dp[k][0] * p * s[i];
			res += dp[k][1] * p * ( (double)( sum - s[i] ) / ( N - 1 ) );
		}

		return res;
	}
};