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

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