首页 > 代码库 > [ACM] hdu 2262 Where is the canteen (高斯消元求期望)
[ACM] hdu 2262 Where is the canteen (高斯消元求期望)
Where is the canteen
‘@‘ 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.
1 2 @$ 2 2 @. .$ 1 3 @#$
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; }
[ACM] hdu 2262 Where is the canteen (高斯消元求期望)