HDOJ 2262 Where is the canteen (高斯消元 + bfs)

先让我吐槽几句，这道题卡了我三天。主要是没有找到好的高斯消元的模板以及某些题解博客的误导。原来以为百度google到的资源是万能的，这次充分体现了团结就是力量，该抱大腿就抱大腿。情人节全献给了这道题，不过也不孤单，一直都是在讨论和头脑风暴啊中度过。希望能顺利地写完这篇题解，然后给感冒重症患者放个小假。^_^

题意：@是起点，$是终点，#不能经过，.可以经过。求从起点到终点的期望步数。注意，可能有多个终点！

从起点开始进行bfs，把能扩展到的点都做好标记。对于做标记的店，每个点都列一个方程，a[x][y] = ……+1 a表示每个点到终点的期望步数（多个终点的平均）。

只要有一个终点被扩展到，那么一定有唯一解，否则无解。

为了避免小数，减小误差，可以把所有系数都乘上一个cnt。那么原来的1/cnt变成1,1变成cnt。

为了使得操作方便，增广矩阵的系数矩阵a大小为（n*m)(n*m)，把没有扩展到的点也都列了式子，并且也有对应的未知数。 

注意系数矩阵a和常数矩阵b都是要在每一次输入后都初始化。b初始化为任意常数，因为没有扩展到的点，不会经过它到达终点，相当于就算可以经过但是到终点的期望步数为无穷大。c/0就是无穷大。

最关键的就是高斯消元部分。一般的高斯消元解方程组对于任意一个方程，不加额外条件的，可以判断这个方程是有唯一解，还是没有唯一解，没有唯一解就跳出子程序。高斯消元解方程组，一定是列出了一个方阵，化简以后，对角线以下一定都要是0，对角线上可以是0也可以是1。只有对角线上全部是1才是有唯一解的情况。否则就一定没有唯一解。后面第二个是普通高斯消元解方程组的模板。

对于这道题，这个方程组的前提是，某一些未知数确定是无解的，其他未知数一定是有唯一解的，那么在高斯消元的过程中，碰到的a[i][i]=0的情况，就是这个未知数一定是无解的意思，可以直接跳到下一行去求下一个未知数而不能跳出子程序，在回代的时候也要跳过。

点击打开链接

#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<cmath>
#include<queue>
#include<algorithm>
#define eps 1e-18
using namespace std;
struct Node{
int x;
int y;
};
const int N = 400;
const int dirx[4] = {0, 0, 1, -1};
const int diry[4] = {-1, 1, 0, 0};
int equ, var, col, max_r, x, y, cnt, sx, sy, ex, ey, n, m, u;
double tmp;
double a[N][N], ans[N], b[N];
int flag[N][N], t[N];
string ch[N];
queue<Node> q;
struct Node v, p;
void Gauss(int n)
{
    int i,j,k,r;
    for(i=0;i<n;++i)
    {
        //数据稳定性优化
        r = i;
        for(j=i+1;j<n;++j) if(fabs(a[j][i]) > fabs(a[r][i])) r=j;
        if(fabs(a[r][i]) < eps) continue;
        if(r!=i)
        {
            for(j=0;j<n;++j)swap(a[r][j],a[i][j]);
            swap(b[r],b[i]);
        }
        //消元
        for(j=n;j>=i;j--)
        for(k=i+1;k<n;++k)
            if(j==n) b[k]-=a[k][i]/a[i][i]*b[i];
            else a[k][j]-=a[k][i]/a[i][i]*a[i][j];
    }
    //回代
    for(i=n-1;i>=0;--i)
    {
        if(fabs(a[i][i]) < eps) continue;
        for(j=i+1;j<n;++j)b[i]-=a[i][j]*b[j];
        b[i]/=a[i][i];
    }
}
int bfs()
{
    memset(flag, 0, sizeof(flag));
    while (!q.empty()) q.pop();
    p.x = sx;
    p.y = sy;
    int t = 0;
    q.push(p);
    flag[sx][sy] = 1;
    while (!q.empty())
    {
        v = q.front();
        q.pop();
        for (int i = 0; i < 4; i++)
        {
            x = v.x + dirx[i];
            y = v.y + diry[i];
            if (x >= 0 && x < n && y >= 0 && y < m && !flag[x][y] && ch[x][y] != '#')
            {
                p.x = x;
                p.y = y;
                if (ch[x][y] == '$') t = 1;
                flag[x][y] = 1;
                q.push(p);
            }
        }
    }
    return t;
}
void work()
{
    memset(a, 0, sizeof(a));
     memset(b, 1, sizeof(b));
    for (int i = 0; i < n; i++)
        for (int j = 0; j < m; j++)
        if (flag[i][j] == 1)
        {
            cnt = 0;

            if (ch[i][j] == '$')
            {
                a[i*m+j][i*m+j] = 1;
                continue;
            }
            for (int k = 0; k < 4; k++)
            {
                x = i + dirx[k];
                y = j + diry[k];
                if (x >= 0 && x < n && y >= 0 && y < m && flag[x][y])
                {
                    cnt++;
                    a[i*m+j][x*m+y] = 1;
                }
            }
            b[i*m+j] = -1*cnt;
            a[i*m+j][i*m+j] = -1*cnt;
        }
}
int main()
{
    while (~scanf("%d%d", &n, &m))
    {
        memset(t, 0, sizeof(t));
        for (int i = 0; i < n; i++)
        {
            cin >> ch[i];
            for (int j = 0; j < m; j++)
            {
                if (ch[i][j] == '@')
                {
                    sx = i;
                    sy = j;
                }
            }
        }
        if (!bfs()) printf("-1\n");
        else
        {
            work();
            Gauss(n*m);
            printf("%.6f\n", b[sx*m+sy]);
        }
    }
    return 0;
}

void Gauss(int n)
{
    int i,j,k,r;
    for(i=0;i<n;++i)
    {
        //数据稳定性优化
        r = i;
        for(j=i+1;j<n;++j) if(fabs(a[j][i]) > fabs(a[r][i])) r=j;
        if(fabs(a[r][i]) < eps) return; //判断是不是唯一解
        if(r!=i)
        {
            for(j=0;j<n;++j)swap(a[r][j],a[i][j]);
            swap(b[r],b[i]);
        }
        //消元
        for(j=n;j>=i;j--)
        for(k=i+1;k<n;++k)
            if(j==n) b[k]-=a[k][i]/a[i][i]*b[i];
            else a[k][j]-=a[k][i]/a[i][i]*a[i][j];
    }
    //回代
    for(i=n-1;i>=0;--i)
    {
        for(j=i+1;j<n;++j)b[i]-=a[i][j]*b[j];
        b[i]/=a[i][i];
    }
}