[SDOI2014]重建

传送门


这题要用到的是矩阵-树定理。

对于有边权的图,其所有生成树边权积的和等于基尔霍夫矩阵的任意余子式(M(i,i))的行列式。

其中,基尔霍夫矩阵(C=D-A).(D)为发出边权和矩阵,对于任意(u eq v),有(D(u,v) = 0).(A)是邻接矩阵。

总而言之,矩阵-树定理能求出的是(sumlimits_{T} prodlimits_{e in T} p_e).


回头看这道题,这道题让我们求的是(sumlimits_{T} [prodlimits_{e in T} p_e prodlimits_{e otin T}(1-p_e)]).
所以我们对上式变形,化成能求的形式:

[egin{align*} sumlimits_{T} [prodlimits_{e in T} p_e prodlimits_{e otin T}(1-p_e)] &= sumlimits_{T} [prodlimits_{e in T} p_e frac{prodlimits_{e in E} (1-p_e)}{prodlimits_{e in T} (1-p_e)}]\ &= prodlimits_{e in E} (1-p_e) * sumlimits_{T} prodlimits_{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
", sum * Gauss(p, n - 1));
	return 0;
}
原文地址:https://www.cnblogs.com/mrclr/p/15115950.html