[BZOJ3534][Sdoi2014]重建（矩阵树定理）

最新推荐文章于 2022-12-14 11:15:40 发布

xyz32768

最新推荐文章于 2022-12-14 11:15:40 发布

阅读量210

点赞数

CC 4.0 BY-SA版权

分类专栏： BZOJ、UOJ、LOJ 文章标签： Matrix-Tree 卡精度

本文链接：https://blog.youkuaiyun.com/xyz32768/article/details/80150532

BZOJ、UOJ、LOJ 专栏收录该内容

341 篇文章

订阅专栏

本文介绍了一种计算洪水后城市间道路连通概率的方法。通过构建概率图模型，结合矩阵树定理，实现对特定条件下城市网络连通性的评估。

Description

T 国有 $N$ 个城市，用若干双向道路连接。一对城市之间至多存在一条道路。
在一次洪水之后，一些道路受损无法通行。虽然已经有人开始调查道路的损毁情况，但直到现在几乎没有消息传回。
幸运的是，此前 T 国政府调查过每条道路的强度，现在他们希望只利用这些信息估计灾情。具体地，给定每条道路在洪水后仍能通行的概率，请计算仍能通行的道路恰有 $N-1$ 条，且能联通所有城市的概率。

Input

输入的第一行包含整数 $N$ 。
接下来 $N$ 行，每行 $N$ 个实数，第 $i+1$ 行，列的数 $G[i][j]$ 表示城市 $i$ 与 $j$ 之间仍有道路联通的概率。
输入保证 $G[i][j]=G[j][i]$ ，且 $G[i][j]=0$ ； $G[i][j]$ 至多包含两位小数。

Output

输出一个任意位数的实数表示答案。
你的答案与标准答案相对误差不超过 $10^{-4}$ 即视为正确。

Sample Input

3
0 0.5 0.5
0.5 0 0.5
0.5 0.5 0

Sample Output

0.375

HINT

$1<N\le50$
数据保证答案非零时，答案不小于 $10^{-4}$

Solution

题意就是求（ $e$ 表示一条边， $T$ 表示原图的一棵生成树， $P(e)$ 表示边 $e$ 出现的概率）：

a n s = \sum_{T} \prod_{e \in T} P (e) \prod_{e \notin T} (1 - P (e))

$ans=\sum_{T}\prod_{e\in T}P(e)\prod_{e\notin T}(1-P(e))$
发现既要考虑生成树上的边，又要考虑生成树外的边，不好做，因此先假设所有的边都不存在：

n o w = \prod e (1 - P (e))

$now=\prod_{e}(1-P(e))$
然后求生成树：

d e l t a = \sum T \prod e \in T P ( e ) 1 - P ( e )

$delta=\sum_{T}\prod_{e\in T}\frac{P(e)}{1-P(e)}$
将上面两者相乘：

n o w \times d e l t a = \sum T \prod e (1 - P (e)) \prod e \in T P ( e ) 1 - P ( e )

$now\times delta=\sum_T\prod_{e}(1-P(e))\prod_{e\in T}\frac{P(e)}{1-P(e)}$

= \sum T \prod e \notin T (1 - P (e)) \prod e \in T (1 - P (e)) \prod e \in T P ( e ) 1 - P ( e )

$=\sum_T\prod_{e\notin T}(1-P(e))\prod_{e\in T}(1-P(e))\prod_{e\in T}\frac{P(e)}{1-P(e)}$

= \sum T \prod e \in T (1 - P (e)) P ( e ) 1 - P ( e ) \prod e \notin T (1 - P (e))

$=\sum_T\prod_{e\in T}(1-P(e))\frac{P(e)}{1-P(e)}\prod_{e\notin T}(1-P(e))$

= \sum T \prod e \in T P (e) \prod e \notin T (1 - P (e)) = a n s

$=\sum_T\prod_{e\in T}P(e)\prod_{e\notin T}(1-P(e))=ans$
而现在就是要求出

deltadelta $delta$ 。
将 Matrix-Tree 矩阵树定理进行扩展，可以得到：
如果定义矩阵

DD $D$ ，

D [i] [i]

$D[i][i]$ 为从

ii $i$ 发出的边权之和，

D [i] [j] = 0 (i \neq j)

$D[i][j]=0(i\ne j)$
矩阵

AA $A$ 为邻接矩阵，

A [i] [j]

$A[i][j]$ 为边

(i,j)(i,j) $(i,j)$ 的边权，
矩阵

C=D−AC=D−A $C=D-A$ ，
那么

CC $C$ 的任意主子式就等于图的所有生成树的边权积之和。
也就是

\sum_{T} \prod_{(i, j) \in T} A [i] [j]

$\sum_T\prod_{(i,j)\in T}A[i][j]$
所以，就可以建一个新图，边

ee $e$ 的权值为

\frac{P (e)}{1 - P (e)}

$\frac{P(e)}{1-P(e)}$ 。
求出

CC $C$ 矩阵后用高斯消元求解行列式就可以得到

d e l t a

$delta$ 。
此外，

P(e)P(e) $P(e)$ 有可能是

00 $0$ 或

1

$1$ ，这样会出现除以

00 $0$ 的计算。
因此要用一个极小值

e p s

$eps$ （一般取

10−810−8 $10^{-8}$ 或

10−910−9 $10^{-9}$ ），
当

P(e)=0P(e)=0 $P(e)=0$ 时，用

epseps $eps$ 代替

00 $0$ ，
当

P (e) = 1

$P(e)=1$ 时，用

1−eps1−eps $1-eps$ 代替

11 $1$ 。
复杂度：

O (n^{3})

$O(n^3)$

Code

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define For(i, a, b) for (i = a; i <= b; i++)
using namespace std; const int N = 55; const double eps = 1e-8;
int n; double p[N][N], mat[N][N], ans0 = 1.0;
double solve(int n) {
    int i, j, k; For (i, 1, n) {
        int p = i; For (j, i + 1, n)
            if (fabs(mat[j][i]) > fabs(mat[p][i])) p = j;
        if (p != i) For (j, i, n) swap(mat[i][j], mat[p][j]);
        if (fabs(mat[i][i]) < eps) return 0; For (j, i + 1, n) {
            double tmp = mat[j][i] / mat[i][i];
            For (k, i, n) mat[j][k] -= mat[i][k] * tmp;
        }
    }
    double ans = 1; For (i, 1, n) ans *= mat[i][i]; return fabs(ans);
}
int main() {
    int i, j; cin >> n; For (i, 1, n) For (j, 1, n) scanf("%lf", &p[i][j]),
        i != j && p[i][j] < eps ? p[i][j] = eps : 0,
        i != j && 1.0 - p[i][j] < eps ? p[i][j] = 1.0 - eps : 0;
    For (i, 1, n) For (j, i + 1, n) ans0 = ans0 * (1.0 - p[i][j]);
    For (i, 1, n) For (j, 1, n) if (i != j) mat[i][j] = -p[i][j] / (1.0 - p[i][j]);
    For (i, 1, n) For (j, 1, n) if (i != j) mat[i][i] -= mat[i][j];
    printf("%.10lf\n", solve(n - 1) * ans0); return 0;
}