CodeForces 24 D.Broken robot（概率DP+高斯消元）

Description

给出一个$n\times m$的矩阵区域，一个机器人初始在第$x$行第$y$列，每一步机器人会等概率的选择停在原地，左移一步，右移一步，下移一步，如果机器人在边界则不会往区域外移动，问机器人到达最后一行的期望步数

Input

第一行两个整数$n,m$表示矩阵行列数，之后输入两个整数$x,y$表示机器人初始位置$(1\le n,m\le 1000,1\le x\le n,1\le y\le m)$

Output

输出机器人走到最后一行的期望步数

Sample Input

10 10
10 4

Sample Output

0.0000000000

Solution

$m=1$时每次$\frac{1}{2}$概率不动，$\frac{1}{2}$概率下移一步，答案显然为$2\cdot (n-x)$

$dp[i]][j]$表示机器人从第$i$行第$j$列出发到达最后一行的期望步数，显然$dp[n][j]=0,1\le j\le m$

由机器人等概率选择停在原地，左移一步，右移一步，下移一步有转移：

$dp[i][1]=1+\frac{1}{3}dp[i][1]+\frac{1}{3}dp[i][2]+\frac{1}{3}dp[i+1][1]$

$dp[i][j]=1+\frac{1}{4}dp[i][j]+\frac{1}{4}dp[i][j-1]+\frac{1}{4}dp[i][j+1]+\frac{1}{4}dp[i+1][j],1<j<m$

$dp[i][m]=1+\frac{1}{3}dp[i][m]+\frac{1}{3}dp[i][m-1]+\frac{1}{3}dp[i+1][m]$

写成矩阵形式

$dp[i][1]-\frac{1}{2}dp[i][2]=\frac{3}{2}+\frac{1}{2}dp[i+1][1]$

$-\frac{1}{4}dp[i][j-1]+\frac{3}{4}dp[i][j]-\frac{1}{4}dp[i][j+1]=1+\frac{1}{4}dp[i+1][j],1<j<m$

$-\frac{1}{2}dp[i][m-1]+dp[i][m]=\frac{3}{2}+\frac{1}{2}dp[i+1][m]$

由于转移矩阵是一个三对角矩阵，$O(n)$即可消元，即通过第一行的$dp[i][1]+a[1]\cdot dp[i][2]=b[1]$和第二行消去$dp[i][1]$得到$dp[i][2]+a[2]\cdot dp[i][3]=b[2]$，以此类推求出$dp[i][j]+a[j]\cdot dp[i][j+1]=b[j],1\le j<m$，用$dp[i][m-1]+a[m-1]\cdot dp[i][m]=b[m-1]$以及最后一个等式求出$dp[i][m]$，然后从后往前通过$dp[i][j]=b[j]-a[j]\cdot dp[i][j+1]$求出所有的$dp[i][j]$即可，时间复杂度$O(n^2)$

Code

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#include<queue>
#include<map>
#include<set>
#include<ctime>
using namespace std;
typedef long long ll;
typedef pair<int,int>P;
const int INF=0x3f3f3f3f,maxn=1005;
int n,m,x,y;
double a[maxn],b[maxn],dp[maxn][maxn];
int main()
{
    scanf("%d%d%d%d",&n,&m,&x,&y);
    if(m==1)printf("%.5f\n",2.0*(n-x));
    else
    {
        for(int j=1;j<=m;j++)dp[n][j]=0;
        for(int i=n-1;i>=x;i--)
        {
            a[1]=-0.5;b[1]=1.5+0.5*dp[i+1][1];
            for(int j=2;j<m;j++)
            {
                b[j]=1.0+0.25*dp[i+1][j]+0.25*b[j-1];
                a[j]=-0.25;
                double temp=0.75+0.25*a[j-1];
                a[j]/=temp;b[j]/=temp;
            }
            dp[i][m]=(3.0+dp[i+1][m]+b[m-1])/(2.0+a[m-1]);
            for(int j=m-1;j>=1;j--)dp[i][j]=b[j]-a[j]*dp[i][j+1];
        }
        printf("%.5f\n",dp[x][y]);
    }
    return 0;
}