[SDOI2014]重建

Posted by Panda2134's Blog on July 6, 2018

题意

给出一个无向图,每个边有个存在概率 $p_{<u, v>}$ ,求所有存在边恰好组成生成树的概率。

思路

第一眼感觉是裸的矩阵树定理,后来发现自己 $\text{Naive}$!

滚去膜拜 dkw 聚聚的题解,学会了一种矩阵树定理的化简技巧。

我们显然可以知道答案是下式:

\[\sum_{T}\left(\prod_{<u, v> \in T} p_{u, v} \times \prod_{<u, v> \notin T} (1-p_{<u, v>})\right)\]

而矩阵树定理求的是这个式子:

\[\sum_T\left(\prod_{<u, v> \in T} A_{<u, v>}\right)\]

考虑如何变形。上面那个式子里面有 2 个 $\prod$,所以我们考虑在下面这个式子外面乘上一个系数,再在化简的时候分开。这样就可以拆出 2 个 $\prod$.

列个方程:

\[\sum_{T}\left(\prod_{<u, v> \in T} p_{u, v} \times \prod_{<u, v> \notin T} (1-p_{<u, v>})\right) = \left(\prod_{<u, v>} (1-p_{<u, v>})\right) \times \sum_T\left(\prod_{<u, v> \in T} A_{<u, v>}\right)\]

于是化简后有

\[\sum_{T}\left(\prod_{<u, v> \in T} p_{u, v} \times \prod_{<u, v> \notin T} (1-p_{<u, v>})\right) =\sum_T\left(\prod_{<u, v> \in T} (1-p_{<u, v>})A_{<u, v>}\times \prod_{<u, v> \notin T} (1-p_{<u, v>})\right)\]

那么显然

\[p_{<u, v>} = (1-p_{<u, v>}) A_{<u, v>}\]

系数化为1:

\[A_{<u, v>} = \frac{p_{<u, v>}}{1-p_{<u, v>}}\]

现在是货真价实的矩阵树定理裸题了。用矩阵树定理算 $\sum_T \left(\prod_{<u, v> \in T} A_{<u, v>}\right)$ ,然后乘上系数即可(没乘系数调试了好久……以后记住,数学相关题目,一定要写清楚待求的式子!)。

$p = 1$ 的时候分母为 0 ,怎么办?考虑到 $\mathtt{nan}$ 不能参与运算,设其值为 $\frac{1}{\epsilon}$ 即可。

(其实精度相关处理是最恶心的地方……)

代码

#include <bits/stdc++.h>
using namespace std;

typedef long double ld;

const double eps = 1e-6;
const int MAXN = 50;

int n;
ld mul = 1.0l, A[MAXN + 10][MAXN + 10], G[MAXN + 10][MAXN + 10];

int dcmp(double x) {
	return fabs(x) < eps ? 0 : (x > 0 ? 1 : -1);
}

int main() {
	scanf("%d", &n);
	for(int i = 1; i <= n; i++)
		for(int j = 1; j <= n; j++) {
			scanf("%Lf", &G[i][j]);
			if(dcmp(G[i][j]) == 0) G[i][j] = eps;
			if(dcmp(G[i][j]-1) == 0) G[i][j] = 1-eps;
			if(i < j) mul *= (1 - G[i][j]);
		}
	
	for(int i = 1; i <= n; i++) {
		for(int j = 1; j <= n; j++)
			if(i != j) {
				if(dcmp(G[i][j]-1) != 0)
					A[i][j] -= G[i][j] / (1-G[i][j]);
				else
					A[i][j] -= 1 / eps;
				A[i][i] -= A[i][j];
			}
	}
	
	n--;
	for(int i = 1; i <= n; i++) {
		int k = i;
		for(int j = n; j >= i; j--)
			if(dcmp(A[j][i]) != 0) k = j;
		
		if(k != i) {
			for(int j = 1; j <= n; j++)
				swap(A[i][j], A[k][j]);
			mul *= -1;
		}
		
		if(dcmp(A[i][i]) == 0)
			break;
		
		for(int j = i+1; j <= n; j++) {
			ld t = A[j][i] / A[i][i];
			for(int k = 1; k <= n; k++)
				A[j][k] -= t * A[i][k];
		}
	}
	
	for(int i = 1; i <= n; i++)
		mul = mul * A[i][i];
	
	printf("%Lf", mul);
	
	return 0;
}