【CF113D】Museum

Portal --> cf113D

Solution

　　额题意的话大概就是给一个无向图然后两个人给两个出发点，每个点每分钟有(p[i])的概率停留，问这两个人在每个点相遇的概率是多少

　　如果说我们知道最后在哪里相遇，处理起来会比较方便。注意到(n)比较小，所以考虑枚举最后相遇的房间，然后可以看成求两个人同时从最后相遇的房间出发，走到(X)和(Y)的概率，这样初始状态什么的比较好表示，就会好求很多了（跟。。hnoi2013游走有点像？）

　　考虑(dp) ，假设我们现在枚举到终点是(t)，(f[i][j])表示第一个人在(i)这个点，第二个人在(j)这个点的概率是多少

　　那么有初始状态：

　　　　1.(i==j==t)的情况，(f[i][j]=1)

　　　　2.(i==j!=t)的情况，(f[i][j]=0)

　　考虑从第(i)个房间走向一个相邻房间的概率(k[i])，记第(i)个点的度数为(du[i])，那么有

[k[i]=frac{1-p[i]}{du[i]} ]

　　（有(1-p[i])的概率离开这个房间，并且走向每个相邻房间的概率是一样的）

　　所以我们可以得到转移：

[egin{aligned} f[i][j]&=p[i]*p[j]*f[i][j]&(都停留）\ &+k[i]*p[j]*sum f[u][j]&（第一个人走一步）\ &+k[j]*p[i]*sum f[i][v]&（第二个人走一步）\ &+k[i]*k[j]*sum f[u][v]&（两个人都走一步） end{aligned} ]

　　其中(u)满足(i)和(u)在原图中有一条边直接相连，(v)满足(v)和(j)在原图中有一条边直接相连

　　然后移下项，对于每一对((i,j))我们都可以得到一条长得像这样的式子：

[sum a_{i,j}*f[i][j]=c ]

　　将(f[i][j])看成未知数，就可以直接用高斯消元来做了

　　注意，移项的时候，因为在初始化中(f[i][i])的值已经定下来了，所以应该当成常数来看

　　

　　代码大概长这个样子

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#define db double
using namespace std;
const int N=25;
int mp[N][N],du[N];
db p[N],a[N*N][N*N],k[N],f[N],num[N][N];
int n,m,X,Y,mx;
db ans[N];
void fillmatrix(int t);
db gauss(int n);
void solve(int t);
void Fill(int x,int y);
int Id(int x,int y){return num[x][y];}

int main(){
#ifndef ONLINE_JUDGE
	freopen("a.in","r",stdin);
#endif
	int x,y;
	scanf("%d%d%d%d",&n,&m,&X,&Y);
	for (int i=1;i<=m;++i){
		scanf("%d%d",&x,&y);
		mp[x][y]=mp[y][x]=1;
		++du[x]; ++du[y];
	}
	if (X==Y){
		for (int i=1;i<=n;++i)
			if (i==X) printf("%.10lf
",1.0);
			else printf("%.10lf
",0.0);
		return 0;
	}
	for (int i=1;i<=n;++i) scanf("%lf",p+i);
	for (int i=1;i<=n;++i) k[i]=(1.0-p[i])/(1.0*du[i]);
	mx=0;
	for (int i=1;i<=n;++i)
		for (int j=1;j<=n;++j)
			if (i!=j) num[i][j]=++mx;
	for (int i=1;i<=n;++i)
		solve(i);
	for (int i=1;i<=n;++i) 
		printf("%.10lf ",ans[i]);
	printf("
");
}

void solve(int t){
	fillmatrix(t);
	ans[t]=gauss(mx);
}

void fillmatrix(int t){
	int id;
	memset(a,0,sizeof(a));
	memset(f,0,sizeof(f));
	f[t]=1;
	for (int i=1;i<=n;++i){
		for (int j=1;j<=n;++j){
			if (i!=j) Fill(i,j);
		}
	}
}

void Fill(int x,int y){
	int id=Id(x,y);
	a[id][id]=p[x]*p[y]-1.0;
	for (int i=1;i<=n;++i){
		if (mp[i][x]){
			if (i!=y)
				a[id][Id(i,y)]+=k[x]*p[y];
			else
				a[id][mx+1]-=k[x]*p[y]*f[y];
		}
		if (mp[i][y]){
			if (i!=x)
				a[id][Id(x,i)]+=k[y]*p[x];
			else
				a[id][mx+1]-=k[y]*p[x]*f[x];
		}
		for (int j=1;j<=n;++j){
			if (mp[i][x]&&mp[j][y]){
				if (i!=j)
					a[id][Id(i,j)]+=k[x]*k[y];
				else
					a[id][mx+1]-=k[x]*k[y]*f[j];
			}
		}
	}
}

db gauss(int n){
	int id;
	db tmp;
	for (int i=1;i<=n;++i){
		id=i;
		for (int j=i+1;j<=n;++j)
			if (fabs(a[j][i])>fabs(a[id][i])) id=j;
		if (id!=i)
			for (int j=1;j<=n+1;++j) swap(a[i][j],a[id][j]);	
		for (int j=i+1;j<=n;++j){
			tmp=a[j][i]/a[i][i];
			for (int k=i;k<=n+1;++k)
				a[j][k]-=tmp*a[i][k];
		}
	}
	for (int i=n;i>=1;--i){
		for (int j=n;j>i;--j)
			a[i][n+1]-=a[i][j]*a[j][n+1];
		a[i][n+1]/=a[i][i];
	}
	return a[Id(X,Y)][mx+1];
}