【CF24D】Broken Robot （DP+高斯消元）

题目链接
题意：给定一个\(n\times m\)的矩阵，每次可以向→↓←移动一格，也可以原地不动，求从\((x,y)\)到最后一行的期望步数。
此题标签\(DP\)
看到上面这个肯定会想到

方法一： \(f[i][j]\)表示表示从\((x,y)\)走到\((i,j)\)的期望步数，正推
方法二： \(f[i][j]\)表示从\((i,j)\)走到最后一行的期望步数，倒推

事实上，方法二更优秀。
因为如果用方法一，我们要求的答案就是\(\frac{\sum f[\text{最后一行}][\text{每个位置}]}{m}\)
而如果我们用方法二，答案就是\(f[x][y]\)，省时又省力。
所以我们就用方法二。

然而上面两种方案都是有后效性的，因为\(f[i][j]\)取决与\(f[i-1][j](\text{或}f[i+1][j],f[i][j-1],f[i][j+1]\)甚至他本身，我们无法确定一个正确的递推顺序。
这时就要用\(DP\)套高斯消元了。
显然\(f[\text{最后一行}][\text{每个位置}]\)都为\(0\)，于是枚举行，从\(n-1\)枚举到\(x\)。
先把状态转移方程写出来。

\(f[i][j]=\begin{cases}\ \frac{1}{3}(f[i][j]+f[i][j+1]+f[i+1][j])\ (j=1)\\ \frac{1}{4}(f[i][j-1]+f[i][j]+f[i][j+1]+f[i+1][j])\ (j!=1,j!=m)\\ \frac{1}{3}(f[i][j]+f[i][j-1]+f[i+1][j])\ (j=m)\end{cases}\)

可以发现，方程中\(f[i+1][j]\)固定出现，而这是我们已知的（最后一行均为\(0\)，一行一行往上倒推），于是考虑移项。
以\(j=1\)为例，

\[f[i][j]=\frac{1}{3}(f[i][j]+f[i][j+1]+f[i+1][j]) \]

移项得

\[-\frac{1}{3}f[i+1][j]=\frac{1}{3}(-2f[i][j]+f[i][j+1]+f[i+1][j]) \]

于是我们便得到了\(m\)个一次方程，用高斯消元即可解出每一行所有的\(f\)值。

高斯消元理论复杂度\(O(n^3)\)，但是我们仔细观察本题，每行事实上只有两个数需要消元，所以我们可以在\(O(m)\)的时间复杂度里完成高斯消元，算法总时间复杂度\(O(nm^2)\)。

#include <cstdio>
#include <cstdlib>
const int MAXN = 1010;
double f[MAXN][MAXN];
double M[MAXN][MAXN];
int n, m, x, y;
void Gauss(){            //高斯消元
    for(int i = 1; i <= m; ++i){
       double tmp = 1.0 / M[i][i];        //系数化一
       M[i][i] *= tmp; M[i][m + 1] *= tmp;
       if(i == m) break;
       M[i][i + 1] *= tmp; 
       tmp = M[i + 1][i] / M[i][i];         //下一行消掉该元
       M[i + 1][i] -= tmp * M[i][i]; M[i + 1][i + 1] -= tmp * M[i][i + 1]; M[i + 1][m + 1] -= tmp * M[i][m + 1];
    }
    for(int i = m - 1; i; --i) M[i][m + 1] -= M[i + 1][m + 1] * M[i][i + 1];   //回代
}
int main(){
    scanf("%d%d%d%d", &n, &m, &x, &y);
    for(int i = n - 1; i >= x; --i){
       M[1][1] = -1.0 + 1.0 / 3;         //
       M[1][2] = 1.0 / 3;
       for(int j = 2; j < m; ++j){
          M[j][m + 1] = (-f[i + 1][j]) / 4.0 - 1;
          M[j][j] = -1.0 + 1.0 / 4;
          M[j][j - 1] = M[j][j + 1] = 1.0 / 4;
       }
       M[m][m] = -1.0 + 1.0 / 3;
       M[m][m - 1] = 1.0 / 3;
       if(m == 1) M[1][1] = -1.0 + 1.0 / 2;
       M[1][m + 1] = (-f[i + 1][1]) / 3.0 - 1;
       M[m][m + 1] = (-f[i + 1][m]) / 3.0 - 1;   //构建矩阵
       if(m == 1) M[m][m + 1] = (-f[i + 1][m]) / 2.0 - 1;   //特判$m=1$的情况
       Gauss();
       for(int j = 1; j <= m; ++j)
          f[i][j] = M[j][m + 1];
    }
    printf("%.10lf", f[x][y]);
    return 0;
}