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

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

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

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

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

(f[i][j]=egin{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;
}