【洛谷-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;
}

cpp

Description\rm{Description}Description#

Solution\rm{Solution}Solution#

矩阵树定理#

分析#

Code\rm{Code}Code#

$\rm{Description}$ #

$\rm{Solution}$ #

$\rm{Code}$ #