挑战程序竞赛例题 4.1 Random Walk（高斯消元求期望值）

给你一幅N*M的地图，地图中有不能到达的障碍物'#'与可以走的点'.'，从（1，1）开始走到（N，M），其中每一次走动均等概率地向周围的可达的格子走去，求到达（N，M）的期望步数。（N，M<=10）

一开始根本不知道这题居然是用高斯消元来做的，感觉非常神奇，高斯消元作用就是你自己列出一系列关于期望的方程，然后求一个$E(1,1)$的变量值即可。

首先可以设每一个格子（X，Y）到达（N，M）的期望值为未知数$E(x,y)$，那么我们有N*M个格子，有N*M个未知数即N*M个变量，然后方程怎么列呢？可以发现每一个变量和周围的变量是存在关系的，由于是等概率地向周围移动，对于不可走的格子'#'，显然$E(x,y)=0$；对于任意一个可以走的格子'.'坐标为（x，y），其四个方向中有k个方向是可行的，显然有：$E(x,y) =1+ Sigma {{1 over k}{E(x',y')}|(x',y')合法}$，比如当前点为（2，3），其中从（2，3）可以到达（2，4），（3，3），（1，3），假设（2，2）是'#'而不可到达（2，2），那么有$E(2,2)={1 over 3}*E(2,4)+{1 over 3}*E(3,3)+{1 over 3}*E(1,3)+1$，两边同时乘以3再移项可以得到$3*E(2,2)-E(2,4)-E(3,3)-E(1,3)=3$，然后就可以得到很多个这样的等式，然后化成矩阵用高斯消元求解，其中答案就是E（1，1）这个点所对应的未知数的解

继昨天把j打成i死活没找到错误之后，今天又把m打成了n，又怀疑了一波人生，我可能用的是假键盘

代码：

#include <bits/stdc++.h>
using namespace std;
#define INF 0x3f3f3f3f
#define LC(x) (x<<1)
#define RC(x) ((x<<1)+1)
#define MID(x,y) ((x+y)>>1)
#define fin(name) freopen(name,"r",stdin)
#define fout(name) freopen(name,"w",stdout)
#define CLR(arr,val) memset(arr,val,sizeof(arr))
#define FAST_IO ios::sync_with_stdio(false);cin.tie(0);
typedef pair<int, int> pii;
typedef long long LL;
const double PI = acos(-1.0);
const int N = 12;
const double eps = 1e-6;

double Mat[N * N][N * N], ans[N * N];
int id[N][N];
int vis[N][N], dir[4][2] = {{1, 0}, {0, 1}, { -1, 0}, {0, -1}};
char pos[N][N];
int n, m;

void init()
{
	CLR(Mat, 0);
	CLR(vis, 0);
	CLR(ans, 0);
}
inline bool check(const int &x, const int &y)
{
	return (x >= 1 && x <= n && y >= 1 && y <= m);
}
void dfs(int x, int y)
{
	vis[x][y] = 1;
	for (int i = 0; i < 4; ++i)
	{
		int vx = x + dir[i][0];
		int vy = y + dir[i][1];
		if (check(vx, vy) && !vis[vx][vy] && pos[vx][vy] == '.')
			dfs(vx, vy);
	}
}
int Gaussian(int ne, int nv)
{
	int ce, cv, i, j;
	// for (i = 1; i <= ne; ++i)
	// 	for (j = 1; j <= nv + 1; ++j)
	// 		printf("%.0f%c", Mat[i][j], "	
"[j == nv + 1]);
	for (ce = 1, cv = 1; ce <= ne && cv <= nv; ++ce, ++cv)
	{
		int te = ce;
		for (i = ce + 1; i <= ne; ++i)
			if (fabs(Mat[i][cv]) > fabs(Mat[te][cv]))
				te = i;
		if (ce != te)
			for (j = cv; j <= nv + 1; ++j)
				swap(Mat[ce][j], Mat[te][j]);
		double bas = Mat[ce][cv];
		for (j = cv; j <= nv + 1; ++j)
			Mat[ce][j] /= bas;
		for (i = ce + 1; i <= ne; ++i)
			for (j = cv + 1; j <= nv + 1; ++j)
				Mat[i][j] -= Mat[i][cv] * Mat[ce][j];
	}
	for (i = ce; i >= 1; --i)
	{
		ans[i] = Mat[i][nv + 1];
		for (j = i + 1; j <= nv; ++j)
			ans[i] -= ans[j] * Mat[i][j];
		ans[i] /= Mat[i][i];
	}
	return 1;
}
int main(void)
{
	int i, j;
	while (~scanf("%d%d", &n, &m))
	{
		init();
		for (i = 1; i <= n; ++i)
		{
			scanf("%s", pos[i] + 1);
			for (j = 1; j <= m; ++j)
				id[i][j] = (i - 1) * m + j;
		}
		dfs(1, 1);
		for (i = 1; i <= n; ++i)
		{
			for (j = 1; j <= m; ++j)
			{
				int I = id[i][j];
				if ((i == n && j == m) || !vis[i][j])
				{
					Mat[I][I] = 1;
					continue;
				}
				int cnt = 0;
				for (int k = 0; k < 4; ++k)
				{
					int vi = i + dir[k][0];
					int vj = j + dir[k][1];
					if (check(vi, vj) && pos[vi][vj] == '.')
					{
						Mat[I][id[vi][vj]] = -1;
						++cnt;
					}
				}
				Mat[I][I] = cnt;
				Mat[I][n * m + 1] = cnt;
			}
		}
		Gaussian(n * m, n * m);
		printf("%.8f
", ans[1]);
	}
	return 0;
}
/*
10 10
..######.#
......#..#
.#.##.##.#
.#........
##.##.####
....#....#
.#######.#
....#.....
.####.####
....#.....

3 10
.#...#...#
.#.#.#.#.#
...#...#..

2 3
.#.
...
*/