題目來自《挑戰程式設計競賽》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][n]=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][n];
    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;
}