一道有后效性的DP,所谓有后效性,也就是后面状态会干扰前面的状态,这道题可以向左走,也可以向右走,显然有后效性。
目标状态是确定的,走到最后一行任意一个位置,而初始状态会随着$x,y$的改变而改变,不妨倒序处理此问题。
第$i$行的状态可以由第$i+1$行推出,第一列只能往左或上走,第$M$列只能往右或上走,需要特殊处理。另外,只有$1$行时,也需要处理。
方程为:
$f_{i,1}=\frac{1}{3}(f_{i,1}+f_{i,2}+f_{i+1,1})+1\\f_{i,m}=\frac{1}{3}(f_{i,m-1}+f_{i,m}+f_{i+1,m})+1\\f_{i,j}=\frac{1}{4}(f_{i,j}+f_{i,j-1}+f_{i,j+1}+f_{i+1,j})+1$

对于有后效性的DP,通常采取高斯消元解方程

然而这道题的高斯消元不用$O(m^3)$,只有$O(m)$。

设$x_1=f_{i,1},x_2=f_{i,2},\cdots,x_m=f_{i,m}$,已经解出了$i+1$层的答案,那么有方程组

$\begin{cases}\frac{1}{3}x_1+\frac{1}{3}x_2+0*x_3+0*x_4+\cdots+0*x_m+\frac{1}{3}f_{i+1,1}+1=x_1\\\frac{1}{4}x_1+\frac{1}{4}x_2+\frac{1}{4}x_3+0*x_4+\cdots+0*x_m+\frac{1}{4}f_{i+1,2}+1=x_2 \\ \cdots\cdots\\0*x_1+\cdots +0*x_{m-2}+\frac{1}{3}x_{m-1}+\frac{1}{3}x_m+\frac{1}{3}f_{i+1,m}+1=x_m \end{cases}$

常量丢右边,变量丢左边之后,就有

$\begin{cases}-\frac{2}{3}x_1+\frac{1}{3}x_2=-\frac{1}{3}f_{i+1,1}-1\\\frac{1}{4}x_1-\frac{3}{4}x_2+\frac{1}{4}x_3=-\frac{1}{4}f_{i+1,2}-1\\\cdots\cdots\end{cases}$

简化为系数矩阵之后,就是一个这样的矩阵:

$-\frac{2}{3}~~\frac{1}{3}~~0~~0~~\cdots~~0\\~~\frac{1}{4}-\frac{3}{4}\frac{1}{4}~~0~~\cdots ~0\\~~0~~~~\frac{1}{4}-\frac{3}{4}\frac{1}{4}\cdots~ 0\\~\cdots\cdots\cdots\cdots\cdots\cdots\\~~0~~~0~~0~~\cdots\frac{1}{3}~~-\frac{2}{3}$

具体就是$m*m$的一个系数矩阵。

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<cstdlib>
using namespace std;
const int N=1e3+5;
double d[N][N],f[N][N];
int n,m,x,y;
void gauss()
{
    for(int i=1;i<=m;i++)
    {
        double w=1.0/d[i][i];
        d[i][i]*=w;//使d[i][i]这个系数为1 
        d[i][m+1]*=w;//扩大倍数 
        if(i==m)break;
        d[i][i+1]*=w;
        /*消掉d[i+1][i]*/
        w=d[i+1][i]/d[i][i];
        d[i+1][i+1]-=w*d[i][i+1];
        d[i+1][i]-=w*d[i][i];
        d[i+1][m+1]-=w*d[i][m+1];
    }
    for(int i=m-1;i;i--)
        d[i][m+1]-=d[i+1][m+1]*d[i][i+1];//d[i][i+1]/d[i+1][i+1](d[i+1][i+1]=1)*d[i+1][m];
    //消去d[i][i+1],由于d[i][i]=1,因此解就出来了。 
}
int main()
{
    scanf("%d%d%d%d",&n,&m,&x,&y);
    for(int i=n-1;i>=x;i--)
    {
        d[1][1]=d[m][m]=-2/3.0;
        d[1][2]=d[m][m-1]=1/3.0;
        for(int j=2;j<m;j++)
        {
            d[j][m+1]=-f[i+1][j]/4.0-1;
            d[j][j]=-3/4.0;
            d[j][j-1]=d[j][j+1]=1/4.0;
        }
        if(m==1)d[1][1]=-1/2.0;
        d[1][m+1]=-f[i+1][1]/3.0-1;
        d[m][m+1]=-f[i+1][m]/3.0-1;
        if(m==1)d[m][m+1]=-f[i+1][m]/2.0-1;
        gauss();
        for(int j=1;j<=m;j++)f[i][j]=d[j][m+1];
    }
    printf("%.4lf\n",f[x][y]);
    return 0;
}