[SDOI2014]重建(矩阵树定理)

[SDOI2014]重建(矩阵树定理)

题面

T 国有 N 个城市，用若干双向道路连接。一对城市之间至多存在一条道路。
在一次洪水之后，一些道路受损无法通行。虽然已经有人开始调查道路的损毁情况，但直到现在几乎没有消息传回。
幸运的是，此前 T 国政府调查过每条道路的强度，现在他们希望只利用这些信息估计灾情。
具体地，给定每条道路在洪水后仍能通行的概率，请计算仍能通行的道路恰有 N-1条，且能联通所有城市的概率。

分析

设图为\((V,E)\),满足条件的生成树为\(T\),每条边连通的概率为\(P_i\)则要求的答案是\(\sum_{T} \prod_{e \in T} P_e \prod_{e \notin T} (1-P_e)\).
考虑矩阵树定理实际上求的是\(\sum_{T} \prod_{e \in T} 1\).(把边权看成连乘,一棵树乘起来就是1),我们考虑对矩阵进行变换。

\[\begin{aligned} \sum_{T} \prod_{e \in T} P_e \prod_{e \notin T} (1-P_e)&=\sum_{T} \prod_{e \in T} P_e \frac{\prod_{e \in E}(1-P_e)}{\prod_{e \in T} (1-P_e)} \\ &= \prod_{e \in E}(1-P_e) \sum_{T} \prod_{e \in T} \frac{P_e}{1-P_e}\end{aligned} \]

\(\prod_{e \in E}(1-P_e)\)在输入时就可以求出。考虑如何求后面部分.考虑原来定义:

\[\bm{L}_{i,j}=\begin{cases} \sum_{j=1}^{|V|}[(i,j)\in E],i=j \ (\text{即}i\text{的度})\\ -\sum_{e \in E}[(i,j)=e],i \neq j \ (\text{即边}i,j\text{的个数})\end{cases} \]

不妨把\(1\)换成\(\frac{P_e}{1-P_e}\),记\(i,j\)间的边的概率为\(P_{i,j}\),定义

\[\bm{L}_{i,j}=\begin{cases} \sum_{j=1}^{|V|}\frac{P_{i,j}}{1-P_{i,j}} ,i=j \ (\text{即}i\text{连出去的所有边边权之和 })\\ -\frac{P_{i,j}}{1-P_{i,j}},i \neq j \ (\text{即边}i,j\text{的值,题目保证无重边})\end{cases} \]

那么所求就是\(\det(\bm{L}_0)\)。注意当\(P=1\)时会出现除0错误,把\(P\)减去一个很小的值\(\varepsilon\)即可。

代码

#include<iostream>
#include<cstdio>
#include<cstring> 
#define maxn 50
#define eps 1e-9
using namespace std;
int n;
double g[maxn+5][maxn+5];
double a[maxn+5][maxn+5];
double gauss(int n){
	double ans=1;
//	int sign=1; 
	for(int i=1;i<=n;i++){
		for(int j=i+1;j<=n;j++){
			double d=a[j][i]/a[i][i]; 
			for(int k=1;k<=n;k++) a[j][k]-=d*a[i][k];
		}
		if(a[i][i]==0) return 0;
		ans*=a[i][i];
	}
	return ans;
}
int main(){
	double s=1;
	scanf("%d",&n);
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++){
			scanf("%lf",&g[i][j]);
			if(i==j) continue;//i!=j时是(i,j)的边权和 
			if(g[i][j]+eps>=1) g[i][j]-=eps;//防止g[i][j]=1导致分母为0 
			a[i][j]=-g[i][j]/(1-g[i][j]);
			
			if(i<j) s=s*(1-g[i][j]);
		}
	}
	for(int i=1;i<=n;i++){
		for(int j=0;j<=n;j++) if(i!=j) a[i][i]+=g[i][j]/(1-g[i][j]);//i=j时定义为g[i][j]/(1-g[i][j]),即i点连出去的所有边边权之和 
	} 
	printf("%.10lf\n",gauss(n-1)*s); 
}