BZOJ 3534: [Sdoi2014]重建 矩阵树定理

3534: [Sdoi2014]重建

Time Limit: 10 Sec  Memory Limit: 512 MBSec  Special Judge
Submit: 804  Solved: 366
[Submit][Status][Discuss]

Description

  T国有N个城市，用若干双向道路连接。一对城市之间至多存在一条道路。
    在一次洪水之后，一些道路受损无法通行。虽然已经有人开始调查道路的损毁情况，但直到现在几乎没有消息传回。
    辛运的是，此前T国政府调查过每条道路的强度，现在他们希望只利用这些信息估计灾情。具体地，给定每条道路在洪水后仍能通行的概率，请计算仍能通行的道路恰有N-1条，且能联通所有城市的概率。

Input

  输入的第一行包含整数N。
  接下来N行，每行N个实数，第i+l行，列的数G[i][j]表示城市i与j之
间仍有道路联通的概率。
    输入保证G[i][j]=G[j][i]，且G[i][j]=0；G[i][j]至多包含两位小数。

Output

    输出一个任意位数的实数表示答案。
    你的答案与标准答案相对误差不超过10^(-4)即视为正确。

Sample Input

3
0 0.5 0.5
0.5 0 0.5
0.5 0.5 0

Sample Output

0.375

HINT

1 < N < =50

数据保证答案非零时，答案不小于10^-4

基尔霍夫矩阵的任意一个代数余子式是所有生成树的边权积的和 

（代数余子式定义：在一个n阶 行列式A中,把

   

元

   

所在的第

   

行和第

   

列划去后，留下来的

   

阶 行列式叫做

 

元

   

的余子式，记作

   

;即

 

，

 

叫做

   

元

   

的代数余子式）

注意：其中所求的

   

为一个数值，并非矩阵。

邻接矩阵中的的权可以是其他权值 比如概率

这样计算出来的就是所有生成树的概率和

然而这并没有考虑不属于树上的边

所以生成一颗生成树T的概率应该是

然而纵使有了式子也不是很会求TAT

式子的变换还是蛮精巧的

设G要求的矩阵,P是给出的矩阵

我们令

    

令

gauss消元求G的n-1阶主子式即有

那么把它乘上tmp

get ans

其中注意不要让分母为零

转载传送门

高斯消元时没判g[i][i]<eps 没return abs GG好几发

#include<cmath>
#include<ctime>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<complex>
#include<iostream>
#include<algorithm>
#include<iomanip>
#include<vector>
#include<string>
#include<bitset>
#include<queue>
#include<map>
#include<set>

using namespace std;

typedef double db;

inline int read()
{
	int x=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch<='9'&&ch>='0'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
	return x*f;
}
void print(int x)
{if(x<0)putchar('-'),x=-x;if(x>=10)print(x/10);putchar(x%10+'0');}

const int N=60;const db eps=1e-9;

db p[N][N],g[N][N];

db det(int num)
{
	register int i,j,k;
	db tmp,t;
	for(i=1;i<=num;++i)
	{
		for(j=i;j<=num;++j)if(abs(g[j][i])>eps)break;
		if(j==num+1)return 0;
		if(j!=i)for(k=1;k<=num;++k)swap(g[i][k],g[j][k]);
		tmp=g[i][i];
		for(j=1;j<=num;++j)
		{
			if(j==i)continue;t=g[j][i];
			for(k=i;k<=num;++k)
			{g[j][k]=g[j][k]-g[i][k]*t/tmp;}
		}
		if(abs(g[i][i])<eps)return 0;
	}
	db res=1;
	for(i=1;i<=num;++i)res=res*g[i][i];
	return abs(res);
}

int n;
int main()
{
	n=read();
	register int i,j;
	for(i=1;i<=n;++i)for(j=1;j<=n;++j)scanf("%lf",&p[i][j]);
	db tmp=1,ans;
	for(i=1;i<=n;++i)for(j=i+1;j<=n;++j)tmp=tmp*((1-p[i][j])<eps?eps:(1-p[i][j]));
	for(i=1;i<=n;++i)for(j=1;j<=n;++j)if(i!=j)
	g[i][j]=p[i][j]/((1-p[i][j])<eps?eps:(1-p[i][j]));
	for(i=1;i<=n;++i)for(j=1;j<=n;++j)if(i!=j)g[i][i]-=g[i][j];
	ans=det(n-1);ans*=tmp;
	printf("%.8lf\n",ans);
	return 0;
}
/*
3
0 0.5 0.5
0.5 0 0.5
0.5 0.5 0

0.375
*/