【BZOJ3534】重建（SDOI2014）-矩阵树定理

测试地址：重建
做法：本题需要用到矩阵树定理。
这两天去学（背）了矩阵树定理，主要就是将度数矩阵D$D$减去邻接矩阵A$A$得到基尔霍夫矩阵，然后将矩阵最后一行和最后一列去掉，剩下的部分求个行列式值，那么这个行列式的值就是这个图的生成树个数。根据行列式的性质，将一行乘上同一个数加到另一行上，所得到的行列式值不变，所以我们可以用类似高斯消元的方法把行列式消成上三角行列式。而上三角行列式的值等于其主对角线上元素的乘积，那么直接运算即可。
矩阵树定理的一个扩展就是，如果把度数换成和该点相连的边权和，把邻接矩阵中的Aij$A_{ij}$改成点i$i$到点j$j$的边的边权和，那么按照上面计算出来的就是图的所有生成树边权积的和。如果不理解上面那句话，上面计算出来的实际上是下面一个东西：
∑T is tree∏(u,v)∈TAuv$\sum _{Tistree}\prod _{(u,v)\in T}{A}_{uv}$
回到这题本身，我们发现要求的是：
ans=∑T is tree∏(u,v)∈TAuv∏(u,v)∉T(1−Auv)$ans=\sum _{Tistree}\prod _{(u,v)\in T}{A}_{uv}\prod _{(u,v)\notin T}(1-A_{uv})$
这就不能直接用上面的方法求了，但是我们可以将边权转化，使得能用上面的方法求出一部分。我们把∏(u,v)(1−Auv)$\prod _{(u,v)}(1-A_{uv})$提出来，得：
ans=∏(u,v)(1−Auv)⋅∑T is tree∏(u,v)∈TAuv1−Auv$ans=\prod _{(u,v)}(1-A_{uv})\cdot \sum _{Tistree}\prod _{(u,v)\in T}\frac{A_{uv}}{1-A_{uv}}$
那么我们把边权转化成Auv1−Auv$\frac{A_{uv}}{1-A_{uv}}$，就可以用矩阵树定理算出后面的和式了，时间复杂度为O(n3)$O(n^3)$。
为了防止出现奇怪的情况，当1−Auv<eps$1-A_{uv}<eps$时，我们令Auv=eps$A_{uv}=eps$，这样就不会出现分母为0$0$等鬼畜情况了。
以下是本人代码：

#include <bits/stdc++.h>
using namespace std;
const long double eps=1e-6;
int n;
long double p[55][55],tot=1.0;

void gauss()
{
    n--;
    for(int i=1;i<=n;i++)
    {
        int now=i;
        for(int j=i+1;j<=n;j++)
            if (fabs(p[j][i])>fabs(p[now][i])) now=j;
        for(int j=1;j<=n;j++)
            swap(p[i][j],p[now][j]);
        if (fabs(p[i][i])<eps) {p[i][i]=0.0;return;}
        for(int j=i+1;j<=n;j++)
        {
            for(int k=i+1;k<=n;k++)
                p[j][k]+=-p[j][i]*p[i][k]/p[i][i];
            p[j][i]=0.0;
        }
    }
}

int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
        {
            scanf("%Lf",&p[i][j]);
            if (i==j) continue;
            if (p[i][j]+eps>1.0) p[i][j]-=eps;
            if (i<j) tot*=1-p[i][j];
            p[i][j]/=p[i][j]-1;
        }
    for(int i=1;i<=n;i++)
    {
        p[i][i]=0.0;
        for(int j=1;j<=n;j++)
            if (i!=j) p[i][i]-=p[i][j];
    }

    gauss();
    long double ans=1.0;
    for(int i=1;i<=n;i++)
        ans*=p[i][i];
    ans*=tot;
    printf("%.6Lf",ans);

    return 0;
}