random walk

题意：给定一个N*M矩阵，'.'表示可以通过，'#'表示不能通过，求从[0][0]走到[N-1][M-1]的期望步数

递推式写出来，再转化为矩阵，然后高斯消元求解

#include<stdio.h>
#include<iostream>
#include<vector>
using namespace std;
const double EPS = 1e-8;
typedef vector<double> vec;
typedef vector<vec> mat;
const int MAX_N=10;
const int MAX_M=10;
//输入

char grid[MAX_N][MAX_M+1];
int N,M;
bool can_goal[MAX_N][MAX_M];	//can_goal[x][y]为true的话，(x,y)可以到达终点
int dx[4]={-1,1,0,0},dy[4]={0,0,-1,1};

//求解Ax-b,其中A是方阵
//当方程组无解或者有无穷多解时，返回一个长度为0的数组
vec gauss_jordan(const mat& A,const vec& b) {
	int i,j,k;
	int n=A.size();
	mat B(n,vec(n+1));
	for(i=0;i<n;i++) 
		for(j=0;j<n;j++) B[i][j]=A[i][j];
	//把b存放在A的右边方便一起处理
	for(i=0;i<n;i++) B[i][n]=b[i];

	for(i=0;i<n;i++) {
		//把正在处理的未知数系数的绝对值最大的式子换到第i行
		int pivot=i;
		for(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(j=i+1;j<=n;j++) B[i][j]/=B[i][i];
		for(j=0;j<n;j++) {
			if (i!=j) {
				//从第j个式子中消去第i个未知数
				for(k=i+1;k<=n;k++) B[j][k]-=B[j][i]*B[i][k];
			}
		}
	}
	vec x(n);
	//存放在右边的b就是答案
	for(i=0;i<n;i++) x[i]=B[i][n];
	return x;
}

//搜索可以达到终点的点
void dfs(int x,int y) {
	can_goal[x][y]=true;
	int i;
	for(i=0;i<4;i++) {
		int nx=x+dx[i],ny=y+dy[i];
		if (0<=nx && nx < N && 0<=ny && ny< M && !can_goal[nx][ny] && grid[nx][ny]!='#') {
			dfs(nx,ny);//这里grid[nx][ny]可能会数组越界
		}
	}
}

void solve() {
	mat A(N*M,vec(N*M,0));
	vec b(N*M,0);
	int x,y,k;
	for(x=0;x<N;x++) {
		for(y=0;y<M;y++) {
			can_goal[x][y]=false;
		}
	}
	dfs(N-1,M-1);	//从终点开始搜，所有可达的点就是能到终点的点
/*
E(x,y)=1/4E(x-1,y)+1/4E(x+1,y)+1/4E(x,y-1)+1/4E(x,y+1)+1
转化为矩阵的形式:
4*E(x,y)-E(x-1,y)-E(x+1,y)-E(x,y-1)-E(x,y+1)=4
*/
	//构建矩阵
	for(x=0;x<N;x++) {
		for(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(k=0;k<4;k++) {
				int nx=x+dx[k],ny=y+dy[k];
				if (0<=nx && nx<N && 0<=ny && 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\n",res[0]);
}