有一个 \(h×w\) 的方格表，有些格子是墙，有些格子是终点，剩下的格子是空地。
你现在要在这里进行随机游走，其中，每一轮向上、下、左、右的概率分别为 \(u\%,d\%,l\%,r\%\)，如果移动后超过方格表或到达墙，则这次移动失效，保持原地不动。
当你到达终点时，游戏就结束了。现在，在方格表第一行的所有空地中均匀随机一个作为起点，求你就进入每个终点的概率。
\(1≤w≤20, 2≤h≤10000\)。
设其中一个终点为 \(t\)，\(f_{i, j}\) 表示 \((i, j)\) 走到 \(t\) 的概率。 \[ f_{i, j} = uf_{i-1, j}+df_{i+1, j}+lf_{i, j-1}+rf_{i, j+1} \]
\[ f_{i, j} - uf_{i-1, j} - df_{i+1, j} - lf_{i, j-1} - rf_{i, j+1}=0 \]
\[ f_t=1 \]
\[ P\cdot\mathbf{f}=\mathbf{e} \]
\(P\) 是 \(wh \times wh\) 的方阵，\(\mathbf{f},\mathbf{e}\) 是列向量。其中 \(e_i=[i=t]\)。 \[ \mathbf{f}=P^{-1}\mathbf{f} \] 设 \(1\dots s\) 为起点 ，答案为 \[ ans = \frac1s \sum_{i=1}^sf_i \] 写成矩阵的形式，设行向量 \(\mathbf{v}^T=\begin{pmatrix} 1 & 1 & 1 & \cdots & 0 & 0 & \cdots \end{pmatrix}\)，前面有 \(s\) 个1，其他都是0。 \[ \begin{aligned} ans &= \frac1s \mathbf{v}^T\mathbf{f}\\ &= \frac1s\mathbf{v}^TP^{-1}\mathbf{e}\\ \end{aligned} \] 相当于求行向量 \(\alpha^T = \dfrac1s\mathbf{v}^TP^{-1}\) 第 \(t\) 维的值，所以只要求出这个行向量即可。 \[ \alpha^T P = \frac1s\mathbf{v} \]
\[ P^T\alpha=\frac1s\mathbf{v}^T \]
解这个线性方程组，直接高斯消元显然不行。
如果按照从上到下从左到右的顺序标号，\(P^T\) 是一个带状矩阵，每行只需存储 \(2w+1\) 个元素。
然后消元即可。
这还有个小问题没解决，留坑！
// Author: HolyK
// Created: Mon Jul 5 11:36:58 2021
#include <bits/stdc++.h>
#define dbg(a...) fprintf(stderr, a)
template <class T, class U>
inline bool smin(T &x, const U &y) {
return y < x ? x = y, 1 : 0;
}
template <class T, class U>
inline bool smax(T &x, const U &y) {
return x < y ? x = y, 1 : 0;
}
using LL = long long;
using PII = std::pair<int, int>;
constexpr int dx[] = {-1, 1, 0, 0}, dy[] = {0, 0, -1, 1};
int main() {
std::ios::sync_with_stdio(false);
std::cin.tie(nullptr);
int n, m, tot;
double p[4];
std::cin >> m >> n;
tot = n * m;
for (int i = 0; i < 4; i++) std::cin >> p[i], p[i] /= 100;
std::vector<std::string> s(n);
for (int i = 0; i < n; i++) std::cin >> s[i];
std::vector<std::vector<double>> a(tot, std::vector<double>(m * 2 + 1));
for (int i = 0, o = 0; i < n; i++) {
for (int j = 0; j < m; j++, o++) {
a[o][m] = 1;
if (s[i][j] != '.') continue;
for (int d = 0; d < 4; d++) {
int x = i + dx[d], y = j + dy[d];
if (x < 0 || y < 0 || x >= n || y >= m || s[x][y] == 'X') {
a[o][m] -= p[d];
} else {
// a[o][m + dx[d] * m + dy[d]] -= p[d ^ 1];
a[x * m + y][m + o - (x * m + y)] -= p[d];
}
}
}
}
// for (int i = 0; i < tot; i++) {
// for (int j = 0; j < m * 2 + 1; j++) std::cout << a[i][j] << " \n"[j == 2 * m];
// }
int cnt = 0;
for (int i = 0; i < m; i++) cnt += s[0][i] == '.';
std::vector<double> b(tot);
for (int i = 0; i < m; i++) b[i] = 1.0 / cnt;
for (int i = 0; i < tot; i++) {
if (fabs(a[i][m]) < 1e-12) {
b[i] = 0;
continue;
}
for (int j = i + 1; j < tot && j <= i + m; j++) {
double c = a[j][m + i - j] / a[i][m];
for (int k = 0; k <= m; k++) a[j][m + k + i - j] -= a[i][m + k] * c;
b[j] -= b[i] * c;
}
}
for (int i = tot - 1; i >= 0; i--) {
if (fabs(a[i][m]) < 1e-12) continue;
for (int j = i + 1; j < tot && j <= i + m; j++) {
b[i] -= b[j] * a[i][m + j - i];
}
b[i] /= a[i][m];
}
for (int i = 0; i < n; i++)
for (int j = 0; j < m; j++)
if (s[i][j] == 'T')
printf("%.9f\n", b[i * m + j]);
return 0;
}最后修改于 2021-08-13