[SDOI2014]重建(矩阵树定理)
题面
T 国有 N 个城市，用若干双向道路连接。一对城市之间至多存在一条道路。
在一次洪水之后，一些道路受损无法通行。虽然已经有人开始调查道路的损毁情况，但直到现在几乎没有消息传回。
幸运的是，此前 T 国政府调查过每条道路的强度，现在他们希望只利用这些信息估计灾情。
具体地，给定每条道路在洪水后仍能通行的概率，请计算仍能通行的道路恰有 N-1条，且能联通所有城市的概率。
分析
设图为((V,E)),满足条件的生成树为(T),每条边连通的概率为(P_i)则要求的答案是(sum_{T} prod_{e in T} P_e prod_{e
otin T} (1-P_e)).
考虑矩阵树定理实际上求的是(sum_{T} prod_{e in T} 1).(把边权看成连乘,一棵树乘起来就是1),我们考虑对矩阵进行变换。
[egin{aligned} sum_{T} prod_{e in T} P_e prod_{e
otin 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))在输入时就可以求出。考虑如何求后面部分.考虑原来定义:
[m{L}_{i,j}=egin{cases} sum_{j=1}^{|V|}[(i,j)in E],i=j ( ext{即}i ext{的度})\ -sum_{e in E}[(i,j)=e],i
eq j ( ext{即边}i,j ext{的个数})end{cases}
]
不妨把(1)换成(frac{P_e}{1-P_e}),记(i,j)间的边的概率为(P_{i,j}),定义
[m{L}_{i,j}=egin{cases} sum_{j=1}^{|V|}frac{P_{i,j}}{1-P_{i,j}} ,i=j ( ext{即}i ext{连出去的所有边边权之和 })\ -frac{P_{i,j}}{1-P_{i,j}},i
eq j ( ext{即边}i,j ext{的值,题目保证无重边})end{cases}
]
那么所求就是(det(m{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
",gauss(n-1)*s);
}