SRM 356 950pt : EscapeTheJail

keyword

連立一次方程式　期待値　C++

問題概要

N*M(N,M<15)のボードがあり、スタート地点、ゴール地点（複数）、壁が指定されている。到達可能なマスを等確率で選んで移動するとき、ゴールにたどり着くまでのターン数の期待値を求める問題。

解法

蟻本に載ってるランダムウォークの問題とほとんど同じ。ゴールが複数あるという違いくらいしかない。

#include <vector>
#include <string>
#include <cmath>
using namespace std;

int dy[] = {1,-1,0,0};
int dx[] = {0,0,1,-1};
vector<string> jail;
int N, M, S;

const double EPS = 1e-10;

double A[15*15][15*15];
double b[15*15];

struct EscapeTheJail {

	double findExit(vector <string> _jail) {
		jail = _jail;
		N = jail.size(); M = jail[0].size();
		S = N*M;
		int sx=-1, sy=-1;
		for(int i=0; i<N; i++)for(int j=0; j<M; j++){
			if(jail[i][j] == '@'){
				sy = i; sx = j;
				jail[i][j] = '.';
			}
			if(jail[i][j] == '$'){
				jail[i][j] = '.';
			}
		}

		dfs(sy, sx);
		bool found = false;
		for(int i=0; i<N; i++)for(int j=0; j<M; j++){
			if(_jail[i][j] == '$' && jail[i][j] == '*'){
				found = true;
			}
		}
		if(!found) return -1.0;

		for(int i=0; i<N; i++)for(int j=0; j<M; j++){
			if(jail[i][j] != '*') jail[i][j] = '#';
		}

		for(int i=0; i<N; i++)for(int j=0; j<M; j++){
			if(_jail[i][j] == '$'){
				A[i*M + j][i*M + j] = 1.0;
				b[i*M + j] = 0.0;
				continue;
			}
			int cnt = 0;
			b[i*M + j] = 1.0;
			A[i*M + j][i*M + j] = 1.0;
			for(int k=0; k<4; k++){
				int ny = i + dy[k], nx = j + dx[k];
				if(0 <= ny && ny < N && 0 <= nx && nx < M && jail[ny][nx] == '*'){
					cnt++;
				}
			}
			for(int k=0; k<4; k++){
				int ny = i + dy[k], nx = j + dx[k];
				if(0 <= ny && ny < N && 0 <= nx && nx < M && jail[ny][nx] == '*'){
					A[i*M + j][ny*M + nx] = -1.0/cnt;
				}
			}
		}

		gauss_jordan();
		return b[sy*M + sx];
	}

	void dfs(int y, int x){
		jail[y][x] = '*';
		for(int k=0; k<4; k++){
			int ny = y + dy[k], nx = x + dx[k];
			if(0 <= ny && ny < N && 0 <= nx && nx < M && jail[ny][nx] == '.'){
				dfs(ny,nx);
			}
		}
	}

	void gauss_jordan(){
		for(int i=0; i<S; i++){
			int pivot = i;
			for(int j=i; j<S; j++){
				if(fabs(A[j][i]) > fabs(A[pivot][i])) pivot = j;
			}
			for(int j=0; j<S; j++){
				swap(A[i][j], A[pivot][j]);
			}
			swap(b[i], b[pivot]);

			if(fabs(A[i][i]) < EPS){
				for(int j=0; j<S; j++){
					b[j] = -1;
				}
				return ;
			}

			for(int j=i+1; j<S; j++) A[i][j] /= A[i][i]; b[i] /= A[i][i];
			for(int j=0; j<S; j++){
				if(i != j){
					for(int k=i+1; k<S; k++){
						A[j][k] -= A[j][i]*A[i][k];
					}
					b[j] -= A[j][i] * b[i];
				}
			}
		}
	}

};