First Knight UVA - 12177 （高斯消元）

题意：给出每个点向四方走的概率，保证和为1，并且向边界外走的概率为0，求从（0, 0）走到（n-1, m-1）的期望步数。

题解：dp[i][j] 表示从 （i，j）走到终点的期望步数， 然后列n * m 个方程， 用高斯消元求解。

注意精度问题，求解过程中可以构成一个下三角形式，然后直接得到在（0, 0）位置的期望步数，可以减少很多误差，并且eps需要开到1e-12。

因为这是一个稀疏矩阵，高斯消元的剪枝很重要！！

#include<iostream>
#include<sstream>
#include<iterator>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<string>
#include<set>
#include<vector>
#include<bitset>
#include<climits>
#include<queue>
#include<iomanip>
#include<cmath>
#include<stack>
#include<map>
#include<ctime>
#include<new>
using namespace std;
#define LL long long
#define ULL unsigned long long
#define MT(a,b) memset(a,b,sizeof(a))
const int INF  =  0x3f3f3f3f;
const int O    =  1e6;
const int mod  =  1e6 + 3;
const int maxn =  1605;
const double PI  =  acos(-1.0);
//const double E   =  2.718281828459;
const long double eps = 1e-12;

struct Matrix {
    int equ, var; // equ 表示等式个数，var 表示未知变量个数
    double a[maxn][maxn]; // a表示增广矩阵
    void clear() { for(int i=0; i<equ; i++) for(int j=0; j<=var; j++) a[i][j] = 0; }

    double Gauss(){
        for(int i=equ-1, col=var-1; i>=0 && col >=0; i--, col --) {
            for(int k=i-1; k>=0; k--) {
                double p = a[k][col] / a[i][col];
                if(fabs(p) < eps) continue; // 重要
                for(int j=col; j>=0; j--) a[k][j] -= a[i][j] * p;
                a[k][var] -= a[i][var] * p;
            }
        }
        return a[0][var] / a[0][0];
    }

} M;

int run[4][2]  = {1, 0, 0, 1, -1, 0, 0 , -1};
double p [maxn][maxn][4];

int main(){
    int n, m ;
    while( cin >> n >> m && n + m) {
        for (int k = 0; k < 4; k++)
            for (int i = 0; i < n; i++)
                for (int j = 0; j < m; j++) cin >> p[i][j][k];

        M.equ = M.var = n * m;
        M.clear();
        for(int i=0; i<n*m; i++) {
            int x = i / m;
            int y = i - i / m * m;
            M.a[i][i] = M.a[i][M.var] = 1.0;
            for(int k=0; k<4; k++) {
                int xx = x + run[k][0];
                int yy = y + run[k][1];
                if( xx < 0 || xx >= n || yy < 0 || yy >= m ) continue;
                int j = xx * m + yy;
                M.a[i][j] -= p[x][y][k];
            }
        }
        for(int i=0; i<=M.var; i++) if(i != M.var -1 ) M.a[M.equ-1][i] = 0;
        double ans = M.Gauss();
        printf("%.10f\n", ans);
    }
    return 0;
}