『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\)，简单做一下移项分类：
同理，可以化简剩下两个方程：
那就可以写出系数矩阵了：
我们发现这个矩阵很特殊，总共只有三条斜列有系数，所以我们可以线性直接消元。具体的说，我们可以从上往下消一遍，将第一斜列的系数消去，同时正确地处理第二第三斜列。再从下往上消一遍，将第三斜列的系数消去，这样就可以直接计算答案了。
利用如上的\(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;
}
<后记>
