CF113D Museum / BZOJ3270 博物馆

以比较抽象的视角来看， 大概就是用停留在状态 \(i\) 的概率乘上从 \(i\) 转移到 \(j\) 的概率取更新停留在状态 \(j\) 的概率，这个可以用矩阵表示。

设初始停留概率列向量为 \(st\)， 转移矩阵为 \(T\)， 那么答案向量就是：

\[\sum_{i\ge 0} T^i * st \\ = \frac{st}{I-T} \]

其中 \(I\) 是单位矩阵， \(T^0=I\)。

求矩阵逆就算了吧， 假设答案向量是 \(ans\)， 上式即为：

\[(I-T)*ans = st \]

这其实就是一个线性方程组， 高斯消元即可。

附：

\[\begin{cases} a_{1,1}*x_1+\cdots +a_{1,n}*x_1 = b_1 \\ \vdots \\ a_{n,1}*x_1+\cdots +a_{n,n}*x_1 = b_n \end{cases} \]

用矩阵表示就是：

\[\left[ \begin{matrix} a_{1,1} & \cdots & a_{1,n} \\ \vdots&\ddots&\vdots \\ a_{n,1} & \cdots &a_{n,n} \end{matrix} \right] * \left[ \begin{matrix} x_1 \\ \vdots \\ x_n \end{matrix} \right] = \left[ \begin{matrix} b_1 \\ \vdots \\ b_n \end{matrix} \right] \]

具体地， 对于转移矩阵 \(T\)， 第 \(i\) 行第 \(j\) 列的数表示从状态 \(j\) 转移到状态 \(i\) 的概率。

参考代码：

#include<bits/stdc++.h>

using namespace std;
const int N = 25, S = 503;
const double eps = 1e-9;


int n,m,a,b,g[N][N],deg[N];
double p[N];

int s;
int id(int x,int y) {	return (x-1)*n+y; 	}
double T[S][S];
void Ins(int x,int y,double k) {
	for(int i=1;i<=s+1;++i) T[x][i] += T[y][i]*k;
}

int main()
{
	scanf("%d%d%d%d",&n,&m,&a,&b);
		s = n*n;
	for(int i=0;i<m;++i)
	{
		int x,y;
		scanf("%d%d",&x,&y);
		g[x][y] = g[y][x] = 1;
		++deg[x], ++deg[y];
	}
	for(int i=1;i<=n;++i)
		scanf("%lf",&p[i]);
	for(int i=1;i<=s;++i) T[i][i]=1;
	for(int i=1;i<=n;++i)
		for(int j=1;j<=n;++j) if(i!=j)
			for(int i_=1;i_<=n;++i_) if(i==i_ || g[i][i_])
				for(int j_=1;j_<=n;++j_) if(j==j_ || g[j][j_])
					{
						int now = id(i,j), nxt = id(i_,j_);
						double P = 1;
						if(i==i_) P *= p[i];
						else P *= (1-p[i])/deg[i];
						if(j==j_) P *= p[j];
						else P *= (1-p[j])/deg[j];
						T[nxt][now] -= P;
					}
	T[id(a,b)][s+1] = 1;
	
	for(int i=1;i<=s;++i) {
		int mx = i;
		for(int j=i+1;j<=s;++j) if(fabs(T[j][i]) > fabs(T[mx][i])) mx=j;
		swap(T[i],T[mx]);
		for(int j=1;j<=s;++j) if(i^j) Ins(j,i,-T[j][i]/T[i][i]);
	}
	
	for(int i=1;i<=n;++i) {
		int q = id(i,i);
		printf("%.10lf ", T[q][s+1]/T[q][q]);
	}
	
	return 0;
}