[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;
}