难度：省选/NOI-

题目描述

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

输入输出格式

输入格式：

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

输出格式：

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

输入输出样例

Sample input

3
0 0.5 0.5
0.5 0 0.5
0.5 0.5 0

Sample ouput

0.375

说明

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

题解

        本质上是求按照给定每条边出现的概率，求得到生成树的概率。本题需要矩阵树前缀知识。
        在构造基比霍夫矩阵时，需要用到邻接矩阵，这里的邻接矩阵在有边时为1，否则为0，同时度数矩阵的对角线数值为相邻边的数目。如果将有边时的值改成权值，同时将度数矩阵的对角线数值改成相邻边的权值之和，再求N-1阶主子式绝对值会得到什么？这里求出的就是：

        $e(i,j)$是边权值，这个结论称为变元矩阵树定理。在有重边的情况下，我们还可以将权值定义为重边数目，即在有重边的情况下变元矩阵树定理同样可以应用。
        回到原问题，现在需要求：

        就是：

        即：

        于是有：

        令权值为：

        这样就可以用矩阵树求了，但还需要注意一个问题：$P(i,j)$为1时分母为0，这时可以将所有值为1的概率减去EPS（通常为$10^{-5}$）再计算。至于行列式的计算可以用高斯消元法化为上三角行列式再处理，时间复杂度$O(n^3)$。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#include <bits/stdc++.h>

using namespace std;
const double EPS = 1e-5;
double p[55][55], K[55][55], op[55][55];
int n;

inline double guess() {
    double ss = 1.0;
    for (int i = 1; i <= n; i++) {
        int maxn = i;
        for (int j = i + 1; j <= n; j++)if (fabs(K[j][i]) > fabs(K[maxn][i]))maxn = j;
        if (i != maxn)swap(K[i], K[maxn]), ss *= -1;
        if (fabs(K[maxn][i]) < EPS)return 0;
        double div = K[i][i];
        for (int j = i; j <= n; j++)K[i][j] /= div;
        ss *= div;
        for (int j = i + 1; j <= n; j++) {
            div = K[j][i];
            for (int z = 1; z <= n; z++)K[j][z] -= K[i][z] * div;
        }
    }
    for (int i = 1; i <= n; i++)ss *= K[i][i];
    return ss;
}

int main() {
    cin >> n;
    for (int i = 1; i <= n; i++)
        for (int j = 1; j <= n; j++) {
            cin >> p[i][j];
            op[i][j] = p[i][j];
            if (fabs(op[i][j] - 1.00) < EPS)op[i][j] = 1.00 - EPS;
        }
    for (int i = 1; i <= n; i++)
        for (int j = 1; j <= n; j++) {
            if (fabs(p[i][j] - 1.00) < EPS)p[i][j] = 1.00 - EPS;
            p[i][j] = p[i][j] / (1 - p[i][j]);
            K[i][i] += p[i][j];
            if (i != j)K[i][j] = -p[i][j];
        }
    --n;
    double ans = fabs(guess());
    ++n;
    for (int i = 1; i <= n; i++) {
        for (int j = i + 1; j <= n; j++)ans *= 1 - op[i][j];
    }
    printf("%.6f", ans);
    return 0;
}