[Luogu P3317] [BZOJ 3534] [SDOI2014] 重建

洛谷传送门

BZOJ传送门

题目描述

$T$国有$N$个城市，用若干双向道路连接。一对城市之间至多存在一条道路。

在一次洪水之后，一些道路受损无法通行。虽然已经有人开始调查道路的损毁情况，但直到现在几乎没有消息传回。

幸运的是，此前T国政府调查过每条道路的强度，现在他们希望只利用这些信息估计灾情。具体地，给定每条道路在洪水后仍能通行的概率，请计算仍能通行的道路恰有$N-1$条，且能联通所有城市的概率。

输入输出格式

输入格式：

输入的第一行包含整数$N$。

接下来$N$行，每行$N$个实数，第$i+l$行，列的数$G[i][j]$表示城市$i$与$j$之间仍有道路联通的概率。

输入保证$G[i][j]=G[j][i]$，且$G[i][j]=0$；$G[i][j]$至多包含两位小数。

输出格式：

输出一个任意位数的实数表示答案。

你的答案与标准答案相对误差不超过$10^{(-4)}$即视为正确。

输入输出样例

输入样例#1：

3
0 0.5 0.5
0.5 0 0.5
0.5 0.5 0

输出样例#1：

0.375

说明

$1 < N < =50$

数据保证答案非零时，答案不小于$10^{(-4)}$

解题分析

首先我们考虑如果将构成基尔霍夫矩阵的邻接矩阵中$A[i][j]$换成$P[i][j]$， $B[i][i]$换为$\sum_{i=1}^nP[i][j]$的话，我们实际上就求出了所有生成树树上边联通概率乘积的和。

但这样是有问题的， 我们实际上要求的是：

$$\sum_{S}\prod_{E(A,B)\in S}P[A][B]\prod_{E(A,B)\notin S}(1-P[A][B])$$
而直接算的话等于只求了
$$\sum_{S}\prod_{E(A,B)\in S}P[A][B]$$
少了后面一堆， 怎么办?首先我们可以乘上这么一团：
$$\prod_{E(A,B)}(1-P[A][B])$$
然后就变成了这个样子：
$$\sum_{S}\prod_{E(A,B)\in S}(P[A][B]\times (1-P[A][B]))\prod_{E(A,B)\notin S}(1-P[A][B])$$
这下前面多出来一团， 但 $P[A][B]$实际上就是我们在矩阵中填入的数， 所以我们把基尔霍夫矩阵中 $T[A][B]$中元素换成 $\frac {P[A][B]}{1-P[A][B]}$即可。

注意到部分$P[A][B]=1$， 所以当这种情况的时候我们把$P[A][B]$换为$P[A][B]-\epsilon$即可保证不出现$INF$的情况。

代码如下：

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <cmath>
#include <algorithm>
#define R register
#define IN inline
#define gc getchar()
#define W while
#define EPS 1e-7
#define db double
#define MX 205
int dot;
db mat[MX][MX], sum;
IN db Gauss()
{
    db ret = 1, mul;
    R int i, j, k;
    for (i = 1; i < dot; ++i)
    {
        for (j = i + 1; j < dot; ++j)
        if(fabs(mat[j][i]) > fabs(mat[i][i])) std::swap(mat[j], mat[i]);
        if(fabs(mat[i][i] < EPS)) return 0;
        for (j = i + 1; j < dot; ++j)
        {
            mul = mat[j][i] / mat[i][i];
            for (k = i; k < dot; ++k) mat[j][k] -= mul * mat[i][k];
        }
        ret *= mat[i][i];
    }
    return ret;
}
int main(void)
{
    scanf("%d", &dot); sum = 1;
    for (R int i = 1; i <= dot; ++i)
    {
        for (R int j = 1; j <= dot; ++j)
        {
            scanf("%lf", &mat[i][j]);
            if(i == j) continue;
            if(mat[i][j] > 1.0 - EPS) mat[i][j] -= EPS;
            if(i < j) sum *= (1.0 - mat[i][j]);//只算一次
            mat[i][j] = mat[i][j] / (1 - mat[i][j]);
        }
    }
    for (R int i = 1; i <= dot; ++i)
    {
        for (R int j = 1; j <= dot; ++j)
        if(i ^ j) mat[i][i] += mat[i][j], mat[i][j] = -mat[i][j];
    }
    printf("%.6lf", sum * Gauss());
}