Random Walk 挑战程序设计竞赛 期望值和方程组

题目来自《挑战程序设计竞赛》4.1更加复杂的数学问题 期望值和方程组

这题其实就是ZJUT 1423，然而ZJUT似乎挂了。。。

1.题目详情

有一个N*M的格子，从(0,0)出发，每一步朝着上下左右四个格子中可以移动的格子等概率的移动。另外有些格子中有石头，因此无法移至这些格子。求第一次到达(N-1,M-1)格子的期望步数。题目假定至少存在一条从(0,0)出发到格子(N-1,M-1)的路径

限制条件

2<=N,M<=10

样例：('#'和'.'分别表示石头和可以移动的格子)

2.解题思路

设E(x,y)表示从格子(x,y)出发，到达终点的期望步数。先考虑从格子(x,y)向上下左右四个方向移动的情况，由于是等概率的，有如下关系:E(x,y)=1/4*E(x-1,y)+1/4*E(x,y+1)+1/4*E(x,y-1)+1/4*E(x+1,y)+1。如果移动不是等概率的，把1/4改成相应的数值就可以了。如果存在不能移动方向，也可以列出类似的式子。此外，当(x,y)=(N-1,M-1)时，有E(N-1,M-1)=0.为了使方程有唯一解，我们令无法到达终点的格子和有石头的格子都有E(x,y)=0。把N*M个方程联立起来就可以求期望步数了。

另外，为了更方便的表示上述状态转移方程，把上式改写成4*E(x,y)-E(x-1,y)-E(x,y+1)-E(x,y-1)-E(x+1,y)=4(把方程整数化)，其中4表示格子(x,y)可以移动的方向的数目(可以上下左右移动，这就很容易理解，代码中的move以及为什么令A[x*M+y][nx*M+ny]=-1了。

3.代码

#include <iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<vector>
using namespace std;
const double EPS= 1e-8 ;
typedef vector<double> vec;
typedef vector<vec> mat;
char grid[15][15];
int N,M;
bool can_goal[15][15];//can_goal[x][y]为true，表示格子(x,y)可以到达终点
int dx[]={-1,0,0,1};
int dy[]={0,1,-1,0};
//求解Ax=b，其中A是方阵
//当方程组无解或有无穷多解时，返回一个长度为0的数组
vec gauss_jordan(const mat& A,const vec& b){
int n=A.size();
mat B(n,vec(n+1));
//把A复制给B
for (int i=0;i<n;i++)
for (int j=0;j<n;j++) B[i][j]=A[i][j];
//把正在处理的未知数系数的绝对值最大的式子换到第i行
for (int i=0;i<n;i++) B[i]
=b[i];
for (int i=0;i<n;i++){
int pivot=i;
for (int j=i;j<n;j++){
if (abs(B[j][i])>abs(B[pivot][i])) pivot=j;
}
swap(B[i],B[pivot]);

//无解或有无穷多解
if (abs(B[i][i])<EPS) return vec();

//把正在处理的未知数系数变为1
for (int j=i+1;j<=n;j++) B[i][j]/=B[i][i];
for (int j=0;j<n;j++){
if (i!=j){
//从第j个式子中消去第i个未知数
for (int k=i+1;k<=n;k++) B[j][k]-=B[j][i]*B[i][k];
}
}
}
vec x(n);
//存放在右边的b就是答案
for (int i=0;i<n;i++) x[i]=B[i]
;
return x;
}
//搜索可以到达的点
void dfs(int x,int y)
{
can_goal[x][y]=true;
for(int i=0;i<4;i++){
int nx=x+dx[i],ny=y+dy[i];
if(nx>=0&&nx<N&&ny>=0&&ny<M&&!can_goal[nx][ny]&&grid[nx][ny]!='#'){
dfs(nx,ny);
}
}
return;
}
void solve()
{
//全部置为0
mat A(N*M,vec(N*M,0));
vec b(N*M,0);
for(int x=0;x<N;x++){
for(int y=0;y<M;y++){
can_goal[x][y]=false;
}
}
dfs(N-1,M-1);
//构建矩阵
for(int x=0;x<N;x++){
for(int y=0;y<M;y++){
//到达终点或者(x,y)是无法到达终点的情况
if((x==N-1&&y==M-1)||!can_goal[x][y]){
A[x*M+y][x*M+y]=1;
continue;
}
//其余情况
int move=0;
for(int k=0;k<4;k++){
int nx=x+dx[k],ny=y+dy[k];
if(nx>=0&&nx<N&&nx>=0&&ny<M&&grid[nx][ny]=='.'){
A[x*M+y][nx*M+ny]=-1;
move++;
}
}
b[x*M+y]=A[x*M+y][x*M+y]=move;
}
}
vec res=gauss_jordan(A,b);
printf("%.8f",res[0]);
return;
}
int main()
{
while(scanf("%d%d",&N,&M)!=EOF){
for(int i=0;i<N;i++){
for(int j=0;j<M;j++){
cin>>grid[i][j];
}
}
solve();
}
return 0;
}