BZOJ3534 [Sdoi2014]重建 【矩陣樹定理】
題目
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;
}BZOJ3534 [Sdoi2014]重建 【矩陣樹定理】