【CF113D】Museum【概率期望】【高斯消元】

题意：一张 $n$ 个点的无向连通图，两个人开始时分别在 $a,b$。每次在 $u$ 时会以 $p$ 的概率原地不动，$1-p$ 的概率等概率随机选择到一个相邻的点，当两人在同一点时停止。分别求在每个点相遇的概率。

$n\le 22$

网上一堆 “从起点走到 $(i,j)$ 的概率”，看得我一脸懵逼……

很容易分析出转移矩阵 $M$，然后相当于求这个东西：

$$\underset{t\to +\infty }{\lim }{M}^{t}V$$

但这个并没有通用的求法，因为矩阵的特征向量有无数多个。

算法一

比较直观的解法。

设 $f(i,j)$ 表示 $(i,j)$ 这个状态到达次数的期望，即这个状态在所有世界线中出现次数的平均值。

由于最终点只有可能出现 $0$ 次或 $1$ 次，所以它的期望次数就是概率。

而这个期望随便消一下就可以算出来。

算法二

比较本质的解法。

考虑枚举一个最终状态 $(s,s)$，在此条件下求 $f(i,j)$ 表示最终到这个状态的概率，令 $f(i,i)=[i=s]$,就可以消元了。但这样是 $O(n^7)$ 的，无法通过。

考虑一次性把每个点作为终点的 $n$ 个答案算出来，即构建出 $(n^2-n)\times (n^2)$ 的矩阵。这样有 $n^2$ 个未知数但只有 $n^2-n$ 个方程，无法解出，但可以求出 $f(a,b)$ 关于 $f(1,1),f(2,2),\dots ,f(n,n)$ 的线性表达，就可以求出答案了。

复杂度 $O(n^6)$

所以算法一算法二以及上面那个假算法写出来都一样的……

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cctype>
#include <cmath>
#include <vector>
using namespace std;
vector<int> e[25];
double p[25],a[505][505];
int n,m,sa,sb;
inline int id(int x,int y)
{
	if (x==y) return n*n-n+x;
	return (x-1)*(n-1)+y-(y>x);
}
void gauss(int n,int m)
{
	for (int i=1;i<=n;i++)
	{
		int pos=i;
		for (int j=i+1;j<=n;j++) if (fabs(a[j][i])>fabs(a[pos][i])) pos=j;
		if (pos>i) swap(a[i],a[pos]);
		for (int j=1;j<=n;j++)
			if (j!=i)
			{
				double t=a[j][i]/a[i][i];
				for (int k=i;k<=m;k++)
					a[j][k]-=t*a[i][k];
			}
	}
}
int main()
{
	scanf("%d%d%d%d",&n,&m,&sa,&sb);
	if (sa==sb)
	{
		for (int i=1;i<=n;i++) 
			if (i==sa) printf("%.10f ",1.0);
			else printf("%.10f ",0.0);
		return 0;
	}
	for (int i=1;i<=m;i++)
	{
		int u,v;
		scanf("%d%d",&u,&v);
		e[u].push_back(v),e[v].push_back(u);
	}
	for (int i=1;i<=n;i++) scanf("%lf",&p[i]);
	for (int u=1;u<=n;u++)
		for (int v=1;v<=n;v++)
			if (u!=v)
			{
				int s=id(u,v);
				double pu=1.0/e[u].size(),pv=1.0/e[v].size();
				for (int i=0;i<=(int)e[u].size();i++)
					for (int j=0;j<=(int)e[v].size();j++)
					{
						double t=(i<(int)e[u].size()? (1-p[u])*pu:p[u])*(j<(e[v].size())? (1-p[v])*pv:p[v]);
						a[s][id(i<(int)e[u].size()? e[u][i]:u,j<(int)e[v].size()? e[v][j]:v)]=t;
					}
				a[s][s]-=1;
			}
	gauss(n*n-n,n*n);
	int s=id(sa,sb);
	for (int i=n*n-n+1;i<=n*n;i++) printf("%.10f ",-a[s][i]/a[s][s]);
	return 0;
}