题目
T国有N个城市,用若干双向道路连接。一对城市之间至多存在一条道路。
在一次洪水之后,一些道路受损无法通行。虽然已经有人开始调查道路的损毁情况,但直到现在几乎没有消息传回。
辛运的是,此前T国政府调查过每条道路的强度,现在他们希望只利用这些信息估计灾情。具体地,给定每条道路在洪水后仍能通行的概率,请计算仍能通行的道路恰有N-1条,且能联通所有城市的概率。
输入格式
输入的第一行包含整数N。
接下来N行,每行N个实数,第i+l行,列的数G[i][j]表示城市i与j之
间仍有道路联通的概率。
输入保证G[i][j]=G[j][i],且G[i][j]=0;G[i][j]至多包含两位小数。
输出格式
输出一个任意位数的实数表示答案。
你的答案与标准答案相对误差不超过10^(-4)即视为正确。
输入样例
3
0 0.5 0.5
0.5 0 0.5
0.5 0.5 0
输出样例
0.375
提示
1 < N < =50
数据保证答案非零时,答案不小于10^-4
题解
矩阵树定理:
一个图的生成树个数等于矩阵G的\(n - 1\)阶行列式的值
其中矩阵G = D - A
其中D为度数矩阵,只有\(i == j\)的地方不为0,为\(i\)的度数
其中A为邻接矩阵
由基尔霍夫定理,我们直接求\(n - 1\)阶行列式实际上得到的是这个东西:
\[\sum\limits_{Tree} \prod\limits_{e \in Tree} w_e\]
即所有可能的生成树中各边权的积之和
也就是说,我们直接求是这个:
\[\sum\limits_{Tree} \prod\limits_{e \in Tree} P_e\]
但我们要求的答案是这个:
\[\sum\limits_{Tree} \prod\limits_{e \in Tree} P_e \prod\limits_{e \notin Tree} (1 - P_e)\]
我们还差后面那一串,怎么办?
考虑对所有边,令答案乘一个
\[\prod\limits_{e} (1 - P_e)\]
我们就得到了这个:
\[\sum\limits_{Tree} \prod\limits_{e \in Tree} P_e * (1 - P_e) \prod\limits_{e \notin Tree} (1 - P_e)\]
我们只需要把\(\prod\limits_{e \in Tree} P_e * (1 - P_e)\)变成\(\prod\limits_{e \in Tree} P_e\)就可以了
即令
\[\prod\limits_{e \in Tree} w_e = \prod\limits_{e \in Tree} \frac{P_e}{1 - P_e}\]
那么我们令原矩阵A中\(A_{i,j} = \frac{P_{i,j}}{1 - P_{i,j}}\)
问题就圆满解决了
考虑除法问题:
如果一个值等于\(0\),就令其等于\(eps\)
#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define eps 1e-9
#define LL long long int
#define Redge(u) for (int k = h[u],to; k; k = ed[k].nxt)
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define BUG(s,n) for (int i = 1; i <= (n); i++) cout<<s[i]<<' '; puts("");
using namespace std;
const int maxn = 55,maxm = 100005,INF = 1000000000;
double A[maxn][maxn],G[maxn][maxn],ans = 1;
int n;
void gause(){
for (int i = 1; i < n; i++){
int j = i;
for (int k = i + 1; k <= n; k++)
if (fabs(A[k][i]) > fabs(A[j][i]))
j = k;
if (j != i) for (int k = i; k <= n; k++) swap(A[i][k],A[j][k]);
if (fabs(A[i][i]) < eps){
ans = 0;
return;
}
for (j = i + 1; j <= n; j++){
double t = -A[j][i] / A[i][i];
for (int k = i; k <= n; k++)
A[j][k] += A[i][k] * t;
}
}
}
int main(){
scanf("%d",&n);
REP(i,n) REP(j,n){
scanf("%lf",&G[i][j]);
if (fabs(G[i][j]) < eps) G[i][j] = eps;
if (fabs(1 - G[i][j]) < eps) G[i][j] = 1 - eps;
}
REP(i,n) REP(j,n){
A[i][j] = G[i][j] / (1 - G[i][j]);
if (i < j) ans *= (1 - G[i][j]);
}
REP(i,n){
A[i][i] = 0;
REP(j,n) if (j != i) A[i][i] -= A[i][j];
}
gause();
REP(i,n - 1) ans *= A[i][i];
printf("%.10lf",fabs(ans));
return 0;
}