[USACO10HOL]赶小猪Driving Out the Pigg：概率+高斯消元

看了一上午，参考了了洛谷Siyuan的两篇博客，终于感觉有点理解了。【HN013】游走，「Luogu 2973」Dotp; 下面主要就是cc其blog了。

题目描述

有一个 n​ 个点 m​ 条边的无向图，节点 1​ 有一个炸弹，在每个单位时间内，这个炸弹有$\frac{P}{Q}$的概率在这个节点炸掉，有 1−$\frac{P}{Q}$ 的概率等概率选择一条与当前节点相连的边到其他的节点上。求炸弹在每个节点上爆炸的概率。

数据范围：2≤n≤300，1≤m≤44850，1≤P,Q≤$10^6$

思路

我们设置第i个点期望经过的次数为$x_i$,那么在第i点爆炸的概率为$x_i\ast \frac{P}{Q}$.我们设$de{g}_i$为i结点的度。 假设有图 1-2,2-3,3-1.初始状态在1点。 $$那么x1次数的由来包括：\{\begin{array}{r}1.初始状态1次；\\ 2.从结点2转移来；\\ 3.从结点3转移来。\end{array}$$ 从结点2到结点1的次数是：在结点2不爆炸的概率（$1-\frac{P}{Q}$)下，有$\frac{1}{de{g}_2}$的概率到达结点1。 同理，从结点3到达结点1的次数是：在结点3不爆炸的概率（1-$\frac{P}{Q}$)下，有$\frac{1}{de{g}_3}$的概率到达结点1。 因此得到式子：$$x_1=1+x_2\ast (1-\frac{P}{Q})\ast \frac{1}{de{g}_2}+x_3\ast (1-\frac{P}{Q})\ast \frac{1}{de{g}_3}$$ 也就是：$$x_1-x_2\ast (1-\frac{P}{Q})\ast \frac{1}{de{g}_2}-x_3\ast (1-\frac{P}{Q})\ast \frac{1}{de{g}_3}=1$$ $$那么x2次数的由来包括：\{\begin{array}{l}1.从结点1转移而来；\\ 2.从结点3转移而来。\end{array}$$ 从结点1到结点2的次数是：在结点1不爆炸的概率（$1-\frac{P}{Q}$)下，有$\frac{1}{de{g}_1}$的概率到达结点2。 同理，从结点3到达结点2的次数是：在结点3不爆炸的概率（1-$\frac{P}{Q}$)下，有$\frac{1}{de{g}_3}$的概率到达结点1。 因此得到式子：$$x_2=x_1\ast (1-\frac{P}{Q})\ast \frac{1}{de{g}_1}+x_3\ast (1-\frac{P}{Q})\ast \frac{1}{de{g}_3}$$ 也就是：$$-x_1\ast (1-\frac{P}{Q})\ast \frac{1}{de{g}_1}+x_2-x_3\ast (1-\frac{P}{Q})\ast \frac{1}{de{g}_3}=0$$

后面的每个结点都可以列出一个式子…… 那么就会得到n个未知数，那个等式的方程组，然后利用高斯消元解方程。 形象的理解之后，可以写出这个方程的通项表达式：

那么 $$x_i=\{\begin{array}{llc}{x}_1={\sum }_{(i,j)\in E}(1-\frac{P}{Q})\ast \frac{x_j}{de{g}_j}+1& & i=1\\ x_i={\sum }_{(i,j)\in E}(1-\frac{P}{Q})\frac{x_j}{de{g}_j}& & 1<i<=n\end{array}$$

所以算法的步骤就是： 1.构造矩阵，,矩阵对角线的元素为1，$x_j$到$x_i$的系数$a_{ij}$为$-(1-\frac{P}{Q})\ast \frac{1}{de{g}_j}$ 方程的解为n*1的行列式（1,0，0……） 2.高斯消元求解方程组。 3.输出答案：$x_i\ast \frac{P}{Q}$ 参考代码：

#include<bits/stdc++.h> using namespace std; const int N=305; int n,m,p,q,e[N][N],deg[N]; double a[N][N],b[N],x[N]; void Gauss(int n) { for(int i=1;i<=n;++i) { int p=i; for(int k=i+1;k<=n;++k) if(fabs(a[k][i])>fabs(a[p][i])) p=k; if(i!=p) swap(a[i],a[p]),swap(b[i],b[p]); for(int k=i+1;k<=n;++k) { double d=a[k][i]/a[i][i]; b[k]-=d*b[i]; for(int j=i;j<=n;++j) a[k][j]-=d*a[i][j]; } } for(int i=n;i>=1;--i) { for(int j=i+1;j<=n;++j) b[i]-=x[j]*a[i][j]; x[i]=b[i]/a[i][i]; } } int main() { scanf("%d%d%d%d",&n,&m,&p,&q); double P=1.0*p/q; while(m--) { int u,v; scanf("%d%d",&u,&v); e[u][v]=e[v][u]=1,++deg[u],++deg[v]; } for(int i=1;i<=n;++i) { a[i][i]=1.0; for(int j=1;j<=n;++j) { if(e[i][j]) a[i][j]=-(1.0-P)/deg[j]; } } /* for(int i=1;i<=n;++i) { for(int j=1;j<=n;++j) { cout<<a[i][j]<<" "; } cout<<endl; }*/ b[1]=1.0; Gauss(n); for(int i=1;i<=n;++i) printf("%.9lf\n",x[i]*P); return 0; }