题目
分析
问题要求
∑
T
(
∏
e
∈
T
p
e
∏
e
∉
T
(
1
−
p
e
)
)
\sum_{T} \left( \prod_{e \in T} p_e \prod_{e \notin T} (1 - p_e)\right)
T∑⎝⎛e∈T∏pee∈/T∏(1−pe)⎠⎞ 由于矩阵树只能处理
e
∈
T
e \in T
e∈T,所以把式子变成
∑
T
(
∏
e
∈
T
p
e
∏
e
(
1
−
p
e
)
∏
e
∈
T
(
1
−
p
e
)
)
=
∏
e
(
1
−
p
e
)
∑
T
(
∏
e
∈
T
p
e
1
−
p
e
)
\begin{aligned} & \sum_{T} \left( \prod_{e \in T} p_e \frac{\prod_{e} (1 - p_e)}{\prod_{e \in T} (1 - p_e)} \right) \\ =& \prod_{e} (1 - p_e) \sum_{T} \left( \prod_{e \in T} \frac{p_e }{1 - p_e} \right)\end{aligned}
=T∑(e∈T∏pe∏e∈T(1−pe)∏e(1−pe))e∏(1−pe)T∑(e∈T∏1−pepe) 于是用矩阵树定理推广计算
∑
T
(
∏
e
∈
T
p
e
1
−
p
e
)
\sum_{T} \left( \prod_{e \in T} \frac{p_e }{1 - p_e} \right)
∑T(∏e∈T1−pepe),最后乘上
∏
e
(
1
−
p
e
)
\prod_{e} (1 - p_e)
∏e(1−pe) 即为答案。
注意
p
e
=
1
p_e = 1
pe=1 会除零而出问题,给
p
e
p_e
pe 减个
ϵ
\epsilon
ϵ 即可。
矩阵树定理推广:计算所有生成树边权乘积的和,即 ∑ T ( ∏ e ∈ T x e ) \sum_{T} \left( \prod_{e \in T} x_e \right) T∑(e∈T∏xe) 其中 x e x_e xe 为一个有理数。
令 x e = p e t x_e = \dfrac{p_e}{t} xe=tpe,其中 p e , t ∈ Z p_e, t \in \Z pe,t∈Z,那么 ∑ T ( ∏ e ∈ T x e ) = 1 t n − 1 ∑ T ( ∏ e ∈ T p e ) \begin{aligned} & \sum_{T} \left( \prod_{e \in T} x_e \right) \\ =& \frac{1}{t^{n - 1}} \sum_{T} \left( \prod_{e \in T} p_e \right)\end{aligned} =T∑(e∈T∏xe)tn−11T∑(e∈T∏pe) 当边权 p p p为整数的时候很好计算:把 u , v u, v u,v 之间连上 p p p 条边,然后计算生成树的个数就是原图所有生成树边权乘积和。
但是 t n − 1 t^{n - 1} tn−1 很显然精度不够,所以我们用 行列式的可乘性 把 t n − 1 t^{n - 1} tn−1 分配给行列式中每个元素,于是我们有了一个新的基尔霍夫矩阵:度为连出去所有边的边权和,边权为原边权的相反数。然后正常矩阵树定理得到的答案就是所求的了。
代码
#include <bits/stdc++.h>
const int MAXN = 50;
const double eps = 1e-5;
int N;
double Mat[MAXN + 5][MAXN + 5];
void AddEdge(int u, int v, double w) {
Mat[u][u] += w, Mat[v][v] += w;
Mat[u][v] -= w, Mat[v][u] -= w;
}
inline double Abs(const double &x) {
return x < 0 ? -x : x;
}
inline int Comp(const double &x, const double &y) {
if (Abs(x - y) < eps)
return 0;
return x > y ? 1 : -1;
}
double Det(int n) {
n--;
double ret = 1;
for (int i = 1; i <= n; i++) {
int p = i;
for (int j = i; j <= n; j++)
if (Comp(Mat[p][i], Mat[j][i]) < 0)
p = j;
if (p != i)
std::swap(Mat[i], Mat[p]);
if (!Comp(Mat[i][i], 0))
return 0;
ret *= Mat[i][i];
for (int j = i + 1; j <= n; j++) {
double tmp = Mat[j][i] / Mat[i][i];
for (int k = i; k <= n; k++)
Mat[j][k] -= Mat[i][k] * tmp;
}
}
return Abs(ret);
}
int main() {
scanf("%d", &N);
double All = 1.0;
for (int i = 1; i <= N; i++)
for (int j = 1; j <= N; j++) {
double p; scanf("%lf", &p);
if (i < j) {
if (Comp(p, 1) == 0)
p -= eps;
AddEdge(i, j, p / (1 - p));
All *= (1 - p);
}
}
printf("%.5f", Det(N) * All);
}