SDOI2014重建-矩阵树定理、概率

标签：int 矩阵 定理 det SDOI2014 double prod rep

link：https://www.luogu.com.cn/problem/P3317
给一张无向图，每条边有一定概率连通，问整张图恰好构成一棵 \(n\) 个点的树的概率。输出实数。
\(1<n\leq 50\)

---

这种问题通常会试着写出来：

\[ans=\sum_{T} (\prod_{e\in T} p_e)(\prod_{e\not\in T} (1-p_e))=\prod_{e\in E} (1-p_e)\sum_T \prod_{e\in T}\frac{p_e}{1-p_e} \]

后面就可以看成是每条边的边权为 \(\frac{p_e}{1-p_e}\) 的图，用矩阵树定理来求（这种对每个生成树求边权积的形式都可以这么做，行列式就是这样算的）

实现的时候需要注意 \(p_e=1\) 的情况，浮点数除0会直接NaN，但是如果给 \(p_e\) 减去一个很小的 \(\epsilon=10^{-8}\) 这样子就可以避免问题。

#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=(a);i<=(b);i++)
using namespace std;
const int N=55;
const double eps=1e-6;
int n;
double p[N][N];
vector<vector<double>> M;
double guass(vector<vector<double>> a,int n){
    double det=1;
    rep(i,1,n){
        int k=i;
        rep(j,i,n)if(fabs(a[j][i])>fabs(a[k][i]))k=j;
        swap(a[i],a[k]);
        if(i!=k)det=-det;
        rep(j,i+1,n){
            double r=a[j][i]/a[i][i];
            rep(k,1,n)a[j][k]=(a[j][k]-a[i][k]*r);
        }
    }
    rep(i,1,n)det=det*a[i][i];
    return det;
}

int main(){
    cin>>n;
    double prod=1;
    rep(i,1,n)rep(j,1,n){
        cin>>p[i][j];
        if(p[i][j]==1)p[i][j]-=eps;
        if(i<j)prod*=1-p[i][j];
        p[i][j]/=(1-p[i][j]);
    }

    M=vector<vector<double>>(n+1,vector<double>(n+1,0));
    rep(i,1,n)rep(j,1,n){
        if(i==j){
            rep(k,1,n)M[i][i]+=p[i][k];
        }else{
            M[i][j]=-p[i][j];
        }
    }
    cout<<fixed<<setprecision(4)<<guass(M,n-1)*prod;
    return 0;
}

标签：int,矩阵,定理,det,SDOI2014,double,prod,rep
From： https://www.cnblogs.com/yoshinow2001/p/18047666