TopCoder SRM472 Div2 Hard RectangularIsland

laycurseさんの放送みて自分もやってみた。
対数使って計算してやらないとどうも誤差死するみたい。

問題概要

幅width高さheightの島のx,yの位置に人が立っていて、ランダムに上下左右のいずれかの方向にstep回移動する。
このとき、島から人が落ちない確率を求めよ。

width,height≤1000
0≤ x ≤ width-1
0≤ y ≤ height-1
1≤ steps ≤ 5000である。

解法

steps回のうちi回上下移動したとする。i回の上下移動で島からはみ出る確率は一次元のランダムウォークなので動的計画法によりO(i*length)で求められる。
左右にも同様の計算を行い、(x方向ではみ出る確率)+(y方向ではみ出る確率)-(x方向ではみ出る確率)*(y方向ではみ出る確率)が、stepsのうちi回が上下移動になる移動ではみでる確率である。
あとはiを0からstepsまで動かして足し合わせればよい。

ただしi回の移動の選び方の計算を対数を使ってやらないと誤差死してしまう。
nCk/2^nを求めるのに毎回logを呼び出すと時間がかかりすぎるので、使う可能性のある5000ほどまで事前に計算しておく。

ソースコード

double l[5005];

double C(int n,int k){
  double r=0;
  rep(i,k)r+=l[n-i],r-=l[i+1];
  r-=l[2]*n;

  return exp(r);
}
void solve(int len,int sx,int s,double *p){
  double dp[2][len+2]; rep(i,2)fill_n(dp[i],len+2,0);
  dp[0][sx+1]=1;
  
  int cur=0,next=1;
  rep(i,s+1){
    p[i]=dp[cur][0]+dp[cur][len+1];

    fill_n(dp[next],len+2,0);
    dp[next][0]=dp[cur][0];
    dp[next][len+1]=dp[cur][len+1];
    for(int j=1;j<=len;j++){
      dp[next][j+1]+=dp[cur][j]/2;
      dp[next][j-1]+=dp[cur][j]/2;
    }
    swap(cur,next);
  }
}


class RectangularIsland {
	public:
	double theProbablity(int width, int height, int x, int y, int steps) {
	  
	  for(int i=1;i<=5004;i++)l[i]=log(i);
	  double p=1,p1[5005],p2[5005];
	  solve(width,x,steps,p1);
	  solve(height,y,steps,p2);

	  rep(i,steps+1){
	    double px=p1[i],py=p2[steps-i];
	    p-=(px+py-px*py)*C(steps,i);
	  }
	  return p;
	}
};