ZigZagK的博客
[期望DP+高斯消元+复杂度分析]Codeforces963E【Circles of Waiting】题解

题目概述

从原点出发,每次往上下左右走都有一定的概率,问第一次走到离原点距离超过 $R$ 的点的期望步数。

解题报告

很显然可以期望DP,令距离超过 $R$ 但最接近原点的一圈的 $f_{x,y}=0$ ,可以写出转移方程:

$$ f_{x,y}=p_0f_{x-1,y}+p_1f_{x,y-1}+p_2f_{x+1,y}+p_3f_{x,y+1}+1 $$

上高斯消元解出 $f$ ,时间复杂度 $O(R^6)$ ,不能承受。


将距离不超过 $R$ 的点按照从上到下从左到右的顺序标号,可以发现矩形大部分都是 $0$ ,每行只有 $2R$ 的长度内至少有一个元素不是 $0$ ,所以每次消元只需要在 $O(R^2)$ 的矩形中做,而不需要全部矩形做,这样复杂度就是 $O(R^4)$ 了。

示例程序

意识模糊石乐志,调样例的时候把一句话删掉了……结果查了半天……

#include<cstdio>
using namespace std;
typedef long long LL;
const int maxr=50,maxn=7845,MOD=1e9+7,df[4][2]={{-1,0},{0,-1},{1,0},{0,1}};

int R,p[4],n,tem[(maxr<<1)+5][(maxr<<1)+5],M[maxn+5][maxn+5],ans[maxn+5];

inline int Pow(int w,int b) {int s=1;for (;b;b>>=1,w=(LL)w*w%MOD) if (b&1) s=(LL)s*w%MOD;return s;}
#define ID(x,y) (tem[(x)+maxr+1][(y)+maxr+1])
inline void AMOD(int &x,int tem) {if ((x+=tem)>=MOD) x-=MOD;}
inline int Solve(){
    for (int i=1;i<=n;i++)
        for (int j=i+1,d=1;j<=n&&d<=(R<<1);j++,d++){
            int now=MOD-(LL)M[j][i]*Pow(M[i][i],MOD-2)%MOD;
            for (int k=i,d=0;k<=n&&d<=(R<<1);k++,d++) AMOD(M[j][k],(LL)M[i][k]*now%MOD);
            AMOD(M[j][0],(LL)M[i][0]*now%MOD);
        }
    for (int i=(ans[0]=1,n);i;i--){
        ans[i]=(LL)M[i][0]*ans[0]%MOD;
        for (int j=i+1;j<=n;j++) AMOD(ans[i],MOD-(LL)M[i][j]*ans[j]%MOD);
        ans[i]=(LL)ans[i]*Pow(M[i][i],MOD-2)%MOD;
    }
    return ans[ID(0,0)];
}
int main(){
    freopen("program.in","r",stdin);freopen("program.out","w",stdout);
    scanf("%d%d%d%d%d",&R,&p[0],&p[1],&p[2],&p[3]);int P=p[0]+p[1]+p[2]+p[3];
    for (int i=0;i<4;i++) p[i]=(LL)p[i]*Pow(P,MOD-2)%MOD;
    for (int x=-R;x<=R;x++) for (int y=-R;y<=R;y++) if (x*x+y*y<=R*R) ID(x,y)=++n;
    for (int x=-R;x<=R;x++)
        for (int y=-R;y<=R;y++){
            if (!ID(x,y)) continue;M[ID(x,y)][ID(x,y)]=M[ID(x,y)][0]=MOD-1;int xx,yy;
            for (int f=0;f<4;f++) if (ID(xx=x+df[f][0],yy=y+df[f][1])) M[ID(x,y)][ID(xx,yy)]=p[f];
        }
    return printf("%d\n",Solve()),0;
}
版权声明:本博客所有文章除特别声明外,均采用 CC BY 4.0 CN协议 许可协议。转载请注明出处!