Hatena::Grouptopcoder

SRM diary(Sigmar)

SigmarのTopcoder SRM参加記録など雑記です。
社会人になってから競技プログラミングを始めました。どこまで行けるか分かりませんが合間を見つけてアルゴリズムの勉強をしています。

2012-05-20SRM543 Div1

SRM543 Div1 250 EllysXors

| 19:29 | SRM543 Div1 250 EllysXors - SRM diary(Sigmar) を含むブックマーク はてなブックマーク - SRM543 Div1 250 EllysXors - SRM diary(Sigmar) SRM543 Div1 250 EllysXors - SRM diary(Sigmar) のブックマークコメント

Problem Statement

コーディングフェーズ

異様に速い提出者がいて若干ライブラリゲーの様相を呈する

ライブラリは持っていなかったので普通に考えた


0~NまでのXORを単に求める問題に置換する

ビットごとに考える

  • 1ビット目は 0,1,0,1,0,1,0,1,...
  • 2ビット目は 0,0,1,1,0,0,1,1,...
  • 3ビット目は 0,0,0,0,1,1,1,1,...

となるから、xビット目は2^xの周期を持つ

1ビット目は、4で割った余りで場合分けが可能

2ビット目以降は、1周期の1のビット数が2の倍数だから、周期で割った余りだけ調べればいい


書いた

合わない


こまごましたoff by one的な間違いがあった

直した

提出


遅っ・・・!


後で

unsigned intにして、ループで計算して通した人がいたと聞いて驚いた

XORって2秒間に40億回も計算できるのね・・・

計算量の見積もりに精通していれば一瞬で書けるというのは、悪くない問題設定だと思う


ソースコード

提出したバージョン

class EllysXors {
public:
	ll calc(ll x) {
		ll res=0;
		if(x%4==1 || x%4==2) res|=1;
		for(int i=1; i<40 && x; i++) {
			ll rem=x%(1LL<<(i+1));
			if(rem>=(1LL<<i)) {
				rem-=(1LL<<i);
				if(rem%2==0) res|=(1LL<<i);
			}
		}
		return res;
	}

	long long getXor(long long L, long long R) {
		long long res;

		res=calc(R)^calc(L-1);
		
		return res;
	}
};

SRM543 Div1 500 EllysRivers

| 19:30 | SRM543 Div1 500 EllysRivers - SRM diary(Sigmar) を含むブックマーク はてなブックマーク - SRM543 Div1 500 EllysRivers - SRM diary(Sigmar) SRM543 Div1 500 EllysRivers - SRM diary(Sigmar) のブックマークコメント

Problem Statement

コーディングフェーズ

明らかに普通のDPは無理だよね


これが船が島のどこからでも(整数の位置でなくても)発着できて、歩くのが一切ないのであれば、単なる光の入射角を決める問題である

高校物理によると光の進む速度によって屈折率が決まり、光が通る道は常にその場所への最短路となるらしい

ということは、光の屈折率を使えば最短路が出せるということである

  • sin(θ1)/sin(θ2)=v1/v2 みたいな感じ

入射角が決まれば、屈折率に従って右端の島のどの位置に着くのかは一意に定まる

そこで、入射角を二分探索することで、目的地に到着するための入射角を決めることができる


walkはどこでやっても同じだから、左端の島だけで歩けばいいんだろう

ということは、どうせwalkの距離に対して到着時間が谷型の関数になるんだろうから、三分探索で出せるんだろう


あとは整数の位置しか港がないから、doubleで計算した港の位置の近くにある適当な港を使って、どれを使えば一番速いか適当に探索すればよさげ


書いてみようとする


・・・


タイムアップ


こんな複雑なの時間内に書けるわけなかった/(^o^)\


チャレンジフェーズ

適当に怪しそうな人にでかいケースをぶつけて+50

たまたま誤差死してくれたみたいで助かった

こんなお茶濁し的なやり方では長くは持たないのでちゃんと問題を解けなければ・・・


後で

Greedyで良かったっていう・・・


どこかで北にlength回進まなければならないので、最初は全ての島で下端にいるとして、1つ北に進むのはどの川 or 歩きを使えば一番速いかっていうのをGreedyにlength回求めるだけだった


DPで計算量が多すぎるときに、部分的にGreedyにしたり、そのために見方を変えるのは常套手段で、今回もそういう意味では典型的なんだけど、全く対応できてない、というか、ちゃんと考えられてない気がする


屈折率による解法でも解けたけど、ちょっと書くのが遅すぎるし、色々いけてなかった。。

まあ高校物理なんて覚えてないんで思い出すのに時間がかかったってのもありますが・・・


ソースコード

せっかくだから屈折率のコードをpracticeに通したので掲載します

港を選ぶ辺りが若干嘘解法っぽいが、lengthの大きさに耐性のある解法だからあながち意味のない考え方でもないかな?

class EllysRivers {
public:
	// 島の長さlengthとしたときの最短となる入射角を求める
	double bestangle(double length, const vector <int> &width, const vector <int> &speed) {
		int n=width.size();
		double hi=M_PI/2, lo=0, mid;
		for(int i=0; i<100; i++) {
			mid=(hi+lo)/2;
			bool ok=true;
			double h=(double)width[0]*tan(mid);
			double theta=mid;
			for(int j=0; j<n-1; j++) {
				double sinn=sin(theta)*(double)speed[j+1]/(double)speed[j];
				if(abs(sinn)>=1) { ok=false; break; } // 全反射
				theta=asin(sinn);
				h+=(double)width[j+1]*tan(theta);
			}
			if(!ok || h>length) hi=mid;
			else lo=mid;
		}
		return lo;
	}

	double getMin(int length, int walk, vector <int> width, vector <int> speed) {
		double res=1e100;
		int n=width.size();
		
		// 歩く距離の最適値を求める(double)
		double hi=length, lo=0, mid[2];
		for(int i=0; i<100; i++) {
			mid[1]=(hi-lo)/3;
			mid[0]=lo+mid[1];
			mid[1]+=mid[0];
			double theta[2];
			theta[0]=bestangle((double)length-mid[0], width, speed);
			theta[1]=bestangle((double)length-mid[1], width, speed);
			double tres[2];
			for(int j=0; j<2; j++) {
				tres[j]=mid[j]/(double)walk;
				tres[j]+=((double)width[0]/cos(theta[j]))/(double)speed[0];
				for(int k=0; k<n-1; k++) {
					theta[j]=asin(sin(theta[j])*speed[k+1]/speed[k]);
					tres[j]+=((double)width[k+1]/cos(theta[j]))/(double)speed[k+1];
				}
			}
			if(tres[0]<tres[1]) hi=mid[1];
			else lo=mid[0];
		}

		// 最適な入射角における各島の到着位置(南北方向)を求める(double)
		double theta=bestangle((double)length-lo, width, speed);
		vector <double> hv;
		double h=lo+(double)width[0]*tan(theta);
		hv.push_back(h);
		for(int j=0; j<n-1; j++) {
			double sinn=sin(theta)*(double)speed[j+1]/(double)speed[j];
			assert(abs(sinn)<1);
			theta=asin(sinn);
			h+=(double)width[j+1]*tan(theta);
			hv.push_back(h);
		}
		
		// 各島の到着位置の前後適当な範囲の港を使って最短路を求める
		const int range=10;
		int curp[range];
		for(int i=0; i<range; i++) curp[i]=(int)lo-(range-1)/2+i;
		double cur[range];
		for(int i=0; i<range; i++) cur[i]=abs((double)curp[i]/walk);
		for(int k=0; k<n; k++) {
			double nxt[range];
			for(int i=0; i<range; i++) nxt[i]=1e100;
			int nxtp[range];
			for(int i=0; i<range; i++) nxtp[i]=(int)hv[k]-(range-1)/2+i;
			for(int i=0; i<range; i++) for(int j=0; j<range; j++) {
				nxt[j]=min(nxt[j], cur[i]+hypot((double)width[k], (double)(nxtp[j]-curp[i]))/(double)speed[k]);
			}
			memcpy(cur, nxt, sizeof(cur));
			memcpy(curp, nxtp, sizeof(curp));
		}
		for(int i=0; i<range; i++) if(curp[i]==length) res=cur[i];
		
		return res;
	}
};
トラックバック - http://topcoder.g.hatena.ne.jp/jackpersel/20120520

2012-04-25SRM541 Div1

SRM541 Div1 250 AntsMeet

| 01:12 | SRM541 Div1 250 AntsMeet - SRM diary(Sigmar) を含むブックマーク はてなブックマーク - SRM541 Div1 250 AntsMeet - SRM diary(Sigmar) SRM541 Div1 250 AntsMeet - SRM diary(Sigmar) のブックマークコメント

Problem Statement

コーディングフェーズ

任意の蟻のペアで衝突判定をして、時間の早い順に並べて、死んだ蟻を取り除きながら真の衝突判定をする

とか思ったけどあまりにも面倒そうだったので他の手段を検討

単純なシミュレーションで間に合うことに気づいたので書く

サンプル合った

提出


チャレンジフェーズ

ボロボロと周りが落とされていくんだが全然何が悪いのか分からない

どうやら座標を2倍していない人がたくさんいたらしい

サンプルにそういう例があるのでまさか2倍せずに(or 0.5刻みにせずに)提出してる人がいるとは思っていなかった・・・


ソースコード

class AntsMeet {
public:
	int countAnts(vector <int> x, vector <int> y, string dir) {
		int res;
		int n=x.size();

		for(int i=0; i<n; i++) {
			x[i]*=2;
			y[i]*=2;
		}

		int dis[52]={0};
		for(int t=0; t<10000; t++) {
			for(int i=0; i<n; i++) {
				if(dir[i]=='N') y[i]++;
				if(dir[i]=='E') x[i]++;
				if(dir[i]=='S') y[i]--;
				if(dir[i]=='W') x[i]--;
			}
			for(int i=0; i<n; i++) {
				if(dis[i]) continue;
				for(int j=i+1; j<n; j++) {
					if(dis[j]) continue;
					if(x[i]==x[j] && y[i]==y[j]) {
						dis[i]=dis[j]=1;
					}
				}
			}
		}

		res=0;
		for(int i=0; i<n; i++) if(!dis[i]) res++;
		
		return res;
	}
};

SRM541 Div1 550 AkariDaisukiDiv1

| 01:12 | SRM541 Div1 550 AkariDaisukiDiv1 - SRM diary(Sigmar) を含むブックマーク はてなブックマーク - SRM541 Div1 550 AkariDaisukiDiv1 - SRM diary(Sigmar) SRM541 Div1 550 AkariDaisukiDiv1 - SRM diary(Sigmar) のブックマークコメント

Problem Statement

コーディングフェーズ

上限が1000万だから単純なループにする問題ですかね

だとすると、Xに対してWaaiとかAkariとかくっつけた時に解が増える増分だけ調べればいいのでは

これって、Xの先頭はWaaiWaaiWaai...になるし、後ろは...DaisukiDaisukiDaisukiになるから、50回以上繰り返せば増分は一定の値に収束するっぽい

Xが短い間は全探索して、ちょっと長くなったらprefixとsuffixだけ覚えて、50回以上になったらあとはループするだけ

やるだけ・・?

なぜ550


書く

増分の計算かなりめんどい

思ったよりは大変だった


できた

提出


ソースコード

const int MOD=1000000007;

class AkariDaisukiDiv1 {
public:
	int countF(string w, string a, string d, string S, string F, int k) {
		int res;

		string x=S;
		string px, sx;
		ll cnt=0;
		int inc=0;
		for(int i=0; i<min(51, k); i++) {
			if(px.empty()) {
				x=w+x+a+x+d;
				if(x.size()>=F.size()) {
					cnt=0;
					for(int j=0; j+F.size()<=x.size(); j++) {
						if(x.substr(j, F.size())==F) cnt++;
					}
				}
				if(x.size()>50) {
					px=x.substr(0, 50);
					sx=x.substr(x.size()-50);
					x.clear();
				}
			} else {
				inc=0;
				string tx=sx+a+px;
				for(int j=sx.size()-F.size()+1; j<=sx.size()+a.size()-1; j++) {
					if(tx.substr(j, F.size())==F) inc++;
				}
				px=w+px;
				sx=sx+d;
				for(int j=0; j<w.size(); j++) {
					if(px.substr(j, F.size())==F) inc++;
				}
				for(int j=sx.size()-d.size()-F.size()+1; j+F.size()<=sx.size(); j++) {
					if(sx.substr(j, F.size())==F) inc++;
				}
				px=px.substr(0, 50);
				sx=sx.substr(sx.size()-50);
				cnt=(cnt*2+inc)%MOD;
			}
		}
		
		for(int i=51; i<k; i++) {
			cnt=(cnt*2+inc)%MOD;
		}
		res=cnt%MOD;
		return res;
	}
};

PraveenPraveen2012/11/14 15:22You've really captured all the esstenials in this subject area, haven't you?

lltwektgdjilltwektgdji2012/11/17 00:38s3e55e <a href="http://weadgecnprbz.com/">weadgecnprbz</a>

zkwvhwizkwvhwi2012/11/17 10:47mzuI7w , [url=http://amjymbwsqtpw.com/]amjymbwsqtpw[/url], [link=http://jtoocgzxmaus.com/]jtoocgzxmaus[/link], http://uoyfejtwsvcz.com/

トラックバック - http://topcoder.g.hatena.ne.jp/jackpersel/20120425

2012-04-20SRM540 Div1

SRM540 Div1 250 ImportantSequence

| 16:55 | SRM540 Div1 250 ImportantSequence - SRM diary(Sigmar) を含むブックマーク はてなブックマーク - SRM540 Div1 250 ImportantSequence - SRM diary(Sigmar) SRM540 Div1 250 ImportantSequence - SRM diary(Sigmar) のブックマークコメント

Problem Statement

コーディングフェーズ

1つ決めると後は全部決まるのはすぐ分かったけど、それでどうすればいんだっけ・・

+がなければ明らかに上限がなくて解が無限にあるから、+のところで区間を区切って考えるのかなあ

何となく+のところで区切った区間ごとに上限と下限を出してみようとしてみる

区間同士はどうやって結合すればいいんだ・・・

あー、実は+のところは上限と下限を反転するだけなのか

じゃあ、区間ごとに区切らなくても、左端の要素から右へ順番になぞって、上限と下限を段々狭めていけばいいのか

書いた

時間かかりすぎた。提出。

あー部分的に0を返す判定が抜けてる・・まあ最後にも判定してるから大丈夫かな。。


チャレンジフェーズで勘違いしてオーバーフローしそうだと思ったけどしなかったんで-50点くらってしまった

ちゃんと見てから投げよう・・・


ソースコード

class ImportantSequence {
public:
	int getCount(vector <int> B, string op) {
		int res;
		int n=B.size();
		int finite=false;
		for(int i=0; i<n; i++) {
			if(op[i]=='+') finite=true;
		}
		if(!finite) return -1;

		ll maxv=1000000000000000000LL, minv=1;
		for(int i=0; i<n; i++) {
			if(op[i]=='-') {
				maxv-=B[i];
				minv-=B[i];
				if(minv<=0) minv=1;
				//ここに if(maxv<=0) return 0; もあったほうが良かったかな
			} else {
				if(minv>=B[i]) return 0;
				if(maxv>=B[i]) maxv=B[i]-1;
				swap(maxv, minv);
				maxv=B[i]-maxv;
				minv=B[i]-minv;
			}
		}
		
		res=maxv-minv+1;
		if(res<0) res=0;
		return res;
	}
};

KelisKelis2012/11/15 02:30Me and this atircle, sitting in a tree, L-E-A-R-N-I-N-G!

vsolmndcevsolmndce2012/11/15 12:352t3kuC <a href="http://jkxotjkfkoqx.com/">jkxotjkfkoqx</a>

lzxslcbfdclzxslcbfdc2012/11/16 11:02lsMjkV , [url=http://pswxdpjynsds.com/]pswxdpjynsds[/url], [link=http://dfyzalsrkwwt.com/]dfyzalsrkwwt[/link], http://rcymvzaqhavf.com/

hzezocyuybrhzezocyuybr2012/11/17 21:21SC9H6Q , [url=http://pczvhtfxmuwm.com/]pczvhtfxmuwm[/url], [link=http://zqyfpyiuonso.com/]zqyfpyiuonso[/link], http://nudkivwjzsgg.com/

トラックバック - http://topcoder.g.hatena.ne.jp/jackpersel/20120420

2012-03-21SRM538 Div1

SRM538 Div1 250 EvenRoute

| 20:26 | SRM538 Div1 250 EvenRoute - SRM diary(Sigmar) を含むブックマーク はてなブックマーク - SRM538 Div1 250 EvenRoute - SRM diary(Sigmar) SRM538 Div1 250 EvenRoute - SRM diary(Sigmar) のブックマークコメント

Problem Statement

コーディングフェーズ

チェッカーボードで同じ色の間を移動するのなら、どんな経路を通っても偶数ステップになるんじゃないか

逆に違う色の間を移動するのなら、どんな経路を通っても奇数ステップになるんじゃないか

正しいっぽい

ということは原点と最後に訪れる点だけ決めればパリティが決まるし、どんな経路を通ってもいいのだから全ての点を通ってきたとすればいい

各点の原点からのマンハッタン距離のパリティを調べれば終わり

久々に結構すぐ解けた


ソースコード

class EvenRoute {
public:
	string isItPossible(vector <int> x, vector <int> y, int wp) {
		string res;
		int n=x.size();

		int par1=0, par0=0;
		for(int i=0; i<n; i++) {
			if((x[i]+y[i])%2==0) par0=1;
			else par1=1;
		}

		if(par0 && par1) return "CAN";
		if(wp==1 && !par1) return "CANNOT";
		if(wp==0 && !par0) return "CANNOT";
		return "CAN";
	}
};

SRM538 Div1 450 TurtleSpy

| 20:26 | SRM538 Div1 450 TurtleSpy - SRM diary(Sigmar) を含むブックマーク はてなブックマーク - SRM538 Div1 450 TurtleSpy - SRM diary(Sigmar) SRM538 Div1 450 TurtleSpy - SRM diary(Sigmar) のブックマークコメント

Problem Statement

コーディングフェーズ

直感的には、forwardとbackwardは全部結合して、forwardとbackwardの間に最適な角度を差し込めばいいように思われる

本当かな・・まあ450だし書いてみよう

適当なDPを書く

意外と時間がかかってしまった・・・実装力低すぎる・・・

とりあえずサンプル通ったんで提出

しかしアルゴリズム的に楽勝すぎる。本当に大丈夫か。これじゃあ450っていうより300くらいのような・・

forwardしてbackwardしてまたforwardするのがいいようなパターンはないのか

  • 一回目のforwardのベクトルとbackwardのベクトルを足すと、新たなベクトルができる。
  • 二回目のforwardのベクトルの方向が、一回目のforwardのベクトルと新たなベクトルの間に入るようにできれば、forwardを分割したほうがいい
  • どうやら検討してみると、そのような二回目のforwardを作れる回転があるなら、backwardする前に使ったほうがいい結果になるっぽい

なんだ。つまらない。しかし解くよりも証明するほうが大変そうな気がする。チキンレース気味だな。


ソースコード

ベクトルの回転にはpolarという便利な関数を使用する

まあ別に使わなくても大した差はないが・・・


あと、leftとrightを分けた辺りがちょっと無駄だった。rightはマイナスの値にしとけば良かったかな・・・

typedef complex <double> pt;

int la[52][370], ra[52][370];

class TurtleSpy {
public:
	double maxDistance(vector <string> com) {
		double res;
		int n=com.size();

		vector <int> l, r;
		int f=0, b=0;
		for(int i=0; i<n; i++) {
			stringstream ss(com[i]);
			string s;
			ss >> s;
			int v;
			ss >> v;
			if(s=="forward") f+=v;
			if(s=="backward") b+=v;
			if(s=="left") l.push_back(v);
			if(s=="right") r.push_back(v);
		}

		memset(la, 0, sizeof(la));
		memset(ra, 0, sizeof(ra));
		int lsz=l.size(), rsz=r.size();

		la[0][0]=ra[0][0]=1;
		for(int i=0; i<lsz; i++) {
			for(int a=0; a<360; a++) {
				if(!la[i][a]) continue;
				la[i+1][a]=1;
				int na=(a+l[i])%360;
				la[i+1][na]=1;
			}
		}
		for(int i=0; i<rsz; i++) {
			for(int a=0; a<360; a++) {
				if(!ra[i][a]) continue;
				ra[i+1][a]=1;
				int na=(a+r[i])%360;
				ra[i+1][na]=1;
			}
		}

		int aa[360];
		memset(aa, 0, sizeof(aa));
		for(int i=0; i<360; i++) {
			for(int j=0; j<360; j++) {
				if(la[lsz][i]&&ra[rsz][j]) aa[((i-j)+360)%360]=1;
			}
		}

		int besta=0;
		for(int i=0; i<360; i++) {
			if(aa[i]) {
				if(abs(180-i)<abs(180-besta)) besta=i;
			}
		}

		pt p(f, 0);
		pt ap=polar((double)b, M_PI*((besta+180)%360)/180.);
		p=p+ap;
		res=abs(p);
		
		return res;
	}
};

EsmeEsme2012/11/15 09:54Fell out of bed feeling down. This has brihgtened my day!

dfcfxtaudfcfxtau2012/11/17 05:21FAXLCr , [url=http://nigkynubluqk.com/]nigkynubluqk[/url], [link=http://sdzcawuwyqac.com/]sdzcawuwyqac[/link], http://eralavivvuvl.com/

fsakkadkfsakkadk2012/11/17 21:48Ln0gXw , [url=http://vgjurqseohhr.com/]vgjurqseohhr[/url], [link=http://gcpmrqxsqmmm.com/]gcpmrqxsqmmm[/link], http://mlzyuaxtrdqy.com/

トラックバック - http://topcoder.g.hatena.ne.jp/jackpersel/20120321

2012-03-17SRM537 Div1

SRM537 Div1 275 KingXNewCurrency

| 03:57 | SRM537 Div1 275 KingXNewCurrency - SRM diary(Sigmar) を含むブックマーク はてなブックマーク - SRM537 Div1 275 KingXNewCurrency - SRM diary(Sigmar) SRM537 Div1 275 KingXNewCurrency - SRM diary(Sigmar) のブックマークコメント

Problem Statement

コーディングフェーズ

よく分からんがmax(A, B)までの数iを候補として、適当にlcm(X, i)までの数で問題ないかチェックすればいいだろう

というのが正しそうかどうかを検討するのに20分くらい要した


後で

そもそもAとBをXとiで表せなければダメだし、逆にそれができれば他の数も全部表せるので、そのチェックだけで良かった

相変わらず不要な検討ばかりして時間がもったいない。。。


ソースコード

提出した版

ll gcd(ll a, ll b) {
	while(b) {
		a=a%b;
		swap(a, b);
	}
	return a;
}

ll lcm(ll a, ll b) {
	return a/gcd(a,b)*b;
}

class KingXNewCurrency {
public:
	int howMany(int A, int B, int X) {
		int res;
		
		if(A>B) swap(A, B);
		int g=gcd(A, B);
		if(g%X==0) return -1;

		res=0;
		for(int i=1; i<=B ;i++) {
			int g2=gcd(X, i);
			if(g%g2!=0) continue;
			int l2=lcm(X, i);
			int can[40010];
			memset(can, 0, sizeof(can));
			for(int j=0; X*j<=l2; j++) {
				for(int k=0; X*j+i*k<=l2; k++) {
					can[X*j+i*k]=1;
				}
			}
			bool ok=true;
			for(int j=0; A*j<=l2; j++) {
				for(int k=0; A*j+B*k<=l2; k++) {
					if(!can[A*j+B*k]) {
						ok=false;
						break;
					}
				}
				if(!ok) break;
			}
			if(ok) res++;
		}

		return res;
	}
};

SRM537 Div1 500 KingXMagicSpells

| 03:57 | SRM537 Div1 500 KingXMagicSpells - SRM diary(Sigmar) を含むブックマーク はてなブックマーク - SRM537 Div1 500 KingXMagicSpells - SRM diary(Sigmar) SRM537 Div1 500 KingXMagicSpells - SRM diary(Sigmar) のブックマークコメント

Problem Statement

コーディングフェーズ

期待値ってことは、また線形性か

XORだとビット毎に変化するので、ビット毎にそのビットが立つ確率を計算して、後で足しあわせればいいだろう

と思ったけど最大10億だったら無理か(←無理でない。30ビットしかない。)


うーん分からん。DPか。

DPも無理ぽい。お手上げ・・・


後で

チャレンジフェーズでようやく30ビットしかないことに気づいた

頭悪すぎる・・


ソースコード

class KingXMagicSpells {
public:
	double expectedNumber(vector <int> ducks, vector <int> s1, vector <int> s2, int K) {
		double res;
		int n=ducks.size();

		vector <vector <double> > p(n, vector <double>(30));
		for(int i=0; i<n; i++) {
			for(int j=0; j<30; j++) {
				if(ducks[i]&(1<<j)) p[i][j]=1;
			}
		}

		for(int i=0; i<K; i++) {
			vector <vector <double> > np(n, vector <double>(30));
			for(int j=0; j<n; j++) {
				for(int k=0; k<30; k++) {
					if(s1[j]&(1<<k)) np[j][k]=(1-p[j][k])/2;
					else np[j][k]=p[j][k]/2;
				}
			}
			for(int j=0; j<n; j++) {
				for(int k=0; k<30; k++) {
					np[s2[j]][k]+=p[j][k]/2;
				}
			}
			p=np;
		}
		
		res=0;
		for(int i=0; i<30; i++) {
			res+=p[0][i]*(1<<i);
		}
		return res;
	}
};
トラックバック - http://topcoder.g.hatena.ne.jp/jackpersel/20120317