题意
一个包括空地、障碍和空洞的 \(H\times W\) 地图，从第一列随机选取空地作为起始位置，到达某个空洞后停止运动。给定向上下左右移动的概率比 \(p_u:p_d:p_l:p_r\)，求在每个空洞停止运动的概率为多少。
sol
由于每次到达空洞后即会停止运动，因此到达空洞的期望次数即为到达每个空洞的概率。记 \(f_{i,j}\) 表示到达 \((i,j)\) 的次数的期望，就可以很容易的推出转移方程：
\[f_{i,j}= f_{i+1,j}\cdot \frac{p_u}{psum} + f_{i-1,j}\cdot \frac{p_d}{psum} + f_{i,j+1}\cdot \frac{p_l}{psum} + f_{i,j -1}\cdot\frac{p_r}{psum} + [i=1]\cdot \frac{1}{stcnt} \]（\(stcnt\) 为第一行空地的数量，\(psum\) 为可以转移到的方向的概率比之和）
需要注意一点：当出现某一列全为 \(0\) 无法消元时，对答案并无影响，因此 continue 即可。
这样，我们只需要对列出的方程进行高斯消元即可解出答案……吗？
我们发现，这样一来，我们所求的未知数数量为 \(HW\) 个，而高斯消元的复杂度是 \(O(n^3)\) 的，因此时间复杂度就达到了 \(O((HW)^3)\)，空间复杂度达到了 \(O((HW)^2)\)，无法接受。
不过，我们的增广矩阵中，最多只会在对角线左右 \(\pm W\) 的位置出现非零数字，这是一个非常经典的优化 band-matrix。在消元和回代时，只向后修改 \(W\) 列即可，即使是矩阵中出现 \(0\)，需要换行时，也最多只需修改 \((2W)\times (2W)\) 的矩阵即可，而不需要修改整个矩阵；而且，我们可以为每一个可能出现非零数字的位置记录一个偏移量，这样时间复杂度就优化到了 \(O(HW^3)\)，空间复杂度优化到了 \(O(HW^2)\)。
注意：如果 \(\color{red}{\text{WA}}\) on #9, Line 13，说明是被卡精度了；如果结果输出了 -nan，请改变写法，不要直接使用 \(i*m+j\) 作为编号存储，而是对每一个非障碍方格进行编号。
代码
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
using namespace std;
const int N = 25, M = 10005;
const long double eps = 1e-10;
int n, m, p[4];
int dx[4] = {-1, 1, 0, 0}, dy[4] = {0, 0, -1, 1};
int pfrom[4] = {1, 0, 3, 2};
long double coeff[N * M][N * 4];
long double ans[N * M];
int psum[N * M];
int id[M][N];
char g[M][N];
int idx;
long double &c(int a, int b){
return coeff[a][b - max(0, a - 2 * m)];
}
void gauss(){
for (int i = 1; i <= idx; i ++ ){
if (abs(c(i, i)) <= eps) continue;
int ed = min(idx, i + 2 * m);
ans[i] /= c(i, i);
for (int j = ed; j >= i; j -- )
c(i, j) /= c(i, i);
for (int j = i + 1; j <= min(i + m, idx); j ++ ){
ans[j] -= c(j, i) * ans[i];
for (int k = ed; k >= i; k -- ) c(j, k) -= c(j, i) * c(i, k);
}
}
for (int i = idx; i; i -- ){
for (int j = i + 1; j <= min(idx, i + 2 * m); j ++ )
ans[i] -= c(i, j) * ans[j];
}
}
int main(){
scanf("%d%d", &m, &n);
for (int i = 0; i < 4; i ++ ) scanf("%d", &p[i]);
for (int i = 0; i < n; i ++ ) scanf("%s", g[i]);
int stcnt = 0;
for (int i = 0; i < m; i ++ ) stcnt += (g[0][i] == '.');
for (int i = 0; i < n; i ++ )
for (int j = 0; j < m; j ++ ){
id[i][j] = (g[i][j] == 'X') ? -1 : ++ idx;
if (id[i][j] == -1) continue;
for (int u = 0; u < 4; u ++ ) {
int sx = i + dx[u], sy = j + dy[u];
if (sx < 0 || sx >= n || sy < 0 || sy >= m || g[sx][sy] == 'X') continue;
psum[id[i][j]] += p[u];
}
}
for (int i = 0; i < n; i ++ )
for (int j = 0; j < m; j ++ ){
if (id[i][j] == -1) continue;
if (!i && g[i][j] == '.') ans[id[i][j]] = (long double) 1 / stcnt;
for (int u = 0; u < 4; u ++ ){
int sx = i + dx[u], sy = j + dy[u];
int px = pfrom[u];
if (sx < 0 || sx >= n || sy < 0 || sy >= m || g[sx][sy] == 'X' || g[sx][sy] == 'T') continue;
int sid = id[sx][sy];
c(id[i][j], sid) = (long double) -p[px] / psum[sid];
}
c(id[i][j], id[i][j]) = 1;
}
gauss();
for (int i = 0; i < n; i ++ )
for (int j = 0; j < m; j ++ )
if (g[i][j] == 'T') printf("%.9Lf\n", ans[id[i][j]]);
return 0;
}