[ACM] hdu 2262 Where is the canteen (高斯消元求期望）

Where is the canteen

Problem Description

After a long drastic struggle with himself, LL decide to go for some snack at last. But when steping out of the dormitory, he found a serious problem : he can't remember where is the canteen... Even worse is the campus is very dark at night. So, each time he move, he check front, back, left and right to see which of those four adjacent squares are free, and randomly walk to one of the free squares until landing on a canteen.

 

Input

Each case begin with two integers n and m ( n<=15,m<=15 ), which indicate the size of the campus. Then n line follow, each contain m characters to describe the map. There are 4 different type of area in the map:

'@' is the start location. There is exactly one in each case.
'#' is an impassible square.
'$' is a canteen. There may be more than one in the campus.
'.' is a free square.

 

Output

Output the expected number of moves required to reach a canteen, which accurate to 6 fractional digits. If it is impossible , output -1.

 

Sample Input

  
  

   
   1 2
@$
2 2
@.
.$
1 3
@#$
  
  

 

Sample Output

  
  

   
   1.000000
4.000000
-1
  
  

解题思路：

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

地图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 <iostream>
#include <queue>
#include <cmath>
#include <stdio.h>
#include <iomanip>
#include <algorithm>
#include <string.h>
using namespace std;
const int maxn=250;
const double eps=1e-12;
char mp[20][20];//地图
bool flag[20][20];//判断是否已访问
double a[maxn][maxn];//每个方程中每个未知数的系数
int n,m;
int sx,sy;//起点
int dx[4]={-1,1,0,0};
int dy[4]={0,0,1,-1};

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

struct Node
{
    int x,y;
};

queue<Node>q;
Node aa,bb;

bool ok(int x,int y,int d)//判断当前坐标是否可行，广搜的时候要求flag[x][y]==0，建立方程时，要求flag[x][y]=1，后继可到达
{
    if(d==0)
    {
        if(x>=0&&x<n&&y>=0&&y<m&&mp[x][y]!='#'&&!flag[x][y])//注意别忘了flag[x][y]
            return true;
    }
    else
    {
        if(x>=0&&x<n&&y>=0&&y<m&&mp[x][y]!='#'&&flag[x][y])
            return true;
    }
    return false;
}

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

void bfs()
{
    while(!q.empty())
    {
        bb=q.front();
        q.pop();
        for(int i=0;i<4;i++)
        {
            aa.x=bb.x+dx[i];
            aa.y=bb.y+dy[i];
            if(ok(aa.x,aa.y,0))
            {
                flag[aa.x][aa.y]=1;
                q.push(aa);
            }
        }
    }
}
// 高斯消元法解方程组(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(cin>>n>>m)
    {
        while(!q.empty())
            q.pop();
        memset(flag,0,sizeof(flag));
        for(int i=0;i<n;i++)
            for(int j=0;j<m;j++)
        {
            cin>>mp[i][j];
            if(mp[i][j]=='@')
            {
                sx=i;
                sy=j;
            }
            if(mp[i][j]=='$')
            {
                aa.x=i;
                aa.y=j;
                q.push(aa);//把每个出口加入到队列中，为下面bfs做准备，注意要从出口开始广搜，不能从起点开始广搜
                flag[i][j]=1;
            }
        }
        bfs();
        memset(a,0,sizeof(a));
        for(int i=0;i<n;i++)//遍历每个格子，为每个格子建立方程，求出每个未知数的系数
            for(int j=0;j<m;j++)
        {
            int cnt=0;//计算后继有多少个状态
            if(mp[i][j]=='#')
                continue;
            if(mp[i][j]=='$')//出口
            {
                a[i*m+j][n*m]=0;
                a[i*m+j][i*m+j]=1;
                continue;
            }
            for(int s=0;s<4;s++)
            {
                int x=i+dx[s];
                int y=j+dy[s];
                if(ok(x,y,1))
                {
                    cnt++;
                    a[i*m+j][x*m+y]=1;
                }
            }
            a[i*m+j][i*m+j]=-1*cnt;
            a[i*m+j][n*m]=-1*cnt;
        }
        if(flag[sx][sy]&&gauss())
            cout<<setiosflags(ios::fixed)<<setprecision(6)<<x[sx*m+sy]<<endl;
        else
            cout<<-1<<endl;
    }
    return 0;
}