$CF24D\ Broken Robot\ DP+$高斯消元

Luogu

 

Description

你收到的礼物是一个非常聪明的机器人,行走在一块长方形的木板上.不幸的是,你知道它是坏的,表现得相当奇怪（随机）.该板由n行和m列的单元格组成.机器人最初是在i行和j列的某个单元格上.然后在每一步机器人可以到另一个单元.目的是去底层（n次）行.机器人可以停留在当前单元,向左移动,向右边移动,或者移动到当前下方的单元.如果机器人在最左边的列不能向左移动,如果它是在最右边的列不能向右移动.在每一步中,所有可能的动作都是同样可能的.返回步的预期数量达到下面的行.

 

Sol

这题和传纸条有点类似,不同的是机器人既能向左走又能向右走,而传纸条只能向右传.

f[i][j]表示从(i,j)走到最后一行所需要的期望步数

f[i][1]=1/3(f[i][1]+f[i][2]+f[i+1][1])+1

f[i][m]=1/3(f[i][m]+f[i][m-1]+f[i+1][m])+1

(j!=1&&j!=m)f[i][j]=1/4(f[i][j]+f[i][j-1]+f[i][j+1]+f[i+1][j])+1

部分状态之间可以互相转移互相影响,并不能满足DP的无后效性

所以不能线性递推,要用高斯消元直接求出状态转移方程的解

需要注意的是m=1的情况要特判

m=3时,系数矩阵如下:

-2/3  1/3    0

1/4    1/4    1/4

0        1/3    -2/3

值得一提的是，我们用f[i][j]表示从(i,j)走到最后一步的期望步数,按照行号倒序进行递推,而不是用f[i][j]表示从(x,y)到(i,j)的期望步数并正序递推.原因是,若正序递推,则还须求出(x,y)到最后一行每一个位置的概率p[n][j],计算Σp[n][j]*f[n][j]才能得到答案,较为复杂.

事实上,很多数学期望DP都会采取倒推的方式执行.

 

Code

关于代码实现其实还有几个值得注意的地方

1.f[i][j]数组可以滚动优化

2.因为要多次解方程组,很多人可能会每次都初始化a数组(系数矩阵).其实没有必要,只要初始化一次,并且同时记录一下c[i]=a[i+1][i]/a[i][i]就好了(没写code这句话可能暂时看不懂,看下code就会懂了鸭QwQ)

 

#include<iostream>
#include<cstdio>
#include<cstring>
#define Rg register
#define il inline
#define db double
#define ll long long
#define mem(a,b) memset(a,b,sizeof(a));
#define go(i,a,b) for(Rg int i=a;i<=b;++i)
#define yes(i,a,b) for(Rg int i=a;i>=b;--i)
using namespace std;
il int read()
{
    int x=0,y=1;char c=getchar();
    while(c<'0'||c>'9'){if(c=='-')y=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=(x<<3)+(x<<1)+c-'0';c=getchar();}
    return x*y;
}
const int N=1001;
int n,m,x,y;
db a[N][N],b[N],f[N],c[N];
il void init_a()
{
    if(m==1){a[1][1]=(db)-1.0/2;return;}
    a[1][1]=(db)-2.0/3;a[1][2]=(db)1.0/3;
    a[m][m-1]=(db)1.0/3;a[m][m]=(db)-2.0/3;
    go(i,2,m-1)a[i][i-1]=(db)1.0/4,a[i][i]=(db)-3.0/4,a[i][i+1]=(db)1/4;
    go(i,1,m)
    {
        c[i]=a[i+1][i]/a[i][i];
        a[i+1][i]-=c[i]*a[i][i];
        a[i+1][i+1]-=c[i]*a[i][i+1];
    }
}
il void init_b()
{
    if(m==1){b[1]=(db)-f[1]/2-1;return;}
    b[1]=-(db)f[1]/3-1;b[m]=-(db)f[m]/3-1;
    go(i,2,m-1)b[i]=-(db)f[i]/4-1;
}
il void calc()
{
    go(i,1,m){db t=c[i];b[i+1]-=t*b[i];}
    f[m]=b[m]/a[m][m];
    yes(i,m-1,1)f[i]=(b[i]-f[i+1]*a[i][i+1])/a[i][i];
}
int main()
{
    n=read(),m=read(),x=read(),y=read();//m==1 special case !
    init_a();
    yes(i,n-1,x){init_b();calc();}
    printf("%.10lf",f[y]);
    return 0;
}