【BZOJ3270】博物馆（概率dp，高斯消元）

一眼看到这题发现和【HNOI2013】游走很像，于是同样的设 dp 然后用高斯消元解：

设 \(f(a,b)\) 表示 \(A\) 在 \(a\) 房间，\(B\) 在 \(b\) 房间的期望次数，\(deg_u\) 表示 \(u\) 的度数。那么由于 \(A,B\) 相遇在某个点 \(i\) 后就不可能再往下走了（即次数只能为 \(1\)），所以 \(f(i,i)\) 就是他们相遇的概率。

容易得到状态转移方程：

\[\begin{aligned} f(a,b)=&P_aP_bf(a,b)\\ &+\sum_{(u,a)}(1-P_u)P_b\frac{f(u,b)}{deg_u}\\ &+\sum_{(v,b)}P_a(1-P_v)\frac{f(a,v)}{deg_v}\\ &+\sum_{(u,a)}\sum_{(v,b)}(1-P_u)(1-P_v)\frac{f(u ,v)}{deg_udeg_v} \end{aligned} \]

需要注意的点：

• \(f(i,i)\) 不能向其他的 \(f\) 转移，因为他们一相遇就停止继续走了。

• \(f(A,B)\) 初始化应为 \(1\)。

最后输出 \(f(i,i)\) 就好了。

代码如下：

#include<bits/stdc++.h>

#define N 25
#define M 195
#define get(a,b) ((a-1)*n+b)

using namespace std;

int n,m,A,B,deg[N];
int cnt,head[N],to[M<<1],nxt[M<<1];
double p[N],f[N*N][N*N],x[N*N];

void adde(int u,int v)
{
	to[++cnt]=v;
	nxt[cnt]=head[u];
	head[u]=cnt;
}

void Gauss()
{
	for(int i=1;i<=n*n;i++)
	{
		int p=i;
		for(int j=i+1;j<=n*n;j++)
			if(fabs(f[j][i])>fabs(f[p][i])) p=j;
		if(i!=p) swap(f[i],f[p]);
		for(int j=i+1;j<=n*n;j++)
		{
			double tmp=f[j][i]/f[i][i];
			for(int k=i;k<=n*n+1;k++) f[j][k]-=f[i][k]*tmp;
		}
	}
	for(int i=n*n;i>=1;i--)
	{
		for(int j=i+1;j<=n*n;j++) f[i][n*n+1]-=x[j]*f[i][j];
		x[i]=f[i][n*n+1]/f[i][i];
	}
}

int main()
{
	scanf("%d%d%d%d",&n,&m,&A,&B);
	for(int i=1;i<=m;i++)
	{
		int u,v;
		scanf("%d%d",&u,&v);
		adde(u,v),adde(v,u);
		deg[u]++,deg[v]++;
	}
	for(int i=1;i<=n;i++)
		scanf("%lf",&p[i]);
	f[get(A,B)][n*n+1]=1;
	for(int a=1;a<=n;a++)
	{
		for(int b=1;b<=n;b++)
		{
			int tmp=get(a,b);
			f[tmp][tmp]=1;
			if(a!=b) f[tmp][tmp]-=p[a]*p[b];
			for(int i=head[a];i;i=nxt[i])
			{
				int u=to[i];
				if(u==b) continue;
				f[tmp][get(u,b)]-=(1-p[u])*p[b]/deg[u];
			}
			for(int i=head[b];i;i=nxt[i])
			{
				int v=to[i];
				if(a==v) continue;
				f[tmp][get(a,v)]-=p[a]*(1-p[v])/deg[v];
			}
			for(int i=head[a];i;i=nxt[i])
			{
				int u=to[i];
				for(int j=head[b];j;j=nxt[j])
				{
					int v=to[j];
					if(u==v) continue;
					f[tmp][get(u,v)]-=(1-p[u])*(1-p[v])/deg[u]/deg[v];
				}
			}
		}
	}
	Gauss();
	for(int i=1;i<=n;i++)
		printf("%.6lf ",x[get(i,i)]);
	return 0;
}