题解 P3317 [SDOI2014]重建

题解

前置芝士：深度理解的矩阵树定理

矩阵树定理能求生成树个数的原因是，它本质上求的是：

[sum_T prod_{ein T} w_e ]

其中 (w_e) 是边权，那么我们会发现其实当边权是 (1) 时，本式所求即为生成树个数。

那么回到这题来，这题让求的是

[sum_{T}prod_{ein T}w_eprod_{e otin T}(1-w_e) ]

很容易看出来，这个式子和上面矩阵树的式子很像。让我们来推一波。

[(sum_T prod_{ein T} w_e) × prod_{e}w_e=sum_Tprod_{ein T}w_e(1-w_e)prod_{e otin T}(1-w_e) ]

此式和上面答案更近了一步，我们只需把 (prod_{ein T}(1-w_e)) 消掉即可。

显然

[w_e=(1-w_e)×frac{w_e}{1-w_e} ]

所以，我们要求的就是 (sum_{T}prod_{ein T}frac{w_e}{1-w_e})　

于是我们在初始化 (Laplace) 矩阵时直接以 (frac{w_e}{1-w_e}) 为边权。

有一个需要注意的地方，因为 (w_ein [0,1])，而 (w_e) 为分子，当 (w_e=0) 是，原式趋于无穷小，所以我们可以将其赋为 (eps) ，在分母上的 (1-w_e) 则反过来，若 (w_e=1) 则原式趋于无穷大，所以可以赋其为 (1-eps) 。

至于为什么这样，是为了防止式子在计算机中浮点数例外。

Code

(ACkern 0.5emCODE:)

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#define ri register signed
#define p(i) ++i
using namespace std;
namespace IO{
    char buf[1<<21],*p1=buf,*p2=buf;
    #define gc() p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++
    inline int read() {
        ri x=0,f=1;char ch=getchar();
        while(ch<'0'||ch>'9') {if (ch=='-') f=-1;ch=getchar();}
        while(ch>='0'&&ch<='9') {x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
        return x*f; 
    }
}
using IO::read;
namespace nanfeng{
    #define cmax(x,y) ((x)>(y)?(x):(y))
    #define cmin(x,y) ((x)>(y)?(y):(x))
    #define FI FILE *IN
    #define FO FILE *OUT
    typedef double db;
    static const int N=55;
    static const db eps=1e-8;
    db G[N][N],ans=1.0,tmp=1.0;
    int n;
    inline void Gauss() {
        int tr=0;
        for (ri i(1);i<=n;p(i)) {
            int k=i;
            for (ri j(i+1);j<=n;p(j)) if (fabs(G[j][i])>fabs(G[k][i])) k=j;
            if (k!=i) swap(G[i],G[k]),tr^=1;
            // for (ri j(1);j<=n;p(j)) printf("%.8lf ",G[i][j]);
            // puts("");
            for (ri j(i+1);j<=n;p(j)) {
                db tmp=G[j][i]/G[i][i];
                for (ri l(i);l<=n;p(l)) G[j][l]-=tmp*G[i][l];
            }
            if (fabs(G[i][i])<eps) {ans=0;return;} 
            ans=ans*G[i][i];
            // printf("ans=%.10lf G[i][i]=%.10lf
",ans,G[i][i]);
        } 
        if (tr) ans=-ans;
    }
    inline int main() {
        // FI=freopen("nanfeng.in","r",stdin);
        // FO=freopen("nanfeng.out","w",stdout);
        n=read();
        for (ri i(1);i<=n;p(i)) {
            for (ri j(1);j<=n;p(j)) scanf("%lf",&G[i][j]); 
        }
        for (ri i(1);i<=n;p(i)) {
            for (ri j(1);j<=n;p(j)) {
                if (G[i][j]<eps) G[i][j]=eps;
                if ((1.0-G[i][j])<eps) G[i][j]=1.0-eps;
                if (i<j) tmp*=(1.0-G[i][j]);
                G[i][j]/=(1.0-G[i][j]);
            }
        }
        // printf("tmp=%.10lf
",tmp);
        for (ri i(1);i<=n;p(i)) {
            G[i][i]=0.0;
            for (ri j(1);j<=n;p(j)) {
                if (i!=j) G[i][i]+=G[i][j],G[i][j]=-G[i][j];
            }
        }
        n-=1;
        Gauss();
        printf("%.10lf
",ans*tmp);
        return 0;
    }
}
int main() {return nanfeng::main();}