[SDOI2014]重建

這題要用到的是矩陣-樹定理。

對於有邊權的圖，其全部生成樹邊權積的和等於基爾霍夫矩陣的任意餘子式\(M(i,i)\)的行列式。

其中，基爾霍夫矩陣\(C=D-A\).\(D\)爲發出邊權和矩陣，對於任意\(u \neq v\)，有\(D(u,v) = 0\).\(A\)是鄰接矩陣。

總而言之，矩陣-樹定理能求出的是\(\sum\limits_{T} \prod\limits_{e \in T} p_e\).


回頭看這道題，這道題讓咱們求的是\(\sum\limits_{T} [\prod\limits_{e \in T} p_e \prod\limits_{e \notin T}(1-p_e)]\).
因此咱們對上式變形，化成能求的形式：

 

\[\begin{align*} \sum\limits_{T} [\prod\limits_{e \in T} p_e \prod\limits_{e \notin T}(1-p_e)] &= \sum\limits_{T} [\prod\limits_{e \in T} p_e \frac{\prod\limits_{e \in E} (1-p_e)}{\prod\limits_{e \in T} (1-p_e)}]\\ &= \prod\limits_{e \in E} (1-p_e) * \sum\limits_{T} \prod\limits_{e \in T} \frac{p_e}{1-p_e} \end{align*}\]

 

所以咱們將每一條邊的邊權改爲\(\frac{p}{1-p}\)，再構造基爾霍夫矩陣求解餘子式的行列式就好了。


可是若是\(p=1\)或\(0\)怎麼辦？

這個我想了好一下子都沒出來，只能看題解了：由於題目說精度在\(10^{-4}\)就算對，因此對於上述狀況，咱們能夠令\(p=1-eps\)和\(eps\)就行的通了！這個着實妙。

#include<cstdio>
#include<iostream>
#include<cmath>
#include<algorithm>
#include<cstring>
using namespace std;
#define In inline
typedef double db;
const db eps = 1e-8;
const int maxn = 55;

int n;
db p[maxn][maxn];

In db Gauss(db f[][maxn], int n)
{
	db ret = 1;
	for(int i = 1; i <= n; ++i)
	{
		int pos = i;
		for(int j = i + 1; j <= n; ++j)
			if(fabs(f[j][i]) > fabs(f[pos][i])) pos = j;
		if(pos != i) swap(f[pos], f[i]), ret = -ret;
		if(fabs(f[i][i]) < eps) return 0;
		for(int j = i + 1; j <= n; ++j)
		{
			db tp = f[j][i] / f[i][i];
			for(int k = i; k <= n; ++k) f[j][k] -= tp * f[i][k];
		}
		ret *= f[i][i];
	}
	return ret;
}

int main()
{
	scanf("%d", &n);
	for(int i = 1; i <= n; ++i)
		for(int j = 1; j <= n; ++j) scanf("%lf", &p[i][j]); 
	db sum = 1;
	for(int i = 1; i <= n; ++i)
		for(int j = 1; j <= n; ++j)
		{
			if(i == j) continue;
			if(1 - p[i][j] < eps) p[i][j] = 1 - eps;
			else if(p[i][j] < eps) p[i][j] = eps;
			if(j > i) sum *= (1 - p[i][j]);
			p[i][j] = -p[i][j] / (1 - p[i][j]);
		}
	for(int i = 1; i <= n; ++i)	
		for(int j = 1; j <= n; ++j) if(i ^ j) p[i][i] -= p[i][j];
	printf("%.9lf\n", sum * Gauss(p, n - 1));
	return 0;
}