题意
一个包括空地、障碍和空洞的 \(H\times W\) 地图,从第一列随机选取空地作为起始位置,到达某个空洞后停止运动。给定向上下左右移动的概率比 \(p_u:p_d:p_l:p_r\),求在每个空洞停止运动的概率为多少。
sol
由于每次到达空洞后即会停止运动,因此到达空洞的期望次数即为到达每个空洞的概率。记 \(f_{i,j}\) 表示到达 \((i,j)\) 的次数的期望,就可以很容易的推出转移方程:
(\(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;
}