做了这道题才知道高斯消元可以这样用
$\texttt{Upd on 2021.8.18}$ 修复 $\LaTeX$ 的锅
$\texttt{Upd on 2022.7.21}$ 改正矩阵系数

题意

这题题意很明了，就是在一个 $N * M$ 的矩阵上走，一开始在位置 $(X, Y)$，每次等概率选择以下四种：

• 不动 $(X, Y)$
• 向左移动一格 $(X, Y - 1)$
• 向右移动一格 $(X, Y + 1)$
• 向下移动一格 $(X + 1, Y)$

不能移出矩阵，问从初始位置移动到最后一行的期望步数。

算法

这种在矩阵上移动求 $Min, Max, Sum$ 的几乎都是 $DP$，比如去年 CSP-J 的第 4 题。

既然确认了是 $DP$ 那就推状态转移方程吧，注意这里使用倒推，即设 $dp[i][j]$ 为从 $(i,j)$ 出发，走到第 $N$ 行的期望步数，而不是表示从 $(x,y)$ 到 $(i,j)$ 的期望步数，不然会很麻烦。

$ \begin{cases} dp[i][j]=\frac{1}{3}(dp[i][j]+dp[i][j+1]+dp[i+1][j])+1 \ \ \ if\ j=1\\\ dp[i][j]=\frac{1}{3}(dp[i][j]+dp[i][j-1]+dp[i+1][j])+1 \ \ \ if\ j=M\\\ dp[i][j]=\frac{1}{4}(dp[i][j]+dp[i][j-1]+dp[i][j+1]+dp[i+1][j])+1 \ \ \ else \end{cases} $

其中，$j=1$ 不能向左移动；$j=M$ 不能向右移动。

特殊情况 $M=1$，只能不动或向下，即 $dp[i][1]=\frac{1}{2}(dp[i][1]+dp[i+1][1])+1$

手动解一下方程，发现有 $dp[i][1]=dp[i+1][1]+2$

那么答案就可以直接算出，为 $2*(N-X)$

但是，发现状态转移方程有后效性，$dp[i][j]$ 会用到 $dp[i][j+1]$，而 $dp[i][j+1]$ 也会用到 $dp[i][j]$，构成了环。

但这题很良心，没有向上移动，行与行之间没有依赖关系，第 $i$ 行只会用到第 $i+1$ 行，不会用到第 $i-1$ 行，所以从 $N$ 开始倒序循环，就可以解决行。

于是，就可以用高斯消元解决列与列之间的依赖关系啦。

高斯消元是用来解决 $N$ 元一次方程组的，需要一个系数矩阵，具体算法过程分三步：

1. 找到第 $i$ 列最大的数并将其所在行挪到当前最靠上的一行，第 $i$ 行。
2. 将该数变为 $1$（将第 $i$ 行所有数除以该数）。
3. 将所有靠下行第 $i$ 列上的数变为 $0$（将所有行数比 $i$ 大的行上的数都减去该行第 $i$ 列数乘以第 $i$ 行上于其在同列的数）。

最后就转化成了一个上三角矩阵，可以轻松求出未知数，且这题不存在无解或有多组解。

$DP$ 状态转移方程在每一行上和高斯消元可以求的一样，是一个 $N$ 元一次方程组。所以可以倒序遍历每一行，每一行再填一下高斯消元系数矩阵就好了。

所以，现在的问题就是，高斯消元的系数矩阵怎么填，其实化简一下状态转移方程就很明了了。

$ \begin{cases} \frac{2}{3}dp[i][j]=\frac{1}{3}dp[i][j+1]+\frac{1}{3}dp[i+1][j]+1 \ \ \ if\ j=1\\\ \frac{2}{3}dp[i][j]=\frac{1}{3}dp[i][j-1]+\frac{1}{3}dp[i+1][j]+1 \ \ \ if\ j=M\\\ \frac{3}{4}dp[i][j]=\frac{1}{4}dp[i][j-1]+\frac{1}{4}dp[i][j+1]+\frac{1}{4}dp[i+1][j]+1 \ \ \ else \end{cases} $

其中 $dp[i+1][j]$ 是已经算出来的，那么把已知数和未知数分开就变成了：

$ \begin{cases} \frac{2}{3}dp[i][j]-\frac{1}{3}dp[i][j+1]=\frac{1}{3}dp[i+1][j]+1 \ \ \ if\ j=1\\\ \frac{2}{3}dp[i][j]-\frac{1}{3}dp[i][j-1]=\frac{1}{3}dp[i+1][j]+1 \ \ \ if\ j=M\\\ \frac{3}{4}dp[i][j]-\frac{1}{4}dp[i][j-1]-\frac{1}{4}dp[i][j+1]=\frac{1}{4}dp[i+1][j]+1 \ \ \ else \end{cases} $

但是高斯消元的时间复杂度是 $O(N^3)$ 级别，乘上行数 $N$，$O(N^4)$ 时间复杂度妥妥爆炸。

不慌，来看一下一个 $5*5$ 的系数矩阵：

$ \begin{bmatrix} \frac{2}{3} & -\frac{1}{3} & 0 & 0 & 0 & \ \vert\ & \frac{1}{3}dp[i+1][1]+1\\\ -\frac{1}{4} & \frac{3}{4} & -\frac{1}{4} & 0 & 0 & \vert & \frac{1}{4}dp[i+1][2]+1\\\ 0 & -\frac{1}{4} & \frac{3}{4} & -\frac{1}{4} & 0 & \vert & \frac{1}{4}dp[i+1][3]+1\\\ 0 & 0 & -\frac{1}{4} & \frac{3}{4} & -\frac{1}{4} & \vert & \frac{1}{4}dp[i+1][4]+1\\\ 0 & 0 & 0 & -\frac{1}{3} & \frac{2}{3} & \vert & \frac{1}{3}dp[i+1][5]+1\\\ \end{bmatrix} $

发现只有主对角线和两旁的对角线以及最后一列上的数非 $0$。

再来看一下三个步骤。

1. 找到第 $i$ 列最大的数并将其所在行挪到当前最靠上的一行，第 $i$ 行。
2. 将该数变为 $1$（将第 $i$ 行所有数除以该数）。
3. 将所有靠下行第 $i$ 列上的数变为 $0$（将所有行数比 $i$ 大的行上的数都减去该行第 $i$ 列数乘以第 $i$ 行上于其在同列的数）。

步骤 $1$ 中不用找，就是当前最靠上的一行，已经在最靠上的一行了，不用换。

步骤 $2$ 中该行被除后有变化，即非 $0$ 的，就只有 $(i,i),(i,i+1),(i,M+1)$ (注意 $(i,i-1)$ 虽然一开始非 $0$，但已经被第 $i-1$ 列做的 $3$ 号操作变为 $0$ 了)，所以每次只要动 $3$ 个数就行。

步骤 $3$ 中一个行数比 $i$ 大的行只有第 $i$ 列上的数非 $0$，减的才有意义，而满足这一条件的，靠下的行就只有第 $i+1$ 行，且第 $i$ 行上非 $0$ 的位置只有 $(i,i),(i,i+1),(i,M+1)$，第 $i+1$ 行上的数减后会改变，也只用动 $3$ 个数。

所以，综上所述，就得到了一个时间复杂度 $O(NM)$ 的算法，$10^3*10^3=10^6$，稳过。

代码

代码里有注释

#include <cmath>
#include <cstdio>

using namespace std;

const double eps = 1e-6;
int N;
int M;
int X;
int Y;
double dp[1009][1009];
double A[1009][1009];
double res;

void In () {
    scanf("%d%d%d%d", &N, &M, &X, &Y);
}

void Gauss () {  //高斯消元
    double t;

    for (int i = 1; i <= M; i++) {
        t = A[i][i];
        A[i][i] /= t; A[i][i + 1] /= t;  //第二步
        if (i < M) {   //这里特判是因为如果 i=M 那 A[i][M+1] 已经在 A[i][i + 1] /= t 时除过 t 了，就没必要再除了
            A[i][M + 1] /= t;
        }
        t = A[i + 1][i];
        A[i + 1][i] -= t * A[i][i]; A[i + 1][i + 1] -= t * A[i][i + 1]; A[i + 1][M + 1] -= t * A[i][M + 1];   //第三步
    }

    for (int i = M; i; i--) {   //推出答案，存在 A[i][M+1] 中
        A[i - 1][M + 1] -= A[i - 1][i] * A[i][M + 1];
        A[i - 1][i] -= A[i - 1][i] * A[i][i];
    }
}

void Solve () {
    if (M == 1) {   //特判 M=1
        res = 2 * (N - X);
        return;
    }

    for (int i = N - 1; i >= X; i--) {
        A[1][1] = A[M][M] = 2.0 / 3;
        A[1][2] = A[M][M - 1] = -1.0 / 3;
        A[1][M + 1] = dp[i + 1][1] / 3 + 1;
        A[M][M + 1] = dp[i + 1][M] / 3 + 1;

        for (int m = 2; m < M; m++) {
            A[m][m - 1] = A[m][m + 1] = -1.0 / 4;
            A[m][m] = 3.0 / 4;
            A[m][M + 1] = dp[i + 1][m] / 4 + 1;
        }

        Gauss();

        for (int m = 1; m <= M; m++) {
            dp[i][m] = A[m][M + 1];
        }
    }

    res = dp[X][Y];   //倒推
}

void Out () {
    printf("%.4lf\n", res);
}

int main () {

    In();
    Solve();
    Out();

    return 0;
}