洛谷 P3317 [SDOI2014]重建【矩阵树定理】【概率期望】
点击量:381
不是那么简单的矩阵树。
题目描述
T 国有 $N$ 个城市,用若干双向道路连接。一对城市之间至多存在一条道路。
在一次洪水之后,一些道路受损无法通行。虽然已经有人开始调查道路的损毁情况,但直到现在几乎没有消息传回。
幸运的是,此前 T 国政府调查过每条道路的强度,现在他们希望只利用这些信息估计灾情。具体地,给定每条道路在洪水后仍能通行的概率,请计算仍能通行的道路恰有 $N-1$ 条,且能联通所有城市的概率。
输入输出格式
输入格式:
输入的第一行包含整数 $N$。
接下来 $N$ 行,每行 $N$ 个实数,第 $i+1$ 行,$j$ 列的数 $G_{i,j}$ 表示城市 $i$ 与 $j$ 之间仍有道路联通的概率。
输入保证 $G_{i,j}=G_{j,i}$,且 $G_{i,i}=0$;$G_{i,j}$ 至多包含两位小数。
输出格式:
输出一个任意位数的实数表示答案。
你的答案与标准答案相对误差不超过 $10^{-4}$ 即视为正确。
输入输出样例
输入样例#1:
3 0 0.5 0.5 0.5 0 0.5 0.5 0.5 0
输出样例#1:
0.375
说明
$1<N\le 50$
数据保证答案非零时,答案不小于 $10^{-4}$。
题解:
注意到这道题的边是有概率存在。
矩阵树定理可以支持重边,因此它所做的贡献会乘上这条边出现的次数。这一点很容易从矩阵行列式的角度来理解。
不过由于这里是实数,不容易直接联想。我们在其他题目中总是通过取模来减小精度误差,但是我们也可以认为取模后的数变成了另外一个数域。
既然模意义下有意义,那么就可以认为原数代表的东西有意义。
因此我们把实数放在边权上,就可以求出某 $n-1$ 条边同时存在的概率。但是我们也要保证剩下的 $m-n+1$ 条边不存在。因此一棵树出现的概率是 $\prod_{i\in\text{出现的边}}p_i\prod_{i\in\text{没有出现的边}}(1-p_i)$。
前面的 $\prod p_i$ 很容易求出,而后面的并不是非常显然。
考虑到我们知道的信息是有关前面“$i\in\text{出现的边}$”的信息,我们可以通过用全集 $-$ 出现过的边 $=$ 没有出现的边来解决这个问题。
那么可以让前面变成
$$
\prod_{i\in\text{出现的边}}\frac{p_i}{1-p_i}\cdot\prod_{i\in U}(1-p_i)
$$
其中 $U$ 是边的全集。
这样的话,把边权改为 $\frac{p_i}{1-p_i}$ 做矩阵树定理,最后再乘上 $\prod_{i\in U}(1-p_i)$ 就可以了。
然而会出现 $p_i=1$ 的情况,这样的话我们可以把边的两端缩成一个点。不过本题是在实数运算下的,因此可以有一定的误差,这样认为 $p_i=1-\varepsilon$ 就可以了。
发现这样做的话数值差异会比较大,最大可能达到 $\varepsilon^n$。为了进一步减小误差,我们用到高斯消元中的一个技巧,当进行到第 $i$ 行时,可以找第 $i+1$ 到第 $n$ 行中第 $i$ 列元素的绝对值最大的一行与第 $i$ 行交换。
注意矩阵树定理要把基尔霍夫矩阵去掉一行一列。
时间复杂度 $O(n^3)$。
Code:
#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define eps 1e-7
double f[55][55];
int main()
{
int n;
scanf("%d",&n);
double P=1;
for(int i=1;i<=n;++i)
for(int j=1;j<=n;++j)
{
if(i==j)
{
scanf("%*lf");
continue;
}
scanf("%lf",&f[i][j]);
if(fabs(f[i][j]-1)<eps)
f[i][j]-=eps;
if(fabs(f[i][j])<eps)
f[i][j]+=eps;
if(i<j)
P*=(1-f[i][j]);
f[i][j]/=(1-f[i][j]);
f[i][i]+=f[i][j];
f[i][j]=-f[i][j];
}
double ans=1;
for(int i=1;i<n;++i)
{
for(int j=i+1;j<n;++j)
if(fabs(f[j][i])>fabs(f[i][i]))
{
std::swap(f[j],f[i]);
ans=-ans;
}
double t=f[i][i];
ans*=t;
for(int j=i;j<n;++j)
f[i][j]/=t;
for(int j=i+1;j<n;++j)
{
double r=f[j][i];
for(int k=i;k<n;++k)
f[j][k]-=r*f[i][k];
}
}
printf("%.8lf\n",ans*P);
return 0;
}
说点什么