矩阵树原来还有一个变元矩阵树定理的么?
简单来说就是要求(G)的所有生成树的权值积的和,可以把基尔霍夫矩阵(K)改一改,(K_{i,i})表示与(i)连的所有边的权值和,(K_{i,j})为((i,j))边权的相反数。那么只要求出它的一个主子式的绝对值就是答案了
考虑本题,就是要我们求$$ans=sum_Tprod_{ein T}p_eprod_{e otin T}(1-p_e)$$
考虑一下那个变元矩阵树定理,可以求出$$sum_Tprod_{ein T}w_e$$
因为$$prod_{e otin T}(1-p_e)=frac{prod_e(1-p_e)}{prod_{ein T}(1-p_e)}$$
然后代入上式,可得$$ans=prod_e(1-p_e)sum_Tprod_{ein T}frac{p_e}{1-p_e}$$
那么我们把(frac{p_e}{1-p_e})设为这一条边的权值,用变元矩阵树带进去搞一搞,再乘上前面那一串东西就好了
然而如果(p_e=1)该咋办
因为(frac{p_e}{0}=infty),又有(frac{p_e}{eps}=infty),那就让(p_e)近似等于(1-eps)即可(反正误差很小了)
//minamoto
#include<bits/stdc++.h>
#define rint register int
using namespace std;
const int N=55;const double eps=1e-8;
int n;double pro=1,G[N][N];
inline int cmp(double x){return fabs(x)<eps?0:(x<0?-1:1);}
double det(int n){
int mx;double t,res=1;
for(rint i=1;i<=n;++i){
mx=i;for(rint j=i+1;j<=n;++j)if(cmp(G[mx][i]-G[j][i])<0)mx=i;
if(mx!=i)for(rint j=i;j<=n;++j)swap(G[i][j],G[mx][j]);
if(!G[i][i])return 0;
for(rint j=i+1;j<=n;++j){
t=G[j][i]/G[i][i];
for(rint k=i;k<=n;++k)G[j][k]-=G[i][k]*t;
}res*=G[i][i];
}
return fabs(res);
}
int main(){
// freopen("testdata.in","r",stdin);
scanf("%d",&n);
for(rint i=1;i<=n;++i)for(rint j=1;j<=n;++j){
scanf("%lf",&G[i][j]);if(i==j)continue;
if(G[i][j]>1-eps)G[i][j]-=eps;
if(i<j)pro*=(1-G[i][j]);G[i][j]/=(1-G[i][j]);
}
for(rint i=1;i<=n;++i)for(rint j=1;j<=n;++j)if(i!=j)
G[i][i]+=G[i][j],G[i][j]=-G[i][j];
printf("%.5lf
",det(n-1)*pro);return 0;
}