cafelier のSRM参加記録です。コンテスト中に考えてたことを執拗に全部書き残すとどうなるだろうかという試み本番中にこういうコードが書きたかったなあ、という後で書いた反省コードを書き残す試み
スパムが来たのでしばらくコメント欄をはてなユーザ限定にしています、すみません、
SRM | |
SRM425 の成績・ソース (要ログイン) : WA/AC/-- : ひどいポカをした
あとで | |
425の1000点問題を考える。
一昨日くらいにようやくd:id:tsukunoの言ってたDPというのを理解したので、あんど、O(3^(N-1) * N) にはなるんじゃないだろうかと思ったので組んでみてました。最悪ケース（16頂点、確率0%や100%の辺が一個もない）で3秒ちょい。少し足りない…
#include <vector> #include <string> #include <map> using namespace std; typedef unsigned int uint; static const uint OUT=0, IN=1, ROOT=1, LEAF=2, HAVELEAF=0xAAAAAAAA; #define is(mm,j,t) (((mm>>(j*2))&3)==(t)) #define isnot(mm,j,t) (((mm>>(j*2))&3)!=(t)) class RoadsOfKingdom { public: map<uint,double> memo; double P[16][16]; double getProbability(vector <string> roads) { // 入力 int N = roads.size(); for(int i=0; i<N; ++i) for(int j=0; j<N; ++j) P[i][j] = (roads[i][j]-'0')/8.0; // メモ化再帰 double ans = 0.0; for(uint m=0; m<(1<<N); m+=2) { // すべての全域木mmに関して... uint mm = ROOT; for(int i=1; i<N; ++i) mm |= ((m>>i)&1 ? IN : LEAF) << i*2; ans += rec(mm, N); // その全域木ができる確率 } return ans; } // 部分木mmを、元グラフの各ノードを // (0: mmに入ってない, 1: mmの内部ノード, 2: mmのリーフ) // というビットパターンで表す。16頂点なので32bitで足りる。 double rec(uint mm, int N) { // trivial cases if( mm == ROOT ) return 1.0; if( isnot(mm,0,ROOT) || !(mm&HAVELEAF) ) return 0.0; // memoized? map<uint,double>::iterator it = memo.find(mm); if( it != memo.end() ) return it->second; // 一番頂点番号の小さいリーフ f をとりのぞく int f = 0; for(; f<N; ++f) if( is(mm,f,LEAF) ) break; // P[f][*] が全部ぶっ壊れる確率 double alldead = 1.0; int numAlive = 0; for(int j=0; j<N; ++j) if( isnot(mm,j,OUT) ) if( P[f][j] == 1.0 ) ++numAlive; else alldead *= 1.0 - P[f][j]; if( numAlive >= 2 ) return 0.0; // 「fの親がjだった場合」をmmの全ての内部ノードjについて調べて足し算 double px = rec( mm^(LEAF<<f*2), N ); // jがf以外にも子供を持ってる場合の確率 double ans = 0.0; for(int j=0; j<N; ++j) if( is(mm,j,IN) && P[f][j] != 0.0 ) { double pdead = alldead; if( P[f][j] != 1.0 ) if( numAlive ) continue; else pdead /= 1.0 - P[f][j]; double prest = rec( mm^(LEAF<<f*2)^(3u<<j*2), N ) + px; // fがjの唯一の子の場合 ans += P[f][j] * pdead * prest; } return memo[mm] = ans; } };
1.5倍くらいだと小手先の技でもどうにかなりそうな気がしなくもないのですが、まあそれよりはアルゴリズムの改良でなんとかしたい…。そうこうしているうちにEditorialが出てしまったけれど見ないぞー。
実際は 3^(N-1) もなくて、memoに入る状態数は高々120万弱なんですよね。常に最初のLeafを除く方法で出てくる木構造はそのくらいしかない。なのでそういう状態だけもっと直接的に列挙する方法とかがあれば良さそうなんだけど。うむむ
あとで | |
class RoadsOfKingdom { public: vector< vector<double> > P; double getProbability(vector <string> roads) { P = read_input(roads); memo = vector<double>(1<<P.size(), -1); memo0 = vector<double>(P.size()<<P.size(), -1); memo1 = vector<double>(P.size()<<P.size(), -1); return prob_tree( (1<<P.size())-1 ); } vector< vector<double> > read_input(const vector<string>& R) { int N = R.size(); vector< vector<double> > P(N, vector<double>(N)); for(int i=0; i<N; ++i) for(int j=0; j<N; ++j) P[i][j] = (R[i][j] - '0') / 8.0; return P; } vector<double> memo; double prob_tree(int S) { if( memo[S] >= 0 ) return memo[S]; // the first and the second smallest IDs in S int v, w; for(v=0; (1<<v)<=S; ++v) if( S & 1<<v ) break; for(w=v+1; (1<<w)<=S; ++w) if( S & 1<<w ) break; // if |S| < 2 then it always forms a tree if( (1<<w) > S ) return memo[S] = 1.0; // Let's consider v as the "root node" of S, and try all possible "subtrees" T containing w. // The situation is (other nodes)=v-(T where w in T) double p = 0.0; for(int T=S; T; T=(T-1)&S) if( (T & 1<<w) && !(T & 1<<v) ) p += prob_tree(T) * prob_tree(S&~T) * prob_separated(T, S&~T&~(1<<v)) * prob_one(v, T); return memo[S] = p; } double prob_separated(int S1, int S2) { double p = 1.0; for(int v=0; (1<<v)<=S1; ++v) if( S1 & 1<<v ) p *= prob_zero(v, S2); return p; } vector<double> memo0; double prob_zero(int v, int S) // 0 connection between v and S { const int key = v<<P.size() | S; if( memo0[key] >= 0 ) return memo0[key]; double p = 1.0; for(int u=0; (1<<u)<=S; ++u) if( S & 1<<u ) p *= 1.0 - P[v][u]; return memo0[key] = p; } vector<double> memo1; double prob_one(int v, int S) // exactly 1 connection between v and S { const int key = v<<P.size() | S; if( memo1[key] >= 0 ) return memo1[key]; double p = 0.0; for(int c=0; (1<<c)<=S; ++c) if( S & 1<<c ) { double q = 1.0; for(int u=0; (1<<u)<=S; ++u) if( S & 1<<u ) q *= u==c ? P[v][u] : 1-P[v][u]; p += q; } return memo1[key] = p; } };
presented by
cafelier/k.inaba
under
