『Broken Robot 后效性dp 高斯消元』

<更新提示>

---

<正文>

Broken Robot

Description

你作为礼物收到一个非常聪明的机器人走在矩形板上。不幸的是，你明白它已经破碎并且行为相当奇怪（随机）。该板由N行和M列单元组成。机器人最初位于第i行和第j列的某个单元格中。然后在每一步，机器人都可以去另一个细胞。目的是走到最底层（N.排。机器人可以停留在当前单元格中，向左移动，向右移动或移动到当前单元格下方的单元格。如果机器人位于最左侧的列中，则它不能向左移动，如果它位于最右侧的列中，则它不能向右移动。在每一步中，所有可能的动作都是同样可能的。返回预期的步数以到达最下面一行。

Input Format

在第一行中你将被提供两个空间隔开的整数N和M(1≤N,M≤?1000)。在第二行中你将得到另外两个空间隔开的整数i和j(1≤i≤N,1≤j≤M) ——初始行的数目和初始列的数量。注意，（1,1）是板的左上角，（N， M）是右下角。

Output Format

在自身的一行上输出预期的步数，小数点后至少有4位数。

Sample Input

10 14
5 14

Sample Output

18.0038068653

解析

显然，这道题看起来很像一道期望\(dp\)，那么我们就用期望\(dp\)的套路设一个状态试试：\(f[i][j]\)代表机器人从\((i,j)\)走到最后一行的步数期望值。

我们可以根据移动规则很容易列出\(dp\)方程：

\(1.\) 当机器人处于第一列时：\(f[i][1]=\frac{1}{3}(f[i][1]+f[i][2]+f[i+1][1])+1\)

\(2.\) 当机器人处于最后一列时：\(f[i][m]=\frac{1}{3}(f[i][m]+f[i][m-1]+f[i+1][m])+1\)

\(3.\) 当机器人处于中间列时：\(f[i][j]=\frac{1}{4}(f[i][j]+f[i][j-1]+f[i][j+1]+f[i+1][j])+1\)

我们发现，在行维度上，这个状态转移方程时没有问题的，可以倒序枚举每一行作为阶段，来进行转移。

但是，在同一行中，这个状态转移方程并不满足无后效性这一动态规划基本原则，于是我们决定使用高斯消元算法来解方程。

总体上，我们还是以行号为阶段，倒序进行转移。在第\(i\)行行内，我们将\(i+1\)行的状态看为常数，剩下的状态看做\(m\)个未知数，\(m\)个状态转移方程看做数学方程，尝试列出增广矩阵。

先看第一类方程：\(f[i][1]=\frac{1}{3}(f[i][1]+f[i][2]+f[i+1][1])+1\)，简单做一下移项分类：

\[2f[i][1]-f[i][2]=f[i+1][1]+3 \]

同理，可以化简剩下两个方程：

\[-f[i][j-1]+3f[i][j]-f[i][j+1]=f[i+1][j]+4 \]

\[-f[i][m-1]+2f[i][m]=f[i+1][m]+3 \]

那就可以写出系数矩阵了：

\[\left[ \begin{array}{ccccccc|c} 2 & -1 & 0 & \cdots & 0 & 0 & 0 & f[i+1][1]+3 \\ -1 & 3 & -1 & \cdots & 0 & 0 & 0 & f[i+1][2]+4 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots& -1 & 3 & -1 & f[i+1][m-1]+4 \\ 0 & 0 & 0 & \cdots& 0 & -1 & 2 & f[i+1][m]+3 \\ \end{array} \right] \]

我们发现这个矩阵很特殊，总共只有三条斜列有系数，所以我们可以线性直接消元。具体的说，我们可以从上往下消一遍，将第一斜列的系数消去，同时正确地处理第二第三斜列。再从下往上消一遍，将第三斜列的系数消去，这样就可以直接计算答案了。

利用如上的\(dp\)以及高斯消元算法，时间复杂度为\(O(nm)\)。

值得注意的是，这三个方程在\(m=1\)是会出现边界问题，简单推导可知\(m-1\)时答案就是\((n-x)*2\)。

\(Code:\)

#include <bits/stdc++.h>
using namespace std;
const int N = 1020;
int n,m,x,y;
double f[N][N],a[N][N],b[N];
inline void input(void)
{
    scanf("%d%d\n%d%d",&n,&m,&x,&y);
}
inline void init(void)
{
    a[1][1] = 2.0 , a[1][2] = -1.0;
    for (int i=2;i<m;i++)
        a[i][i-1] = -1.0 , a[i][i] = 3.0 , a[i][i+1] = -1.0;
    a[m][m-1] = -1.0 , a[m][m] = 2.0;
}
inline void gauss(void)
{
    double rate = a[2][1] / a[1][1];
    a[2][1] = 0.0;
    a[2][2] -= a[1][2] * rate , b[2] -= b[1] * rate;
    for (int i=2;i<m;i++)
    {
        rate = a[i+1][i] / a[i][i];
        a[i+1][i] = 0.0;
        a[i+1][i+1] -= a[i][i+1] * rate , b[i+1] -= b[i] * rate;
    }
    for (int i=m-1;i>=1;i--)
    {
        rate = a[i][i+1] / a[i+1][i+1];
        a[i][i+1] = 0.0 , b[i] -= b[i+1] * rate; 
    }
}
inline void dp(void)
{
    for (int i=1;i<=m;i++)
        f[n][i] = 0;
    for (int i=n-1;i>=1;i--)
    {
        init();
        b[1] = f[i+1][1] + 3.0;
        for (int j=2;j<m;j++) 
            b[j] = f[i+1][j] + 4.0;
        b[m] = f[i+1][m] + 3.0;
        gauss();
        for (int j=1;j<=m;j++) 
            f[i][j] = b[j] / a[j][j];
    }
}
int main(void)
{
    freopen("in.in","r",stdin);
    freopen("out.out","w",stdout);
    input();
    if ( m != 1 ) dp();
    else f[x][y] = ( n - x ) * 2.0;
    printf("%.10lf\n",f[x][y]);
    return 0;
}

---

<后记>