详情看这位神犇的blog
剩下的注释在code里吧.......
#include<iostream> #include<cstdio> #include<cstring> #include<cmath> using namespace std; typedef double db; const db eps=1e-7; db d[55][55],ans=1.0; int n; db det(){//矩阵树定理板子,计算行列式 db re=1.0,div; for(int i=1;i<n;++i){ int x=i; for(int j=i+1;j<n;++j) if(d[j][i]>d[x][i]) x=j; if(x!=i) re=-re,swap(d[x],d[i]); for(int j=i+1;j<n;++j){ div=d[j][i]/d[i][i]; for(int u=i;u<n;++u) d[j][u]-=d[i][u]*div; }re*=d[i][i]; if(fabs(re)<eps) break; }return re; } int main(){ scanf("%d",&n); for(int i=1;i<=n;++i) for(int j=1;j<=n;++j){ scanf("%lf",&d[i][j]); if(i==j) continue; if(fabs(d[i][j])<eps) d[i][j]=eps; if(fabs(1.0-d[i][j])<eps) d[i][j]=1.0-eps;//防止分母出现inf的情况 if(i<j) ans*=(1.0-d[i][j]);//不要重复算 d[i][j]=d[i][j]/(1.0-d[i][j]); } for(int i=1;i<=n;++i){//构造基尔霍夫矩阵 d[i][i]=0.0; for(int j=1;j<=n;++j) if(i!=j) d[i][i]+=d[i][j],d[i][j]=-d[i][j];//度数矩阵+=边权,邻接矩阵为负边权 }ans*=det(); printf("%.10lf",ans); return 0; }