题解 P3317 【[SDOI2014]重建】
玫葵之蝶
2017-09-21 13:34:26
没人发题解？那我就来水一篇（逃）
首先，大家应该都能看出来这是矩阵树定理，然后大部分人应该就会把概率直接带进去算，然后就愉快地WA掉了（我当时就是这么想的，幸亏没交）
然后就来讲这个题的正解思路。
首先我们来看答案应该是怎样的：
$$ans=\sum_{Tree}\prod_{(u,v)\in E} P_{(u,v)}\prod_{(u,v)\notin E}\big(1-P_{(u,v)}\big)$$然后我们来想一下怎么来构造这个答案：首先，我们直接矩阵树用高斯算出来的结果应该是这个：
$$now=\sum_{Tree}\prod_{(u,v)\in E} P_{(u,v)}$$那我们怎么让它变成答案那个样子呢？直观地，我们可以这样：
$$now=now*\prod_{(u,v)}\big(1-P_{(u,v)}\big)$$然后就变成了这样子：
$$now=\sum_{Tree}\prod_{(u,v)\in E} P_{(u,v)}*\big(1-P_{(u,v)}\big)\prod_{(u,v)\notin E}\big(1-P_{(u,v)}\big)$$现在我们发现now与答案有着中间那一部分的差距，那么如何消除那个差距呢？
我们可以注意到式子中的第一个$P_{(u,v)}$其实就是矩阵的初值，那么也就是这个：
$$now=\sum_{Tree}\prod_{(u,v)\in E} A_{u,v}*\big(1-P_{(u,v)}\big)\prod_{(u,v)\notin E}\big(1-P_{(u,v)}\big)$$那么我们现在就是要让：
$$A_{u,v}*\big(1-P_{(u,v)}\big)\to P_{(u,v)}$$想必我说到现在大家应该就懂了吧，我们只要这样：
$$A_{u,v}=\frac {P_{(u,v)}}{\big(1-P_{(u,v)}\big)}$$这样我们现在的now值就是答案了！
$$tmp=\prod_{(u,v)}\big(1-P_{(u,v)}\big)$$ $$now=tmp*\sum_{Tree}\prod_{(u,v)\in E} \frac {P_{(u,v)}}{\big(1-P_{(u,v)}\big)}$$
还有就是几个技巧：
当矩阵中出现$|a|<eps$时就$a=eps$
当矩阵中出现$|1-a|<eps$时就$a=1-eps$
自己想一想为什么。
看我打公式这么辛苦，就来[我的blog](http://blog.csdn.net/stone41123/article/details/78050473)看看吧^\_^
代码如下：
```cpp
/**************************************************************
Problem: 3534
User: stone41123
Language: C++
Result: Accepted
Time:16 ms
Memory:1308 kb
****************************************************************/
#include<bits/stdc++.h>
#define ll long long
using namespace std;
int n;
double a[51][51];
double ans;
double eps=1e-8;
void gauss(){
for(int i=1;i<n;i++){
int mx=i;
for(int j=i+1;j<n;j++){
if(fabs(a[j][i])>fabs(a[mx][i]))mx=j;
}
if(mx!=i)for(int j=1;j<n;j++)swap(a[i][j],a[mx][j]);
for(int k=i+1;k<n;k++){
double mul=a[k][i]/a[i][i];
for(int j=i;j<n;j++){
a[k][j]-=a[i][j]*mul;
}
}
if(fabs(a[i][i])<eps){
ans=0;
return;
}
}
for(int i=1;i<n;i++){
ans*=a[i][i];
}
ans=fabs(ans);
}
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++){
scanf("%lf",&a[i][j]);
}
}
ans=1;
double tmp=1;
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++){
if(fabs(a[i][j])<eps)a[i][j]=eps;
if(fabs(1.0-a[i][j])<eps)a[i][j]=1-eps;
if(i<j)tmp*=1.0-a[i][j];
a[i][j]=a[i][j]/(1.0-a[i][j]);
}
}
for(int i=1;i<=n;i++){
a[i][i]=0;
for(int j=1;j<=n;j++){
if(i!=j)
a[i][i]-=a[i][j];
}
}
gauss();
ans*=tmp;
printf("%.10lf",ans);
return 0;
}
```