$\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;
}