[Luogu P2973] [BZOJ 1778] [USACO10HOL]赶小猪Driving Out the Piggies

#### 洛谷传送门

BZOJ传送门

题目描述

奶牛们建立了一个随机化的臭气炸弹来驱逐猪猡。猪猡的文明包含$1$到$N$ ($2 \leq N \leq 300$)一共$N$个猪城。这些城市由$M$ ($1 \leq M \leq 44,850$)条由两个不同端点$A_j$和$B_j$ ($1 \leq A_j\leq N$; $1 \leq B_j \leq N$)表示的双向道路连接。保证城市$1$至少连接一个其它的城市。一开始臭气弹会被放在城市$1$。每个小时（包括第一个小时），它有$\frac{P}{Q}$ ($1 \leq P \leq 1,000,000$; $1 \leq Q \leq 1,000,000$)的概率污染它所在的城市。如果这个小时内它没有污染它所在的城市，那麽它随机地选择一条道路，在这个小时内沿着这条道路走到一个新的城市。可以离开这个城市的所有道路被选择的概率均等。因为这个臭气弹的随机的性质，奶牛们很困惑哪个城市最有可能被污染。给定一个猪猡文明的地图和臭气弹在每个小时内爆炸的概率。计算每个城市最终被污染的概率。如下例，假设这个猪猡文明有两个连接在一起的城市。臭气炸弹从城市1出发，每到一个城市，它都有$\frac{1}{2}$的概率爆炸。

可知下面这些路径是炸弹可能经过的路径（最后一个城市是臭气弹爆炸的城市）： $1$: 1 $2$: 1-2 $3$: 1-2-1 $4$: 1-2-1-2 $5$: 1-2-1-2-1 … 要得到炸弹在城市$1$终止的概率，我们可以把上面的第$1$，第$3$，第$5$……条路径的概率加起来，（也就是上表奇数编号的路径）。上表中第k条路径的概率正好是$(\frac{1}{2})^k$，也就是必须在前$k-1$个回合离开所在城市（每次的概率为$1 - \frac{1}{2} = \frac{1}{2}$）并且留在最后一个城市（概率为$\frac{1}{2}$）。所以在城市$1$结束的概率可以表示为$\frac{1}{2} + (\frac{1}{2})^3 + (\frac{1}{2})^5 + ...$。当我们无限地计算把这些项一个个加起来，我们最后会恰好得到$\frac{2}{3}$，也就是我们要求的概率，大约是$0.666666667$。这意味着最终停留在城市$2$的概率为$\frac{1}{3}$，大约为$0.333333333$。

输入输出格式

输入格式

第一行四个正整数$N,M,P,Q$。
以下$M$行， 每行两个正整数$A_i,B_i$， 表示编号为$A_i$的城市和$B_i$的城市之间有一条边。

输出格式

一共$N$行， 每行一个保留$9$位的小数，第$i$行的数表示炸弹在第$i$个城市爆炸的概率。

输入输出样例

输入样例#1

2 1 1 2 
1 2 

输出样例#1

0.666666667 
0.333333333 

解题分析

推一波式子： 设$deg_A$为点$A$的度数， 那么大概递推式是这样的：

$$P_A=\sum_{Edge(A,B)}\frac{1-\frac{P}{Q}}{deg_B}P_B$$
然后发现这样搞会有环， 那么就高斯消元搞搞就好了。 因为一开始在$1$号点，所以一号点的值预先设为$1$。

代码如下：

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <cctype>
#include <algorithm>
#include <vector>
#define R register
#define IN inline
#define W while
#define MX 305
#define gc getchar()
#define db double
template <class T>
IN void in(T &x)
{
    x = 0; R char c = gc;
    for (; !isdigit(c); c = gc);
    for (;  isdigit(c); c = gc)
    x = (x << 1) + (x << 3) + c - 48;
}
std::vector <int> to[MX];
int deg[MX];
db mat[MX][MX], pos;
int dot, line, p, q, bd;

void Gauss()
{
    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]);
        for (j = i + 1; j <= bd; ++j) mat[i][j] /= mat[i][i];
        for (j = 1; j <= dot; ++j)
        {
            if(j ^ i)
            {
                for (k = i + 1; k <= bd; ++k)
                mat[j][k] -= mat[j][i] * mat[i][k];
            }
        }
    }
}
int main(void)
{
    int a, b;
    in(dot), in(line), in(p), in(q); bd = dot + 1;
    pos = 1.0 * p / (1.0 * q);
    for (R int i = 1; i <= line; ++i)
    in(a), in(b), to[a].push_back(b), to[b].push_back(a), ++deg[a], ++deg[b];
    mat[1][bd] = 1;
    for (R int i = 1; i <= dot; ++i)
    {
        mat[i][i] = 1;
        for (R int j = 0; j < deg[i]; ++j)
        mat[i][to[i][j]] = (pos - 1.0) / deg[to[i][j]];
    }
    Gauss();
    for (R int i = 1; i <= dot; ++i) printf("%.9lf\n", mat[i][bd] * pos);

}