TopCoder SRM 614 Div2 Hard TorusSailingEasy

Div1 Hardの制限を緩めた問題。
http://community.topcoder.com/stat?c=problem_statement&pm=13085

問題

トーラス状がグリッド上に分割されており、NxMの格子点がある。
左右および上下は循環しており、X座標(N-1)の右はX座標0、Y座標(M-1)の上はY座標0となっている。

プレイヤーが(x,y)にいると、各ターンにおいて(x+1,y+1)または(x-1,y-1)のどちらかの座標に確率半々で遷移する。
原点にいるプレイヤーが(goalx,goaly)に到達するまでのターン数の期待値を求めよ。

解法

N,M<=10なので格子点は高々100個である。
そのため、まずゴールにたどり着けるかどうかは(x+1,y+1)及び(x-1,y-1)を100回以上繰り返してみればわかる。

ここで、原点から(x+1,y+1)をA回か(x-1,y-1)をB回連続で行うとゴールにたどり着くとする。
するとここからはトーラスのことは忘れて1次元で考えることができる。

すなわち、1次元座標上で原点から+1または-1ずつ移動してAか-Bにたどりつけばいいことになる。
座標zからAまたは-Bにたどりつくまでの期待値をP(z)とすると、P(A)=P(-B)=0であり、P(0)を求めればよい。
ここで、-B < i < AにおいてはP(i)=1+(P(i-1)+P(i+1))/2の関係が成り立つ。
よってこれらの式から(A+B+1)元連立方程式ができるので、これを解けばよい。
A+B+1は高々101。ガウスの掃出し法でO(N^3)で解けば時間は余裕。

追記：
P(x)=(x+A)*(B-x)が成り立つのでP(0)=A*Bが答えとなるようです。
P(x)=1+(P(x-1)+P(x+1))/2も確かに成り立ちますね。
教えて頂いたおくら氏に感謝いたします。

const int MAT=301;
int Gauss(int n,double mat_[MAT][MAT],double v_[MAT],double r[MAT]) {
	int i,j,k;
	double mat[MAT][MAT],v[MAT];
	memmove(mat,mat_,sizeof(mat));
	memmove(v,v_,sizeof(v));
	
	FOR(i,n) {
		if(mat[i][i]==0) {
			for(j=i+1;j<n;j++) if(mat[j][i]) break;
			if(j>=n) return -1;
			FOR(k,n) swap(mat[i][k],mat[j][k]);
			swap(v[i],v[k]);
		}
		v[i]/=mat[i][i];
		for(k=n-1;k>=i;k--) mat[i][k]/=mat[i][i];
		for(j=i+1;j<n;j++) {
			v[j]-=v[i]*mat[j][i];
			for(k=n-1;k>=i;k--) mat[j][k]-=mat[i][k]*mat[j][i];
		}
	}
	
	for(i=n-1;i>=0;i--) {
		for(j=n-1;j>i;j--) v[i]-=mat[i][j]*v[j],mat[i][j]=0;
		r[i]=v[i];
	}
	return 0;
}


class TorusSailingEasy {
	public:
	double expectedTime(int N, int M, int goalX, int goalY) {
		int ps[2]={1000000,1000000};
		int x[2]={0,0},y[2]={0,0},i,j;
		
		double mm[MAT][MAT];
		double v[MAT],r[MAT];
		
		FOR(j,10000) {
			x[0]+=1;
			y[0]+=1;
			x[1]+=N-1;
			y[1]+=M-1;
			x[0]%=N;
			y[0]%=M;
			x[1]%=N;
			y[1]%=M;
			if(x[0]==goalX && y[0]==goalY) ps[0]=min(ps[0],j+1);
			if(x[1]==goalX && y[1]==goalY) ps[1]=min(ps[1],j+1);
		}
		if(ps[0]>=100000) return -1;
		ZERO(mm);
		ZERO(v);
		for(i=1;i<ps[0]+ps[1];i++) mm[i][i-1]=mm[i][i+1]=-0.5,v[i]=1;
		FOR(i,ps[0]+ps[1]+1) mm[i][i]=1;
		Gauss(ps[0]+ps[1]+1,mm,v,r);
		
		return r[ps[0]];
	}
};

まとめ

今回は101元方程式なので掃出し法で解いたけど、式がかなり疎なのでもっと効率的に解く方法もありそうだ。