洛谷 P3317 [SDOI2014]重建（矩阵树定理+数学推导） [bzoj3534]

传送门
首先，大家应该都能看出来这是矩阵树定理，然后大部分人应该就会把概率直接带进去算，然后就愉快地WA掉了（我当时就是这么想的，幸亏没交）
然后就来讲这个题的正解思路。
首先我们来看答案应该是怎样的：

ans=∑Tree∏(u,v)∈EP(u,v)∏(u,v)∉E(1−P(u,v))

然后我们来想一下怎么来构造这个答案：首先，我们直接矩阵树用高斯算出来的结果应该是这个：

now=∑Tree∏(u,v)∈EP(u,v)

那我们怎么让它变成答案那个样子呢？直观地，我们可以这样：

now=now∗∏(u,v)(1−P(u,v))

然后就变成了这样子：

now=∑Tree∏(u,v)∈EP(u,v)∗(1−P(u,v))∏(u,v)∉E(1−P(u,v))

现在我们发现now与答案有着中间那一部分的差距，那么如何消除那个差距呢？
我们可以注意到式子中的第一个P(u,v)其实就是矩阵的初值，那么也就是这个：

now=∑Tree∏(u,v)∈EAu,v∗(1−P(u,v))∏(u,v)∉E(1−P(u,v))

那么我们现在就是要让：

Au,v∗(1−P(u,v))→P(u,v)

想必我说到现在大家应该就懂了吧，我们只要这样：

Au,v=P(u,v)(1−P(u,v))

这样我们现在的now值就是答案了！

tmp=∏(u,v)(1−P(u,v))

now=tmp∗∑Tree∏(u,v)∈EP(u,v)(1−P(u,v))

还有就是几个技巧：
当矩阵中出现|a|<eps时就a=eps
当矩阵中出现|1−a|<eps时就a=1−eps
自己想一想为什么。
代码如下：

/**************************************************************
    Problem: 3534
    User: stone41123
    Language: C++
    Result: Accepted
    Time:16 ms
    Memory:1308 kb
****************************************************************/
#include<bits/stdc++.h>
#define ll long long
using namespace std;
int n;
double a[51][51];
double ans;
double eps=1e-8;
void gauss(){
    for(int i=1;i<n;i++){
        int mx=i;
        for(int j=i+1;j<n;j++){
            if(fabs(a[j][i])>fabs(a[mx][i]))mx=j;
        }
        if(mx!=i)for(int j=1;j<n;j++)swap(a[i][j],a[mx][j]);
        for(int k=i+1;k<n;k++){
            double mul=a[k][i]/a[i][i];
            for(int j=i;j<n;j++){
                a[k][j]-=a[i][j]*mul;
            }
        }
        if(fabs(a[i][i])<eps){
            ans=0;
            return;
        }
    }
    for(int i=1;i<n;i++){
        ans*=a[i][i];
    }
    ans=fabs(ans);
}
int main(){
    scanf("%d",&n);
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++){
            scanf("%lf",&a[i][j]);
        }
    }
    ans=1;
    double tmp=1;
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++){
            if(fabs(a[i][j])<eps)a[i][j]=eps;
            if(fabs(1.0-a[i][j])<eps)a[i][j]=1-eps;
            if(i<j)tmp*=1.0-a[i][j];
            a[i][j]=a[i][j]/(1.0-a[i][j]);
        }
    }
    for(int i=1;i<=n;i++){
        a[i][i]=0;
        for(int j=1;j<=n;j++){
            if(i!=j)
            a[i][i]-=a[i][j];
        }
    }
    gauss();
    ans*=tmp;
    printf("%.10lf",ans);
    return 0;
}