TopCoder SRM 583 Div1 Hard RandomPaintingOnABoard

TopCoder 包除原理 動的計画法 確率・場合の数

問題

h行w列の行列があって、(i, j)にはp[i][j]の数字が書かれている。
この行列の一つのセルを、p[i][j]に比例する確率で選ぶことを繰り返す。


全ての行と列について、その行、列に属するセルが一度以上選ばれる状態になったら、
操作を終了する。
操作が終了するまでの回数の期待値を求めよ。

制約条件

h, w≦21
h*w≦150
0≦p[i][j]≦9
各行、列には少なくとも一つの0でない数字がある

方針

editorial見た。
愚直なDPでやろうとするとO(2^(w+h) * poly(w,h))とかかかってしまって無理。


こんな式変形ができる。
E[i] = Σi * (i回目で終了する確率）
= Σi * {(i-1回で終了しない確率) - (i回で終了しない確率)}
= Σi * (i-1回で終了しない確率) - Σi * (i回で終了しない確率)
= Σi * (i-1回で終了しない確率) - Σ(i-1) * (i-1回で終了しない確率)
= Σ(i回で終了しない確率)


なので、i回で終了しない確率を0から無限まで足せばよい。
各iについては包除原理を使うと計算できる。
ある行、列の集合をSとする。
s∈Sについて、sがi回選ばれない確率は、一回でsが選ばれる確率をP(s)とすれば、
(1 - P(s))^iとなる。


i回後に終了しない確率は、包除原理から、
Σ{(1 - P(s))^1 | |s| = 1}
- Σ{(1 - P(s))^1 | |s| = 2}
+ Σ{(1 - P(s))^1 | |s| = 3}……
と計算できるから、Σ{(1 - P(s))^i | s is odd} - Σ{(1 - P(s))^i | s is even}となる。
P(s)はΣ{p[i][j] | (i, j)∈s} / Σp[i][j]である。


ので、Σ{p[i][j]}の値ごとに、奇数・偶数にわけて場合の数を数えればよさそうである。
列のほうで使う成分を固定すれば、後はナップサックDPになって数えられる。


editorialのところではもう少し上手くやっていて、
偶数の場合の数と奇数の場合の数の差分だけを持てばよいので、差分だけでやっている。
このほうが桁落ちが発生しないし定数倍高速。

ソースコード

double dp[2][2][1500];

class RandomPaintingOnABoard {
	public:
	double expectedSteps(vector <string> prob) {
		
		int h = prob.size();
		int w = prob[0].size();
		
		if(w > h){
			vector<string> v(w, string(h, ' '));
			
			rep(i, h) rep(j, w) v[j][i] = prob[i][j];
			swap(h, w);
			swap(prob, v);
		}
		
		int p[22][22];
		int sum = 0;
		double ans = 0;
		
		rep(i, h) rep(j, w){
			p[i][j] = prob[i][j] - '0';
			sum += p[i][j];
		}
		
		rep(mask, 1 << w){
			
			memset(dp, 0, sizeof(dp));
			int parity = __builtin_parity(mask);
			int cur = 0, next = 1;
			
			/*これをやれば簡単だし桁落ちが発生しにくい…
			dp[0] = parity ? -1 : 1;
			
			rep(i, h){
				int s = 0;
				rep(j, w) if(mask & 1 << j) s += p[i][j];
				
				for(int j = sum; j >= 0; j--){
					
					dp[j] *= -1;
					if(j + s <= sum) dp[j + s] -= dp[j];
				}
			}
			rep(i, sum) ans += dp[i] * sum / (sum - i);
			*/
			
			
			dp[0][parity][0] = 1;
			
			rep(i, h){
				int s = 0;
				rep(j, w) if(mask & 1 << j) s += p[i][j];
				
				rep(j, sum+1) rep(k, 2) dp[next][k][j] = dp[cur][k][j];
				
				for(int j = sum - s; j >= 0; j--) rep(k, 2){
					
					dp[next][k^1][j + s] += dp[cur][k][j];
				}
				swap(cur, next);
			}
			rep(i, sum){
				ans += dp[cur][0][i] * sum / (sum - i);
				ans -= dp[cur][1][i] * sum / (sum - i);
			}
		}
		return abs(ans);
	}
};