TopCoder SRM 583 Div1 Hard RandomPaintingOnABoard
問題
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); } };