TopCoder SRM 425 Div1 Hard RoadsOfKingdom

問題

n頂点からなる無向完全グラフで表される都市と道がある。
ここに雨がふって、それぞれの道がroads[i][j] / 8の確率で残り、
それ以外の確率で破壊される。

残った道が、木になっている確率を求めよ。

制約条件

n≦16

方針

接続関係のDPがしたくなる。
けれど、この接続関係は頂点を区別するものなので、
状態数が爆発する。（n=16で10億とかそういうの）

なので上手くやる方法を考えなければならない。
木を作っていくDPでできるか見てみると、
dp[bit]を、bitの集合の頂点が、その範囲では木になるような確率として、
それを求めるDPができないかとなる。

結論から言うと遷移をがんばればできて、
まず、1頂点を、番号の一番小さいものAに固定する。
これをこの木の根と考える。

二番目に小さい番号の頂点Bのほうの部分木を切り離す。
これは、根からBのほうへ向けて、じかに出ている辺できってしまうようにすると、一意に定まり、
なおかつ全ての場合を尽くすことができる。

なので、その辺を切って、なおかつきちんと木になっている確率をメモ化再帰で求めればいい。
きちんと木になる確率は、

Aから、Bのほうの部分木へ、一つだけ辺が出ている
Bの部分木と、Aのそれ以外の部分木（Aは除く）に辺がない
Bの部分木が木になる
Bを除いた部分が木になる

が全て起こる確率。
で、計算量削減のために部分集合列挙をO(3^n)でやり、あと必要な部分を前計算しておけばよい。

Editorialの擬似コードは嘘なので注意。

Aから、Bのほうの部分木へ、一つだけ辺が出ている

が、AからBのほうの部分木へ、Bに辺が出ているになっている。

ソースコード

int n;
double e[16][16], one[16][1 << 16], dis[16][1 << 16];
double dp[1 << 16];

double rec(int bit){
	double &res = dp[bit];
	if(res >= 0) return res;
	if(__builtin_popcount(bit) <= 1) return res = 1;
	
	int a = 0;
	for(; !(bit & 1 << a); a++);
	int b = a + 1;
	for(; !(bit & 1 << b); b++);
	int bit2 = bit ^ 1 << a ^ 1 << b;
	res = 0;
	
	for(int mask = bit2; ; mask = (mask - 1) & bit2){
		
		//B以下の部分木と、全体からAとそれを除いた部分森がdisjointである確率 = p
		double p = 1;
		rep(i, n) if((mask ^ 1 << b) & 1 << i) p *= dis[i][bit ^ mask ^ 1 << a ^ 1 << b];
		
		//p
		//かつAからBへの部分木へは1つしか辺がない
		//かつB以下がきちんと木になっている
		//かつB以下の部分木を除いた部分木がきちんと木になっている
		res += p * one[a][mask ^ 1 << b] * rec(bit ^ mask ^ 1 << b) * rec(mask ^ 1 << b);
		
		if(mask == 0) break;
	}
	return res;
}

class RoadsOfKingdom {
	public:
	double getProbability(vector <string> roads) {
		
		n = roads.size();
		rep(i, n) rep(j, n) e[i][j] = (roads[i][j] - '0') / 8.0;
		
		rep(i, n){
			rep(j, 1 << n) if(!(j & 1 << i)){
				double &sum = one[i][j];
				double prod = 1;
				int zero = 0;
				sum = 0;
				
				//one
				rep(k, n) if(j & 1 << k){
					if(e[i][k] == 1) zero++;
					else prod *= 1 - e[i][k];
				}
				rep(k, n) if(j & 1 << k){
					if(e[i][k] == 1){
						if(zero == 1) sum += prod;
					}
					else{
						if(zero == 0) sum += prod / (1 - e[i][k]) * e[i][k];
					}
				}
				
				//dis
				dis[i][j] = zero ? 0 : prod;
			}
		}
		
		rep(i, 1 << n) dp[i] = -1;
		return rec((1 << n) - 1);
	}
};