_WA自动机的blog

过期高中生。

[BZOJ3534][SDOI2014]重建

发表于 更新于 分类于 阅读次数： Valine：

题目大意：

你有一个无向完全图。在一场洪水（？）以后，每条边都可能会损毁。给定每条边没有被损毁的概率，求最后剩下的边正好组成原来的图的一棵生成树的概率。

分析

设一条边\(w\)损毁的概率为\(E(w)\)。由题意，我们要求的其实是 \[\sum_{T \ is \ a \ spanning \ tree \ of \ G}\prod_{e\in T}E(e)\prod_{e\in \complement_G^T}(1-E(e))\ \ \ (*)\] 把这个柿子写出来，我们发现它很像变元矩阵树定理..所以考虑怎么化成矩阵树定理的形式.
我们是要把它化为\(\sum_T \prod_{e\in T} w(e)\)的形式...那么必须把后面一项向T中的边靠拢。怎么做呢？
我们将后面一个连乘柿子变换一下形式... \[ \prod_{e\in \complement_G^T}(1-E(e))=\prod_{e \in G}(1-E(e))*\prod_{e \in T}\frac{1}{(1-E(e))} \] 则.. \[(*)=\prod_{e \in G}(1-E(e)) *\sum_T \prod_{e\in T}\frac{E(e)}{1-E(e)} \] 这就是矩阵树定理板子了.. 注意一个细节：当概率等于1的时候，我们会算出来nan..于是我想缩点...其实并不需要复杂的分类讨论、缩点之类的，只需要将1减去一个eps就行了。造成的精度误差不会太大的^_^。 然后我们就可以愉快地AC辣！

参考代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
#include <cstdio>
#include <cmath>
#include <algorithm>

using std::swap;
using std::fabs;

const int maxn=1e3+1;
const double eps=1e-6;

double G[maxn][maxn];

inline double determinant(double (*A)[maxn],int n)
{
    for (int i=1,c=1,j;i<=n;++i)
    {
        for (j=c;j<=n && fabs(A[j][i])<eps;++j);
        if (j==n+1) continue;
        for (int k=1;k<=n;++k) swap(A[c][k],A[j][k]);
        for (int j=c+1;j<=n;++j)
            if (fabs(A[j][i])>eps)
            {
                double t=A[j][i]/A[c][i];
                for (int k=i;k<=n;++k)
                    A[j][k]-=A[c][k]*t;
            }
        ++c;
    }
    double ans=1;
    for (int i=1;i<=n;++i) ans*=A[i][i];
    return fabs(ans);
}

int main()
{
    int n;
    scanf("%d",&n);
    for (int i=1;i<=n;++i)
        for (int j=1;j<=n;++j)
        {
            scanf("%lf",&G[i][j]);
            if (G[i][j]==1) G[i][j]-=eps;
        }
    double prod=1;
    for (int i=1;i<=n;++i)
        for (int j=1;j<i;++j)
            prod*=1-G[i][j];
    for (int i=1;i<=n;++i)
        for (int j=1;j<=n;++j)
            if (i!=j) G[i][j]=G[i][j]/(1-G[i][j]),G[i][i]+=G[i][j],G[i][j]=-G[i][j];
    printf("%.10lf",prod*determinant(G,n-1));// BZOJ SPJ好像挂了，必须输出10位小数才能过...
}