高斯消元

高斯消元

高斯消元法通常用于求解如下的 \(n\) 元线性方程组

\[\begin{cases}a_{1, 1} x_1 + a_{2, 2} x_2 + \cdots + a_{1, n} x_n = b_1\\a_{2, 1} x_1 + a_{2, 2} x_2 + \cdots + a_{2, n} x_n = b_2\\\cdots\\a_{n, 1} x_1 + a_{n, 2} x_2 + \cdots + a_{n, n} x_n = b_n\end{cases} \]

高斯消元解线性方程组

考虑求解如下方程组：

\[\begin{cases}x_1 + 3x_2 + 4x_3 = 5\\x_1 + 4x_2 + 7x_3 = 3\\9x_1 + 3x_2 + 2x_3 = 2\end{cases} \]

把系数写到矩阵里：

\[\begin{bmatrix} 1 & 3 & 4 & 5 \\ 1 & 4 & 7 & 3 \\ 9 & 3 & 2 & 2 \\ \end{bmatrix} \]

把未知项的系数单独拿出来考虑：

\[\begin{bmatrix} 1 & 3 & 4 \\ 1 & 4 & 7 \\ 9 & 3 & 2 \\ \end{bmatrix} \]

如果能够将其化为：

\[\begin{bmatrix} ? & ? & ? \\ 0 & ? & ? \\ 0 & 0 & ? \\ \end{bmatrix} \]

就可以从最后一行逆推上去。设第一行第一列为主元，利用加减消元法配系数消掉后两行的第一列：

\[\begin{bmatrix} 1 & 3 & 4 \\ 0 & 1 & 3 \\ 0 & -24 & -34 \\ \end{bmatrix} \]

再以第二行第二列为主元，加减消元得到：

\[\begin{bmatrix} 1 & 3 & 4 \\ 0 & 1 & 3 \\ 0 & 0 & 38 \\ \end{bmatrix} \]

即得目标状态，实现时把 \(b\) 的那一列一起带入消元即可。

每次选取系数最大的为主元可以有效降低精度误差。

P3389 【模板】高斯消元法

#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-9;
const int N = 1e2 + 7;double g[N][N], ans[N];int n;inline bool Gauss() {for (int i = 1; i <= n; ++i) {int mxpos = i;for (int j = i + 1; j <= n; ++j)if (fabs(g[j][i]) > fabs(g[mxpos][i]))mxpos = j;swap(g[i], g[mxpos]);if (fabs(g[i][i]) < eps)return false;for (int j = i + 1; j <= n; ++j) {double div = g[j][i] / g[i][i];for (int k = i + 1; k <= n + 1; ++k)g[j][k] -= g[i][k] * div;}}return true;
}signed main() {scanf("%d", &n);for (int i = 1; i <= n; ++i)for (int j = 1; j <= n + 1; ++j)scanf("%lf", g[i] + j);if (!Gauss())return puts("No Solution"), 0;for (int i = n; i; --i) {ans[i] = g[i][n + 1];for (int j = i + 1; j <= n; ++j)ans[i] -= g[i][j] * ans[j];ans[i] /= g[i][i];}for (int i = 1; i <= n; ++i)printf("%.2lf\n", ans[i]);return 0;
}

约旦消元法解线性方程组

其为高斯消元的进阶版，能避免回带来求出答案，也就是要把矩阵变为：

\[\begin{bmatrix} ? & 0 & 0 & ? \\ 0 & ? & 0 & ? \\ 0 & 0 & ? & ? \\ \end{bmatrix} \]

首先，我们依照朴素的高斯消元不难得到：

\[\begin{bmatrix} ? & ? & ? & ? \\ 0 & ? & ? & ? \\ 0 & 0 & ? & ? \\ \end{bmatrix} \]

观察上下两个矩阵，不难得到，我们在消元时不仅仅消去主元所在式子下方的式子，而对于上方的式子也应当予以处理。

注意一下特殊情况：

• 无解：当左边系数为 \(0\) 而右边系数不为 \(0\) 时即为无解。

• 多解：当左边系数为 \(0\) 而右边系数也为 \(0\) 时即为多解。

inline bool Gauss() {for (int i = 1; i <= n; ++i) {int mxpos = i;for (int j = i + 1; j <= n; ++j)if (fabs(g[j][i]) > fabs(g[mxpos][i]))mxpos = j;if (i != mxpos)swap(g[i], g[mxpos]);if (fabs(g[i][i]) < eps)return false;for (int j = 1; j <= n; ++j)if (j != i) {double div = g[j][i] / g[i][i];for (int k = i + 1; k <= n + 1; ++k)g[j][k] -= g[i][k] * div;}}for (int i = 1; i <= n; ++i)ans[i] = g[i][n + 1] / g[i][i];return true;
}

约旦消元法解异或方程组

异或方程组是指形如：

\[\begin{cases} a_{1, 1} x_1 \oplus a_{2, 2} x_2 \oplus \cdots \oplus a_{1, n} x_n = b_1 \\ a_{2, 1} x_1 \oplus a_{2, 2} x_2 \oplus \cdots \oplus a_{2, n} x_n = b_2 \\ \cdots \\ a_{n, 1} x_1 \oplus a_{n, 2} x_2 \oplus \cdots \oplus a_{n, n} x_n = b_n \end{cases} \]

的方程组，其中 \(a_{i, j}, b_i\) 的值为 \(0\) 或 \(1\) 。

由于 \(\oplus\) 符合交换律与结合律，故可以按照高斯消元法逐步消元求解。值得注意的是，我们在消元的时候应使用异或消元而非加减消元，且不需要进行乘除改变系数（因为系数均为 \(0\) 和 \(1\) ）。

注意到异或方程组的增广矩阵是 \(01\) 矩阵，可以用 bitset 优化到 \(O(\dfrac{n^2 m}{\omega})\) 。

inline bool Gauss() {for (int i = 1; i <= n; ++i) {int mxpos = i;while (mxpos <= n && !g[mxpos].test(i))++mxpos;if (mxpos > n)return false;if (i != mxpos)swap(g[i], g[mxpos]);for (int j = 1; j <= m; ++j)if (j != i && g[j].test(i))g[j] ^= g[i];}return true;
}

矩阵求逆

对于矩阵 \(A\) ，若存在矩阵 \(A^{-1}\) 使得 \(A \times A^{-1} = A^{-1} \times A = I\) ，则称矩阵 \(A\) 可逆，\(A^{-1}\) 为其逆矩阵。

给出 \(n\) 阶方阵 \(A\) ，求解其逆矩阵的方法如下：

• 构造一个 \(n \times 2n\) 的矩阵 \((A, I_n)\)
• 用高斯消元法将其化简为最简形 \((I_n, A^{-1})\) ，即可得到 \(A\) 的逆矩阵 \(A^{-1}\) 。
• 如果最终最简形的左半部分不是单位矩阵 \(I_n\) ，则矩阵 \(A\) 不可逆。

inline bool Gauss() {for (int i = 1; i <= n; ++i)g[i][i + n] = 1;for (int i = 1; i <= n; ++i) {int mxpos = i;while (mxpos <= n && !g[mxpos][i])++mxpos;if (mxpos > n)return false;if (i != mxpos)swap(g[i], g[mxpos]);int inv = mi(g[i][i], Mod - 2);for (int j = 1; j <= n; ++j) {if (j == i)continue;int div = 1ll * g[j][i] * inv % Mod;for (int k = i; k <= n * 2; ++k)g[j][k] = dec(g[j][k], 1ll * g[i][k] * div % Mod);}}for (int i = 1; i <= n; ++i) {int inv = mi(g[i][i], Mod - 2);for (int j = 1; j <= 2 * n; ++j)g[i][j] = 1ll * g[i][j] * inv % Mod;}return true;
}

band-matrix

长这个样子：

空白部分都为 \(0\) ，橙色部分可以为任何数，这样中间就形成了一个宽度大约为 \(d\) 的带。

可以发现任意一个 \(i\) 满足从 \((i, i)\) 向右或向下拓展都有不超过 \(d - 1\) 个非零数字，即很多位置根本不需要消。

具体地，假设现在要消第 \(i\) 列，那么从第 \(i\) 行开始往下枚举 \(d - 1\) 行，每行往右消 \(d\) 个数字即可，最后仍能得到一个上三角矩阵。

与普通高斯消元有点不一样的地方在于当主元为 \(0\) 的时候的处理方法。在 band-matrix 中，若直接交换行会破坏 band-matrix 。注意到每次交换完后交换的行右边最多多出 \(d\) 个数，于是每次往右消元 \(2d\) 个数即可。

时间复杂度 \(O(nd^2)\) 。

inline bool Gauss(int n, int band) {for (int i = 1; i <= n; ++i) {int mxpos = i;while (mxpos <= min(i + band, n) && fabs(g[mxpos][i]) < eps)++mxpos;if (mxpos > min(i + band, n))return false;if (i != mxpos)swap(g[i], g[mxpos]);for (int j = i + 1; j <= min(i + band, n); ++j) {double div = g[j][i] / g[i][i];for (int k = i; k <= min(i + 2 * band, n); ++k)g[j][k] -= div * g[i][k];g[j][n + 1] -= div * g[i][n + 1];}}for (int i = n; i; --i) {ans[i] = g[i][n + 1];for (int j = i + 1; j <= min(i + 2 * band, n); ++j)ans[i] -= g[i][j] * ans[j];ans[i] /= g[i][i];}return true;
}

还有一种解决主元为 \(0\) 的方法，普通高消是交换行，这里只要交换列就可以保持 band-matrix 的性质了，时间复杂度也是 \(O(nd^2)\) 。

inline void Gauss(int n, int band) {iota(id + 1, id + 1 + n, 1);for (int i = 1; i <= n; ++i) {if (!g[i][i]) {for (int j = i + 1; j <= min(n, i + band); ++j)if (g[i][j]) {for (int k = 1; k <= n; ++k)swap(g[k][i], g[k][j]);swap(id[i], id[j]);break;}}int Inv = mi(g[i][i], Mod - 2);for (int j = i + 1; j <= n; ++j) {int div = 1ll * g[j][i] * Inv % Mod;for (int k = i; k <= min(n, i + band); ++k)g[j][k] = dec(g[j][k], 1ll * g[i][k] * div % Mod);g[j][n + 1] = dec(g[j][n + 1], 1ll * g[i][n + 1] * div % Mod);}}for (int i = n; i; --i) {ans[id[i]] = g[i][n + 1];for (int j = i + 1; j <= min(n, i + band); ++j)ans[id[i]] = dec(ans[i], 1ll * ans[id[j]] * g[i][j] % Mod);ans[id[i]] = 1ll * ans[id[i]] * mi(g[i][i], Mod - 2) % Mod;}
}

应用

Broken robot

\(n\) 行 \(m\) 列的矩阵，现在在 \((x,y)\)，每次等概率向左、右、下走或原地不动，但不能走出去，求走到最后一行期望的步数。

\(n, m \leq 10^3\)

记 \(f_{i, j}\) 表示机器人在 \((i, j)\) 时走到最后一行的期望步数，则：

\(m = 1\) 时有（省略第二维）：

\[f_i = 1 + \dfrac{1}{2} (f_{i + 1} + f_i) \]

即：

\[f_i = f_{i + 1} + 2 \]

\(m > 1\) 时有：

\[f_{i, j} = 1 +\begin{cases}\dfrac{1}{3} (f_{i, j} + f_{i + 1, j} + f_{i, j + 1}), & j = 1 \\\dfrac{1}{4} (f_{i, j} + f_{i + 1, j} + f_{i, j - 1} + f_{i, j + 1}), & 1 < j < m \\\dfrac{1}{3} (f_{i, j} + f_{i + 1, j} + f_{i, j - 1}), & j = m\end{cases} \]

注意到这是一个 \(d = 2\) 的 band-matrix ，直接使用 band-matrix 消元即可，时间复杂度 \(O(nmd^2)\) 。

#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-12;
const int N = 1e3 + 7;double g[N][N], ans[N][N];int n, m, x, y;inline void Gauss(int n, int band, double *ans) {for (int i = 1; i <= n; ++i) {int mxpos = i;while (mxpos <= min(i + band, n) && fabs(g[mxpos][i]) < eps)++mxpos;if (mxpos > min(i + band, n))continue;if (i != mxpos)swap(g[i], g[mxpos]);for (int j = i + 1; j <= min(i + band, n); ++j) {double div = g[j][i] / g[i][i];for (int k = i; k <= min(i + 2 * band, n); ++k)g[j][k] -= div * g[i][k];g[j][n + 1] -= div * g[i][n + 1];}}for (int i = n; i; --i) {ans[i] = g[i][n + 1];for (int j = i + 1; j <= min(i + 2 * band, n); ++j)ans[i] -= g[i][j] * ans[j];ans[i] /= g[i][i];}
}signed main() {scanf("%d%d%d%d", &n, &m, &x, &y);if (m == 1)return printf("%.10lf\n", 2.0 * (n - x)), 0;for (int i = n - 1; i; --i) {g[1][1] = 2, g[1][2] = -1, g[1][m + 1] = 3 + ans[i + 1][1];for (int j = 2; j < m; ++j)g[j][j] = 3, g[j][j - 1] = g[j][j + 1] = -1, g[j][m + 1] = 4 + ans[i + 1][j];g[m][m] = 2, g[m][m - 1] = -1, g[m][m + 1] = 3 + ans[i + 1][m];Gauss(m, 1, ans[i]);}printf("%.10lf\n", ans[x][y]);return 0;
}

Circles of Waiting

平面直角坐标系上有一个点，一开始在 \((0, 0)\) ，每秒钟这个点都会随机移动，如果它在 \((x, y)\) ，下一秒：

• 在 \((x - 1, y)\) 的概率是 \(p_1\)
• 在 \((x, y - 1)\) 的概率是 \(p_2\)
• 在 \((x + 1, y)\) 的概率是 \(p_3\)
• 在 \((x, y + 1)\) 的概率是 \(p_4\)

保证 \(p_1 + p_2 + p_3 + p_4 = 1\) ，求该点移动至距离原点距离为大于 \(R\) 的点的期望步数

\(0 \leq R \leq 50\)

把所有满足 \(i^2 + j^2 \leq R^2\) 的点依次编号，显然有 \(O(R^2)\) 个点。

设 \(f_{id(i, j)}\) 表示 \((i, j)\) 走出圆的期望步数，\((i, j)\) 只要转移到 \((i, j - 1), (i - 1, j), (i, j + 1), (i + 1, j)\) 。因为是依次编号，所以建出来的矩阵带宽 \(\leq 2R + 1\) 。套用 band-matrix ，可以做到 \(O(R^2 d^2) = O(R^4)\) 。

#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int Mod = 1e9 + 7;
const int R = 5e1 + 7, SIZE = 8e3 + 7;vector<pair<int, int> > Pos;int g[SIZE][SIZE];
int id[R << 1][R << 1];
int ans[SIZE];int r, p1, p2, p3, p4;template <class T = int>
inline T read() {char c = getchar();bool sign = c == '-';while (c < '0' || c > '9')c = getchar(), sign |= c == '-';T x = 0;while ('0' <= c && c <= '9')x = (x << 1) + (x << 3) + (c & 15), c = getchar();return sign ? (~x + 1) : x;
}inline int add(int x, int y) {x += y;if (x >= Mod)x -= Mod;return x;
}inline int dec(int x, int y) {x -= y;if (x < 0)x += Mod;return x;
}inline int mi(int g, int b) {int res = 1;for (; b; b >>= 1, g = 1ll * g * g % Mod)if (b & 1)res = 1ll * res * g % Mod;return res;
}inline int getid(int x, int y) {return id[x + R][y + R];
}inline void Gauss(int n, int band) {for (int i = 1; i <= n; ++i) {int mxpos = i;while (mxpos <= min(i + band, n) && !g[mxpos][i])++mxpos;if (mxpos > min(i + band, n))continue;if (i != mxpos)swap(g[i], g[mxpos]);int inv = mi(g[i][i], Mod - 2);for (int j = i + 1; j <= min(i + band, n); ++j) {int div = 1ll * g[j][i] * inv % Mod;for (int k = i; k <= min(i + 2 * band, n); ++k)g[j][k] = dec(g[j][k], 1ll * div * g[i][k] % Mod);g[j][n + 1] = dec(g[j][n + 1], 1ll * div * g[i][n + 1] % Mod);}}for (int i = n; i; --i) {ans[i] = g[i][n + 1];for (int j = i + 1; j <= min(i + 2 * band, n); ++j)ans[i] = dec(ans[i], 1ll * g[i][j] * ans[j] % Mod);ans[i] = 1ll * ans[i] * mi(g[i][i], Mod - 2) % Mod;}
}signed main() {r = read(), p1 = read(), p2 = read(), p3 = read(), p4 = read();int invsum = mi(add(add(p1, p2), add(p3, p4)), Mod - 2);p1 = 1ll * p1 * invsum % Mod, p2 = 1ll * p2 * invsum % Mod;p3 = 1ll * p3 * invsum % Mod, p4 = 1ll * p4 * invsum % Mod;for (int i = -r; i <= r; ++i)for (int j = -r; j <= r; ++j)if (i * i + j * j <= r * r)Pos.emplace_back(make_pair(i, j)), id[i + R][j + R] = Pos.size();int n = Pos.size(), band = 0;for (auto it : Pos) {int x = it.first, y = it.second;g[getid(x, y)][getid(x, y)] = g[getid(x, y)][n + 1] = 1;if (getid(x - 1, y)) {g[getid(x, y)][getid(x - 1, y)] = Mod - p1;band = max(band, abs(getid(x, y) - getid(x - 1, y)));}if (getid(x, y - 1)) {g[getid(x, y)][getid(x, y - 1)] = Mod - p2;band = max(band, abs(getid(x, y) - getid(x, y - 1)));}if (getid(x + 1, y)) {g[getid(x, y)][getid(x + 1, y)] = Mod - p3;band = max(band, abs(getid(x, y) - getid(x + 1, y)));}if (getid(x, y + 1)) {g[getid(x, y)][getid(x, y + 1)] = Mod - p4;band = max(band, abs(getid(x, y) - getid(x, y + 1)));}}Gauss(n, band);printf("%d", ans[getid(0, 0)]);return 0;
}

P4457 [BJOI2018] 治疗之雨

你有 \(p\) 滴血量，血量上限为 \(n\) 。每轮操作如下：

• 先以 \(\dfrac{1}{m + 1}\) 的概率增加 \(1\) 滴血，满血时则概率为 \(0\) 。
• \(k\) 次判定，每次以 \(\dfrac{1}{m + 1}\) 的概率减少一滴血，死了则概率为 \(0\) 。

求期望几轮死亡。

\(n \leq 1500\) ，\(m, k \leq 10^9\)

设 \(f_i\) 表示血量为 \(i\) 时期望多少轮结束，则：

\[f_i = 1 + \sum_{j = 0}^{\min(i, k)} (\dfrac{m}{m + 1} f_{i - j} + \dfrac{1}{m + 1} f_{i - j + 1}) \times \dbinom{k}{j} (\dfrac{1}{m + 1})^j (\dfrac{m}{m + 1})^{k - j} \]

注意一下 \(i = n\) 时的情况即可。

观察到 \(f_i\) 只与 \(f_{0 \sim i + 1}\) 有关，所以矩阵应该是一个类似下三角矩阵的东西。因为我们高斯消元的时候是拿自己这行去减下面的，所以每一行中只有 \(2\) 个系数要去和下面的相减，时间复杂度就可以做到 \(O(n^2)\) 了。

#include <bits/stdc++.h>
using namespace std;
const int Mod = 1e9 + 7;
const int N = 1.5e3 + 7;int g[N][N];
int inv[N], d[N], id[N], ans[N];int n, p, m, k;template <class T = int>
inline T read() {char c = getchar();bool sign = (c == '-');while (c < '0' || c > '9')c = getchar(), sign |= (c == '-');T x = 0;while ('0' <= c && c <= '9')x = (x << 1) + (x << 3) + (c & 15), c = getchar();return sign ? (~x + 1) : x;
}inline int add(int x, int y) {x += y;if (x >= Mod)x -= Mod;return x;
}inline int dec(int x, int y) {x -= y;if (x < 0)x += Mod;return x;
}inline int mi(int a, int b) {int res = 1;for (; b; b >>= 1, a = 1ll * a * a % Mod)if (b & 1)res = 1ll * res * a % Mod;return res;
}inline void prework() {inv[0] = inv[1] = 1;for (int i = 2; i < N; ++i) inv[i] = 1ll * (Mod - Mod / i) * inv[Mod % i] % Mod;
}inline void Gauss(int n, int band) {iota(id + 1, id + 1 + n, 1);for (int i = 1; i <= n; ++i) {if (!g[i][i]) {for (int j = i + 1; j <= min(n, i + band); ++j)if (g[i][j]) {for (int k = 1; k <= n; ++k)swap(g[k][i], g[k][j]);swap(id[i], id[j]);break;}}int Inv = mi(g[i][i], Mod - 2);for (int j = i + 1; j <= n; ++j) {int div = 1ll * g[j][i] * Inv % Mod;for (int k = i; k <= min(n, i + band); ++k)g[j][k] = dec(g[j][k], 1ll * g[i][k] * div % Mod);g[j][n + 1] = dec(g[j][n + 1], 1ll * g[i][n + 1] * div % Mod);}}for (int i = n; i; --i) {ans[id[i]] = g[i][n + 1];for (int j = i + 1; j <= min(n, i + band); ++j)ans[id[i]] = dec(ans[i], 1ll * ans[id[j]] * g[i][j] % Mod);ans[id[i]] = 1ll * ans[id[i]] * mi(g[i][i], Mod - 2) % Mod;}
}signed main() {prework();int T = read();while (T--) {n = read(), p = read(), m = read(), k = read();if (!k) {puts("-1");continue;} else if (!m) {if (k == 1)puts("-1");elseprintf("%d\n", (p + k - 2 - (p == n)) / (k - 1));continue;}d[0] = 1ll * mi(m, k) * mi(mi(m + 1, k), Mod - 2) % Mod;for (int i = 1, invm = mi(m, Mod - 2); i <= min(n, k); ++i)d[i] = 1ll * d[i - 1] * invm % Mod * (k - i + 1) % Mod * inv[i] % Mod;memset(g, 0, sizeof(g));for (int i = 1, invm1 = mi(m + 1, Mod - 2); i < n; ++i) {g[i][n + 1] = g[i][i] = 1;for (int j = 0; j <= min(i, k); ++j) {g[i][i - j] = dec(g[i][i - j], 1ll * dec(1, invm1) * d[j] % Mod);g[i][i - j + 1] = dec(g[i][i - j + 1], 1ll * invm1 * d[j] % Mod);}}g[n][n + 1] = g[n][n] = 1;for (int i = 0; i <= min(n, k); ++i)g[n][n - i] = dec(g[n][n - i], d[i]);Gauss(n, 1);printf("%d\n", ans[p]);}return 0;
}

P6899 [ICPC2014 WF] Pachinko

有一个宽度为 \(w\) 高度为 \(h\) 的方格纸， $ w \times h$ 的格子中，有一些是空的，有一些是洞，有一些是障碍物。从第一行的空的格子中随机选一个放置一个球，向上下左右移动的概率比为 \(p_u : p_d : p_l : p_r\) ，不能移动到有障碍物的格子上。对于每个洞，输出落入该洞的概率。

\(2 \leq 20, h \leq 10^4, p_u + p_d + p_l + p_r = 100\) 。

和上面比较类似。

#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-12;
const int dx[] = {-1, 1, 0, 0};
const int dy[] = {0, 0, -1, 1};
const int N = 1e4 + 7, M = 2e1 + 7;double g[N * M][M << 1], ans[N * M];
int id[N][M];
char str[N][M];int n, m, p[4], tot;inline bool check(int x, int y) {return x && x <= n && y && y <= m;
}inline void Gauss(int n, int band) {for (int i = 1; i <= n; ++i) {if (fabs(g[i][i]) < eps)continue;for (int j = i + 1; j <= min(i + band, n); ++j)g[i][j] /= g[i][i];ans[i] /= g[i][i], g[i][i] = 1;for (int j = i + 1; j <= min(i + band, n); ++j) {double div = g[j][i] / g[i][i];for (int k = i; k <= min(i + band, n); ++k)g[j][k] -= g[i][k] * div;ans[j] -= ans[i] * div;}}for (int i = n; i; --i)for (int j = i + 1; j <= min(i + band, n); ++j)ans[i] -= g[i][j] * ans[j];
}signed main() {scanf("%d%d", &m, &n);for (int i = 0; i < 4; ++i)scanf("%d", p + i);for (int i = 1; i <= n; ++i) {scanf("%s", str[i] + 1);for (int j = 1; j <= m; ++j)if (str[i][j] != 'X')id[i][j] = ++tot;}int sum = m - count(id[1] + 1, id[1] + 1 + m, 0);for (int i = 1; i <= n; ++i)for (int j = 1; j <= m; ++j) {if (!id[i][j])continue;if (i == 1)ans[id[i][j]] = 1.0 / sum;if (str[i][j] == 'T')continue;g[id[i][j]][id[i][j]] = 1;int base = 100;for (int k = 0; k < 4; ++k) {int x = i + dx[k], y = j + dy[k];if (!check(x, y) || !id[x][y])base -= p[k];}for (int k = 0; k < 4; ++k) {int x = i + dx[k], y = j + dy[k];if (check(x, y) && id[x][y])g[id[x][y]][id[i][j]] = -1.0 * p[k] / base;}}Gauss(tot, m);for (int i = 1; i <= n; ++i)for (int j = 1; j <= m; ++j)if (str[i][j] == 'T')printf("%.9lf\n", ans[id[i][j]]);return 0;
}