codeforece 24D. Broken robot dp+高斯消元

时间:2022-11-16 07:54:37

You received as a gift a very clever robot walking on a rectangular board. Unfortunately, you understood that it is broken and behaves rather strangely (randomly). The board consists of N rows and M columns of cells. The robot is initially at some cell on the i-th row and the j-th column. Then at every step the robot could go to some another cell. The aim is to go to the bottommost (N-th) row. The robot can stay at it’s current cell, move to the left, move to the right, or move to the cell below the current. If the robot is in the leftmost column it cannot move to the left, and if it is in the rightmost column it cannot move to the right. At every step all possible moves are equally probable. Return the expected number of step to reach the bottommost row.

Input
On the first line you will be given two space separated integers N and M (1 ≤ N, M ≤ 1000). On the second line you will be given another two space separated integers i and j (1 ≤ i ≤ N, 1 ≤ j ≤ M) — the number of the initial row and the number of the initial column. Note that, (1, 1) is the upper left corner of the board and (N, M) is the bottom right corner.

Output
Output the expected number of steps on a line of itself with at least 4 digits after the decimal point.

Examples
input
10 10
10 4
output
0.0000000000
input
10 14
5 14
output
18.0038068653

大意:在(n,m)的矩阵里有一个点(x,y),每次随机向下,向右,向左,或者不走,不会走出边界,问走到最后一行的期望步数。

分析:
一道比较显然的dp……

我们设f[i][j]为从第i行第j列走到最后一行行动次数的期望值,注意是以(i,j)为起点,那么f[x][y]就是答案。

在第一列时,不能往左走,所以各有1/3的概率原地不动,往下走,往右走,也就是从这三个状态转移过来,

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

第m列同理,不能往右边走,

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

其他的位置就有4个方向,也就是4个状态,就是

f[i][j]=(1/4)*(f[i][j]+f[i][j-1]+f[i][j+1]+f[i+1][j])+1

显然,f[n][j]=0 (1<=j<=m)

然后我们发现,其实每行是不满足无后效性的。

于是就可以用高斯消元来搞。

为了好懂,我们设M=f[i][j],L=f[i][j-1],R=f[i][j+1],D=f[i+1][j]。显然D是已知的,因为这个是上一轮的状态,所以考虑转移,有

M=(1/4) M+(1/4) L+(1/4) R+(1/4) D+1

移项可得,

(1/4) L+(1/4) R+(-3/4) M=-(1/4) D-1

左右两边稍微有点不同,应该能想到。我们把f[i+1][j]看作已知数,于是对于
f[i][1]~f[i][m]这些未知量,可以列出
m个方程,可以用高斯消元求出f[i][1]~f[i][m]。

但是高斯消元是O(n^3)的,我们可以观察系数矩阵,我们发现,每行只有不多于3个系数不为0,而且这些系数都在对角线上,其实这个消元是O(n)的。
比如,当M=4时的系数矩阵。

-2/3 1/3 0 0

1/4 -3/4 1/4 0

0 1/4 -3/4 1/4

0 0 1/3 -2/3

大概就这样的,看看就好。我们消元求出上三角矩阵,每次只需要消掉下一行的系数,而且只需消掉2个系数。而求和时,由于一行只有两个系数不为0,可以直接求出解。

代码:

#include <cstdio>
#include <iostream>
#include <cmath>

const int maxn=1005;

using namespace std;

double f[maxn],a[maxn][maxn],b[maxn];

int x,y,n,m,i;

void build(int m)
{
    int i;
    if (m==1)
    {
        a[1][1]=(double)-1/2;
        b[1]=(double)-f[1]/2-1;
        return;
    }
    a[1][1]=(double)-2/3; a[1][2]=(double)1/3;
    a[m][m-1]=(double)1/3; a[m][m]=(double)-2/3;
    b[1]=(double)-f[1]/3-1; b[m]=(double)-f[m]/3-1;
    for (i=2;i<m;i++)
    {
        a[i][i-1]=a[i][i+1]=(double)1/4;
        a[i][i]=(double)-3/4;
        b[i]=(double)-f[i]/4-1;
    }
}

void guass()
{
    int i,j;    
    for (i=1;i<m;i++)
    {
        double rate=a[i+1][i]/a[i][i];
        a[i+1][i]-=rate*a[i][i];
        a[i+1][i+1]-=rate*a[i][i+1];
        b[i+1]-=rate*b[i];
    }   
    f[m]=b[m]/a[m][m];  
    for (i=m-1;i>0;i--)
    {
        f[i]=(b[i]-f[i+1]*a[i][i+1])/a[i][i];
    }
}
int main()
{
    scanf("%d%d",&n,&m);
    scanf("%d%d",&x,&y);    
    for (i=n-1;i>=x;i--)
    {
        build(m);
        guass();
    }
    printf("%.10lf",f[y]);
}