【BZOJ】【P3534】【Sdoi2014】【重建】【题解】【矩阵树定理】

最新推荐文章于 2021-08-19 19:52:55 发布

iamzky

最新推荐文章于 2021-08-19 19:52:55 发布

阅读量2.1k

点赞数

分类专栏： OI 文章标签： bzoj

本文链接：https://blog.csdn.net/iamzky/article/details/41317333

版权

OI 专栏收录该内容

324 篇文章 0 订阅

订阅专栏

传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=3534

dt学了矩阵树定理

邻接矩阵中的的权可以不是1,而是其他权值,比如概率

这样计算出来的就是所有生成树的概率和,即

$\sum_T\prod_{e\in T}P_e$

但是这样不对……

生成一颗生成树T的概率应该是

$\prod_{e\in T}P_e\prod_{e\notin T}(1-P_e)$

接着就是神奇的转换

设G要求的矩阵,P是给出的矩阵

我们令

$G(i,j)=\frac{P(i,j)}{1-P(i,j)}$

$G(i,i)=-\sum_{j\neq i}G(i,j)$

$tmp=\prod_e(1-P_e)$

对G计算n-1阶主子式,即有

$\sum_T\prod_{e\in T}\frac{P_e}{1-P_e}$

那么把它乘上tmp

$\prod_e(1-P_e) \ \sum_T\prod_{e\in T}\frac{P_e}{1-P_e}\\~~ =\sum_T \prod_e(1-P_e)\prod_{e\in T}\frac{P_e}{1-P_e}\\~~ =\sum_T\prod_{e\in T}P_e\prod_{e \notin T}(1-P_e)$

答案就这么出来了!!!!

当P=1时处理需要一点小技巧,把它当做1-eps就可以了

Code:

#include<bits/stdc++.h>
using namespace std;
int n;
const double eps=1e-10;
double A[55][55];
int dcmp(double x){return x<-eps?-1:x>eps;}
double Gauss(){
    double ans=1;
    for(int i=1;i<n;i++){
        int r=i;
        for(int j=i+1;j<n;j++)if(fabs(A[r][i])<fabs(A[j][i]))r=j;
        if(r!=i)for(int j=1;j<n;j++)swap(A[r][j],A[i][j]);
        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]-=A[i][k]*t;
        }
        if(!dcmp(A[i][i]))return 0;
    }for(int i=1;i<n;i++)ans*=A[i][i];
    return fabs(ans);
}
int main(){
    cin>>n;
    double tmp=1;
    for(int i=1;i<=n;i++)
    for(int j=1;j<=n;j++){
        cin>>A[i][j];
        if(i==j)continue;
        if(A[i][j]>1-eps)A[i][j]-=eps;
        if(i<j)tmp*=1-A[i][j];
        A[i][j]=A[i][j]/(1-A[i][j]);
    }for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)if(i!=j)A[i][i]-=A[i][j];
    printf("%.10lf\n",Gauss()*tmp);
    return 0;
}