Luogu – P3317 – [SDOI2014]重建

$\rm{Description}$

给定无向图,每条边以 $p_i$ 的概率存在。求存在的边构成一棵生成树的概率。

点数 $2\leq N\leq 50$,无重边。

$\rm{Solution}$

前置知识:矩阵树定理 – OI Wiki行列式计算 – OI Wiki,题解同步发表于 TonyYin’s Blog

看到 $N\leq 50$ 的范围和生成树,容易想到矩阵树定理

矩阵树定理

对于一个无向图 $G$,设边 $e_i=(u_i, v_i)$,其边权为 $v_i = w(u_i, v_i)$.

定义邻接矩阵 $A$,$A_{i, j}=w(i, j)$.

定义度数矩阵 $D$,$D_{i, i}=\operatorname{deg}_{i}$,存储每个点的度数。

定义基尔霍夫矩阵 $L$ ,$L=D-A$.

任取 $i\in[1, N]$,将基尔霍夫矩阵删去第 $i$ 行和第 $i$ 列,构成一个 $(N-1)\times (N-1)$ 的子矩阵,不妨记为 $L’$.

根据无向图的矩阵树定理,有:

$$
\det L’ = \sum_{T}{\prod_{e\in T}v_e}
$$

对矩阵 $L’$ 计算行列式,得到的结果即为:每棵生成树中,所有边权乘积的和

分析

由概率的性质,我们最终求的答案可以表示为:

$$
\operatorname{Ans}=\sum_T[{\prod_{e\in T}p_e \prod_{e\notin T}(1-p_e)}]
$$

这里一定不要忘记乘上:$\prod_{e\notin T}(1-p_e)$.

注意到 $\operatorname{Ans}$ 的表示式与矩阵树定理很像,考虑如何将其转化成矩阵树定理的标准形式。

要转化为标准形式,我们需要去掉 $\prod_{e\notin T}$ 这个部分,具体转化步骤:

$$
\begin{align}
\operatorname{Ans}&=
\sum_T[{\prod_{e\in T}p_e \prod_{e\notin T}(1-p_e)}]\\
&=\sum_T[{\prod_{e\in T}p_e \frac{\prod_{e}{(1-p_e)}}{\prod_{e\in T}(1-p_e)}}]\\
\end{align}
$$

注意到分子部分为常数,于是:

$$
\begin{align}
\operatorname{Ans}&=
\sum_T[{\prod_{e\in T}p_e \frac{\prod_{e}{(1-p_e)}}{\prod_{e\in T}(1-p_e)}}]\\
&=\prod_{e}{(1-p_e)}\sum_T{\prod_{e\in T}\frac{p_e}{1-p_e}}\\\end{align}
$$

这样,只需要将边权定为 $\dfrac{p_e}{1-p_e}$,之后利用矩阵树定理求解即可。

$\rm{Code}$

实现时,为了防止边权的分母为 $0$,若 $p_e=1$,将 $p_e\leftarrow p_e+\varepsilon$ 即可避免问题发生。

#include <bits/stdc++.h>
using namespace std;
const int MAXN = 100;
const double eps = 1e-10;
int n;
double g[MAXN][MAXN], A[MAXN][MAXN];
double det() {
	for(int col = 2; col <= n; col++) {
		int max_row = col;
		for(int i = col + 1; i <= n; i++)
			if(A[i][col] > A[max_row][col]) max_row = i;
		if(max_row != col) swap(A[col], A[max_row]);
		for(int row = col + 1; row <= n; row++) {
			double t = A[row][col] / A[col][col];
			for(int i = 2; i <= n; i++) {
				A[row][i] -= A[col][i] * t;
			}
		}
	}
	double ans = 1.0;
	for(int i = 2; i <= n; i++) ans *= A[i][i];
	return ans;
}
int main() {
	scanf("%d", &n);
	double prod = 1.0;
	for(int i = 1; i <= n; i++) {
		for(int j = 1; j <= n; j++) {
			scanf("%lf", &g[i][j]); g[i][j] -= eps;
			A[i][j] -= g[i][j] / (1.0 - g[i][j]);
		}
	}
	for(int i = 1; i <= n; i++) {
		for(int j = i + 1; j <= n; j++) {
			A[i][i] += g[i][j] / (1.0 - g[i][j]);
			A[j][j] += g[i][j] / (1.0 - g[i][j]);
			prod *= (1.0 - g[i][j]);
		}
	}
	cout << prod * det() << endl;
	return 0;
}
暂无评论

发送评论 编辑评论

|´・ω・)ノ
ヾ(≧∇≦*)ゝ
(☆ω☆)
(╯‵□′)╯︵┴─┴
 ̄﹃ ̄
(/ω\)
∠( ᐛ 」∠)_
(๑•̀ㅁ•́ฅ)
→_→
୧(๑•̀⌄•́๑)૭
٩(ˊᗜˋ*)و
(ノ°ο°)ノ
(´இ皿இ`)
⌇●﹏●⌇
(ฅ´ω`ฅ)
(╯°A°)╯︵○○○
φ( ̄∇ ̄o)
ヾ(´・ ・`。)ノ"
( ง ᵒ̌皿ᵒ̌)ง⁼³₌₃
(ó﹏ò。)
Σ(っ °Д °;)っ
( ,,´・ω・)ノ"(´っω・`。)
╮(╯▽╰)╭
o(*////▽////*)q
>﹏<
( ๑´•ω•) "(ㆆᴗㆆ)
😂
😀
😅
😊
🙂
🙃
😌
😍
😘
😜
😝
😏
😒
🙄
😳
😡
😔
😫
😱
😭
💩
👻
🙌
🖕
👍
👫
👬
👭
🌚
🌝
🙈
💊
😶
🙏
🍦
🍉
😣
Source: github.com/k4yt3x/flowerhd
颜文字
Emoji
小恐龙
花!
上一篇
下一篇