CF24D Broken robot 题解
lyxeason
·
·
题解
做了这道题才知道高斯消元可以这样用
附 yxc 大佬讲题视频
题意
这题题意很明了，就是在一个 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 元一次方程组的，需要一个系数矩阵，具体算法过程分三步：
- 找到第 i 列最大的数并将其所在行挪到当前最靠上的一行，第 i 行。
- 将该数变为 1（将第 i 行所有数除以该数）。
- 将所有靠下行第 i 列上的数变为 0（将所有行数比 i 大的行上的数都减去该行第 i 列数乘以第 i 行上于其在同列的数）。
最后就转化成了一个上三角矩阵，可以轻松求出未知数，且这题不存在无解或有多组解。
所以，现在的问题就是，高斯消元的系数矩阵怎么填，其实化简一下状态转移方程就很明了了。
\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{2}{3} & -\frac{1}{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)$ 的算法，$1e3*1e3=1e6$，稳过。
------------
### 代码
###### 代码里有注释
```cpp
#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;
}
```