测试地址:重建
做法:本题需要用到矩阵树定理。
这两天去学(背)了矩阵树定理,主要就是将度数矩阵
D
D
减去邻接矩阵得到基尔霍夫矩阵,然后将矩阵最后一行和最后一列去掉,剩下的部分求个行列式值,那么这个行列式的值就是这个图的生成树个数。根据行列式的性质,将一行乘上同一个数加到另一行上,所得到的行列式值不变,所以我们可以用类似高斯消元的方法把行列式消成上三角行列式。而上三角行列式的值等于其主对角线上元素的乘积,那么直接运算即可。
矩阵树定理的一个扩展就是,如果把度数换成和该点相连的边权和,把邻接矩阵中的
Aij
A
i
j
改成点
i
i
到点的边的边权和,那么按照上面计算出来的就是图的所有生成树边权积的和。如果不理解上面那句话,上面计算出来的实际上是下面一个东西:
∑T is tree∏(u,v)∈TAuv
∑
T
i
s
t
r
e
e
∏
(
u
,
v
)
∈
T
A
u
v
回到这题本身,我们发现要求的是:
ans=∑T is tree∏(u,v)∈TAuv∏(u,v)∉T(1−Auv)
a
n
s
=
∑
T
i
s
t
r
e
e
∏
(
u
,
v
)
∈
T
A
u
v
∏
(
u
,
v
)
∉
T
(
1
−
A
u
v
)
这就不能直接用上面的方法求了,但是我们可以将边权转化,使得能用上面的方法求出一部分。我们把
∏(u,v)(1−Auv)
∏
(
u
,
v
)
(
1
−
A
u
v
)
提出来,得:
ans=∏(u,v)(1−Auv)⋅∑T is tree∏(u,v)∈TAuv1−Auv
a
n
s
=
∏
(
u
,
v
)
(
1
−
A
u
v
)
⋅
∑
T
i
s
t
r
e
e
∏
(
u
,
v
)
∈
T
A
u
v
1
−
A
u
v
那么我们把边权转化成
Auv1−Auv
A
u
v
1
−
A
u
v
,就可以用矩阵树定理算出后面的和式了,时间复杂度为
O(n3)
O
(
n
3
)
。
为了防止出现奇怪的情况,当
1−Auv<eps
1
−
A
u
v
<
e
p
s
时,我们令
Auv=eps
A
u
v
=
e
p
s
,这样就不会出现分母为
0
0
<script type="math/tex" id="MathJax-Element-2006">0</script>等鬼畜情况了。
以下是本人代码:
#include <bits/stdc++.h>
using namespace std;
const long double eps=1e-6;
int n;
long double p[55][55],tot=1.0;
void gauss()
{
n--;
for(int i=1;i<=n;i++)
{
int now=i;
for(int j=i+1;j<=n;j++)
if (fabs(p[j][i])>fabs(p[now][i])) now=j;
for(int j=1;j<=n;j++)
swap(p[i][j],p[now][j]);
if (fabs(p[i][i])<eps) {p[i][i]=0.0;return;}
for(int j=i+1;j<=n;j++)
{
for(int k=i+1;k<=n;k++)
p[j][k]+=-p[j][i]*p[i][k]/p[i][i];
p[j][i]=0.0;
}
}
}
int main()
{
scanf("%d",&n);
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
{
scanf("%Lf",&p[i][j]);
if (i==j) continue;
if (p[i][j]+eps>1.0) p[i][j]-=eps;
if (i<j) tot*=1-p[i][j];
p[i][j]/=p[i][j]-1;
}
for(int i=1;i<=n;i++)
{
p[i][i]=0.0;
for(int j=1;j<=n;j++)
if (i!=j) p[i][i]-=p[i][j];
}
gauss();
long double ans=1.0;
for(int i=1;i<=n;i++)
ans*=p[i][i];
ans*=tot;
printf("%.6Lf",ans);
return 0;
}