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();//把正在处理的未知数的系数变为1for(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]);}