SDOI2014 重建

题目链接:戳我

变元矩阵树定理,具体可以去看蒟蒻学习笔记中[树上问题]这一篇的整理QAQ

原式子为\(\sum_T\prod_{(u,v)\in T}p_{u,v}\prod_{(u,v)\notin T}(1-p_{u,v})\)

\(=\prod_{(u,v)\in G}(1-p_{u,v})\sum_T\prod_{(u,v)\in T}\frac{p_{u,v}}{1-p_{u,v}}\)

所以直接上变元矩阵树定理就可以了.

但是要注意我们中间涉及到除法操作,所以记得把概率为1的设置成为1-eps,概率为0的设置成eps

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define MAXN 110
#define eps 1e-8
using namespace std;
int n,p;
double all=1.0;
double a[MAXN][MAXN];
inline double solve()
{
    double cur_ans=1.0;
    for(int i=2;i<=n;i++)
    {
        int pos=i;
        for(int j=i+1;j<=n;j++)
            if(fabs(a[j][i])>fabs(a[pos][i]))
            pos=j;
        if(pos!=i) swap(a[i],a[pos]),cur_ans*=-1.0;
        for(int j=i+1;j<=n;j++)
        {
            double t=a[j][i]/a[i][i];
            for(int k=i;k<=n;k++) a[j][k]-=t*a[i][k];
        }
    }
    for(int i=2;i<=n;i++)
    cur_ans*=a[i][i];
    return cur_ans;
}
int main()
{
    #ifndef ONLINE_JUDGE
    freopen("ce.in","r",stdin);
    #endif
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=n;j++)
        {
            scanf("%lf",&a[i][j]);
            if(a[i][j]<eps) a[i][j]=eps;
            if(a[i][j]>1.0-eps) a[i][j]=1.0-eps;
            if(i<j) all*=(1.0-a[i][j]);
            a[i][j]=a[i][j]/(1.0-a[i][j]);
        }
    }
    /*
    for(int i=1;i<=n;i++)
    {
    for(int j=1;j<=n;j++)
        printf("%.3lf ",a[i][j]);
    printf("\n");
    }
    */
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            if(i!=j)
                a[i][i]-=a[i][j];
    for(int i=1;i<=n;i++)
    for(int j=1;j<=n;j++)
        a[i][j]=-a[i][j];
    //rintf("all=%.3lf\n",all);
    printf("%.6lf\n",     all*solve( )); 
    return 0;
}