洛谷 P2973 [USACO10HOL]Driving Out the Piggies G ，概率论+高斯消元

洛谷 P2973 [USACO10HOL]Driving Out the Piggies G

题目描述

The Cows have constructed a randomized stink bomb for the purpose of driving away the Piggies. The Piggy civilization consists of N (2 <= N <= 300) Piggy cities conveniently numbered 1..N connected by M (1 <= M <= 44,850) bidirectional roads specified by their distinct endpoints A_j and B_j (1 <= A_j <= N; 1 <= B_j <= N). Piggy city 1 is always connected to at least one other city.

The stink bomb is deployed in Piggy city 1. Each hour (including the first one), it has a P/Q (1 <= P <= 1,000,000; 1 <= Q <=

1,000,000; P <= Q) chance of polluting the city it occupies. If it does not go off, it chooses a random road out of the city and follows it until it reaches a new city. All roads out of a city are equally likely to be chosen.

Because of the random nature of the stink bomb, the Cows are wondering which cities are most likely to be polluted. Given a map of the Piggy civilization and the probability that the stink bomb detonates in a given hour, compute for each city the probability that it will be polluted.

For example, suppose that the Piggie civilization consists of two cities connected together and that the stink bomb, which starts in city 1, has a probability of 1/2 of detonating each time it enters a city:

1--2
We have the following possible paths for the stink bomb (where the last entry is the ending city):

1: 1
2: 1-2
3: 1-2-1

4: 1-2-1-2

5: 1-2-1-2-1

etc.
To find the probability that the stink bomb ends at city 1, we can add up the probabilities of taking the 1st, 3rd, 5th, ... paths above (specifically, every odd-numbered path in the above list). The probability of taking path number k is exactly (1/2)^k - the bomb must not remain in its city for k - 1 turns (each time with a probability of 1 - 1/2 = 1/2) and then land in the last city

(probability 1/2).

So our probability of ending in city 1 is represented by the sum 1/2 + (1/2)^3 + (1/2)^5 + ... . When we sum these terms infinitely, we will end up with exactly 2/3 as our probability, approximately 0.666666667. This means the probability of landing in city 2 is 1/3, approximately 0.333333333.

Partial feedback will be provided for your first 50 submissions.

一个无向图，节点1有一个炸弹，在每个单位时间内，有p/q的概率在这个节点炸掉，有1-p/q的概率随机选择一条出去的路到其他的节点上。问最终炸弹在每个节点上爆炸的概率。

输入格式

* Line 1: Four space separated integers: N, M, P, and Q

* Lines 2..M+1: Line i+1 describes a road with two space separated integers: A_j and B_j

输出格式

* Lines 1..N: On line i, print the probability that city i will be destroyed as a floating point number. An answer with an absolute error of at most 10^-6 will be accepted (note that you should output at least 6 decimal places for this to take effect).

样例 #1

样例输入 #1

2 1 1 2 
1 2

样例输出 #1

0.666666667 
0.333333333

---

---

题意

给出一张\(n\)个点，\(m\)条边的无向图，节点\(1\)有一个炸弹，在每个单位时间内，有\(\dfrac {p}{ q}\)的概率在这个节点炸掉，有\(1-\dfrac{p}q\)的概率随机选择一条出去的路到其他的节点上，去每个节点的概率相等。问最终炸弹在每个节点上爆炸的概率。

题解

由于只要一直不爆炸，炸弹就可以在节点间反复横跳，不断对其爆炸概率产生贡献，故通过模拟顺推的方式求概率不太现实(除非可以归纳出无限次数下的表达式)。

这里考虑概率最终稳态下，通过相邻点之间的概率关系来解出每个点的爆炸概率

令

​ 爆炸概率\(P=\dfrac {p}{ q}\)，

​ \(f(i)\)为炸弹在\(i\)节点爆炸的概率，

​ \(arrve(i)\)为炸弹到达i节点的概率

​ \(deg(i)\)为与\(i\)相连节点个数

则： $$f(i)=arrve(i)*P $$

​ $$\quad\quad\quad\quad =P*\sum_{i,j相连}arrve(j)*(1-P)/deg(j) $$

​ $$\quad\quad\quad\quad =P*\sum_{i,j相连}(f(j)/P)*(1-P)/deg(j) $$

​ $$ =\sum_{i,j相连}f(j)*(1-P)/deg(j) $$

而对于\(f(1)\),除加上与之相连的节点为其贡献的概率外，由于炸弹初始于点\(1\)，故还应加上第一轮爆炸概率\(P\)

共\(n\)个点，可写出\(n\)个方程，然后用高斯消元模板求解即可（需注意控制精度）

//https://www.luogu.com.cn/problem/P2973
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef double db;
const db eps=1e-14,one=1.00000000000;
#define rep(i_,from_,to_) for(register int i_=(from_);i_<=(to_);++i_)
#define nep(i,from,to) for(register int i=from;i>=to;--i)
int n,m,p,q,deg[310];
double P,a[310][310];
vector<int> e[310];
inline void addE(int u,int v){e[u].push_back(v);deg[u]++;}
int gauss(int n) {
	int c, r;
	for (c = 0, r = 0; c < n; c ++ ) {
		int t = r;
		for (int i = r; i < n; i ++ ) // 找到绝对值最大的行
			if (fabs(a[i][c]) > fabs(a[t][c]))
				t = i;
		if (fabs(a[t][c]) < eps) continue;
		for (int i = c; i <= n; i ++ ) swap(a[t][i], a[r][i]); // 将绝对值最大的行换到最顶端
		for (int i = n; i >= c; i -- ) a[r][i] /= a[r][c]; // 将当前上的首位变成1
		for (int i = r + 1; i < n; i ++ ) // 用当前行将下面所有的列消成0
			if (fabs(a[i][c]) > eps)
				for (int j = n; j >= c; j -- )
					a[i][j] -= a[r][j] * a[i][c];
		r ++ ;
	}
	if (r < n) {
		for (int i = r; i < n; i ++ )
			if (fabs(a[i][n]) > eps)
				return 2; // 无解
		return 1; // 有无穷多组解
	}
	for (int i = n - 1; i >= 0; i -- )
		for (int j = i + 1; j < n; j ++ )
			a[i][n] -= a[i][j] * a[j][n];
	return 0; // 有唯一解
}
void solve() {
	cin>>n>>m>>p>>q;
	P=((double)p)/q;
	rep(i,1,m){
		int u,v;
		scanf("%d%d",&u,&v);
		addE(u,v);addE(v,u);
	}
	a[0][n]=P;
	rep(i,1,n){
		a[i-1][i-1]=1;
		for(auto j:e[i]){
			a[i-1][j-1]=-(one-P)/deg[j];
		}
	}
	int ret=gauss(n);
	if(ret==0){
		rep(i,1,n){
			printf("%.10lf\n",a[i-1][n]);
		}
	}
}
int main() {
	solve();
	return 0;
}