HDU 2262 Where is the canteen （高斯消元、概率）

题目连接：http://acm.hdu.edu.cn/showproblem.php?pid=2262

题意：n*m的地图，有一个起点，有多个出口，上下左右走，有的格子不能走，求从起点走到一个出口的期望步数是多少。

第一次做浮点数高斯消元求期望的题。   算是复习了一下高斯消元的知识，发现对于以前学过的东西掌握的还是不好...很多地方都想不起来了...要多复习啊。

关于本题推荐一篇博客，写得比较详细。

http://blog.csdn.net/sr_19930829/article/details/38895163

思路如下：

地图n*m看成一个个格子，编号0,1,2,3...n*m-1， 那么坐标i,j的格子编号为 i*m+j
EK表示编号为K的格子到一个出口的期望步数
那么EK=0, K这里为出口的编号
我们要求的则是EK, K这里为起点的编号
由一个状态可以走向其他状态，假设当前K可以向三个方向走，
那么 EK= (EK(a) + EK(b) +EK(c)) /3 +1
一般的 EK=(Enext1+Enext2+Enext3+......Ecnt)/ cnt+1
整理得  Enext1+Enext2+Enext3+...Ecnt - EK*cnt= - cnt
对每个EK(K=0,1,2,3...n*m-1)，都建立这样的一个等式，那么可以列出 n*m个方程，然后采用高斯消元，求出E（起点）
其中有 n*m个方程，n*m个变量
a[i][j]代表第i个式子（从0开始）,第j个未知数的系数（从0开始）,其中EK为未知数
a[0][0] *E0+a[0][1]*E1+.........a[0][n*m-1]*E(n*m-1)= a[0][n*m]
a[1][0]* E0+a[1][1]*E1+.............................................  = a[1][n*m]
.....
.....
a[n*m-1][0]*E0+a[n*m-1][1]*E1+...........................   =a[n*m-1][n*m]
所以关键要求a[i][j]，如果a数组全部求出来了，然后直接带入高斯消元模板就可以了。
预处理：从每个出口进行bfs，把能到达的位置用flag[i][j]=1标记。
遍历每个格子，对每个格子建立一个方程，求出每个方程每个未知数的系数
用高斯消元求解。
如果起点可以访问到(flag[i][j]==1且高斯消元有解），那么输出唯一解，即E（sx*m+sy）其中，sx,sy为起点的坐标
否则输出-1.

#include <bits/stdc++.h>
using namespace std;
const int maxn = 250;
const double eps = 1e-12;
char mp[20][20];
bool book[20][20];
int n, m;
int sx, sy;
int dir[4][2] = {0, 1, 1, 0, -1, 0, 0, -1};
double a[maxn][maxn];
queue<pair<int, int> > q;

int equ,var;//equ个方程,var个变量
double x[maxn];//解集
bool free_x[maxn];

int sgn(double x)  
{  
    return (x>eps)-(x<-eps);  
} 

bool check1(int x, int y) {
    if(x >= 0 && x < n && y >= 0 && y < m && book[x][y] == 0 && mp[x][y] != '#') 
        return 1;
    return 0;
}
bool check2(int x, int y) {
    if(x >= 0 && x < n && y >= 0 && y < m && book[x][y] && mp[x][y] != '#') 
        return 1;
    return 0;
}
void bfs() {
    pair<int, int> now;
    while(!q.empty()) {
        now = q.front();
        q.pop();
        for(int i = 0; i < 4; i++) {
            int tx = now.first + dir[i][0];
            int ty = now.second + dir[i][1];
            if(check1(tx, ty)) {
                q.push(make_pair(tx, ty));
                book[tx][ty] = 1;
            }
        }
    }
}

// 高斯消元法解方程组(Gauss-Jordan elimination).(0表示无解，1表示唯一解，大于1表示无穷解，并返回自由变元的个数)
int gauss()
{
    equ=n*m,var=n*m;
    int i,j,k;
    int max_r; // 当前这列绝对值最大的行.
    int col; // 当前处理的列.
    double temp;
    int free_x_num;
    int free_index;
    // 转换为阶梯阵.
    col=0; // 当前处理的列.
    memset(free_x,true,sizeof(free_x));
    for(k=0;k<equ&&col<var;k++,col++)
    {
        max_r=k;
        for(i=k+1;i<equ;i++)
        {
            if(sgn(fabs(a[i][col])-fabs(a[max_r][col]))>0)
                max_r=i;
        }
        if(max_r!=k)
        { // 与第k行交换.
            for(j=k;j<var+1;j++)
                swap(a[k][j],a[max_r][j]);
        }
        if(sgn(a[k][col])==0)
        { // 说明该col列第k行以下全是0了，则处理当前行的下一列.
            k--; continue;
        }
        for(i=k+1;i<equ;i++)
        { // 枚举要删去的行.
            if (sgn(a[i][col])!=0)
            {
                temp=a[i][col]/a[k][col];
                for(j=col;j<var+1;j++)
                {
                    a[i][j]=a[i][j]-a[k][j]*temp;
                }
            }
        }
    }

    for(i=k;i<equ;i++)
    { 
        if (sgn(a[i][col])!=0)
            return 0;
    }
    if(k<var)  //多解的情况，本题无此情况 
    {
        for(i=k-1;i>=0;i--)
        {
            free_x_num=0;
            for(j=0;j<var;j++)
            {
                if (sgn(a[i][j])!=0&&free_x[j])
                    free_x_num++,free_index=j;
            }
            if(free_x_num>1) continue;
            temp=a[i][var];
            for(j=0;j<var;j++)
            {
                if(sgn(a[i][j])!=0&&j!=free_index)
                    temp-=a[i][j]*x[j];
            }
            x[free_index]=temp/a[i][free_index];
            free_x[free_index]=0;
        }
        return var-k;
    }

    for (i=var-1;i>=0;i--)
    {
        temp=a[i][var];
        for(j=i+1;j<var;j++)
        {
            if(sgn(a[i][j])!=0)
                temp-=a[i][j]*x[j];
        }
        x[i]=temp/a[i][i];
    }
    return 1;
}


int main() {
    while(~scanf("%d %d", &n, &m)) {
        while(!q.empty()) q.pop();
        memset(book, 0, sizeof book);
        memset(a, 0, sizeof a);
        getchar();
        for(int i = 0; i < n; i++) {
            for(int j = 0; j < m; j++) {
                mp[i][j] = getchar();
                if(mp[i][j] == '@') {
                    sx = i, sy = j;
                }
                else if(mp[i][j] == '$'){
                    q.push(make_pair(i, j));
                    book[i][j] = 1;
                }
            }
            getchar();
        }
        bfs();
        for(int i = 0; i < n; i++) {
            for(int j = 0; j < m; j++) {
                if(mp[i][j] == '#') continue;
                    int cnt = 0;
                    if(mp[i][j] == '$') {
                        a[i*m+j][i*m+j] = 1;
                        a[i*m+j][m*n] = 0;
                    }
                    else {
                        for(int k = 0; k < 4; k++) {
                            int tx = i + dir[k][0];
                            int ty = j + dir[k][1];
                            if(check2(tx, ty)) {
                                a[i*m+j][tx*m+ty] = 1;
                                cnt++;
                            }
                        }
                        a[i*m+j][i*m+j] = -cnt;
                        a[i*m+j][m*n] = -cnt;
                    }
                
            }
        }
        if(book[sx][sy] && gauss()) {
            printf("%.6lf\n", x[sx*m+sy]);
        }
        else puts("-1");
    }
    return 0;
}