プロジェクトオイラー問329 - sinapusu2002 @ ウィキ

http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20329
Problem 329 「素数カエル」 †

Susan は素数カエルを飼っている.
このカエルは, 1 から 500 までの番号がつけられた 500 個の正方形の上を跳びまわる. 左右のいずれかの正方形に等確率で跳ぶことができ, 範囲 [1;500] の外に跳ぶことはできない.
（いずれかの端に着地したなら, 次は自動的に移動可能な正方形に跳ぶ. ）
素数である数が書かれた正方形にいる場合, 次の正方形に跳ぶ直前に, カエルは 2/3 の確率で 'P' (PRIME) と鳴き, 1/3 の確率で 'N' (NOT PRIME) と鳴く.
素数でない数が書かれた正方形にいる場合, 次の正方形に跳ぶ直前に, カエルは 1/3 の確率で 'P' と鳴き, 2/3 の確率で 'N' と鳴く.
カエルの開始位置は全ての正方形に対し等確率でランダムであり, また最初の 15 回の鳴き声を聞いたとすると, PPPPNNPPPNPPNPN と聞こえる確率は何か.
約分済みの分数 p/q の形で答えを入力せよ.

両端はそれ以外のマスの２倍として、漸化式で度数分布を数えて比率を計算するだけです。
計算時間数十秒。

#include<stdio.h>
#include<map>
#include<iostream>
const int LIMIT=500;
bool isPrime[LIMIT+1];
__int64 gcd( __int64 a, __int64 b )
{
 __int64 c;
 while ( a != 0 ) {
    c = a; a = b%a;  b = c;
  }
  return b;
}

void set_prime(){
	isPrime[1]=false;
	for(int n=2;n<=LIMIT;n++){
		bool isP=true;
		for(int i=2;i*i<=n;i++){
			if(n%i==0)isP=false;
		}
		isPrime[n]=isP;
	}
} 


int main(){
 	std::map<int,__int64> memo[LIMIT+1];
	set_prime();
for(int i=1;i<=500;i++){
		memo[i][0]=1;
	}

 	
	for(int i=0;i<14;i++){
		std::map<int,__int64> next[LIMIT+1];
		std::map<int,__int64>::iterator it;
		int addPerm=(1<<i);
		for(int j=1;j<=500;j++){
			for(int k=-1;k<=1;k+=2){
				int p=j+k;
				if(p<1||p>LIMIT)continue;
				for(it=memo[p].begin();it!=memo[p].end();it++){
					int perm=(*it).first;
					__int64 c=(*it).second;
					if(isPrime[p]==true){
 						next[j][perm+addPerm]+=c*2;
						next[j][perm]+=c;
				}else{
 						next[j][perm+addPerm]+=c;
						next[j][perm]+=2*c;
						if(p==1||p==500){
							next[j][perm+addPerm]+=c;
							next[j][perm]+=2*c;
						}
					}
				}
			}
		}
 		for(int j=1;j<=LIMIT;j++){
			memo[j].clear();
			memo[j].insert(next[j].begin(),next[j].end());
		}
	}
	
	int ansSeed[15]={1,1,1,1,0,0,1,1,1,0,1,1,0,1,0};
	int ansPerm=0;
 	for(int i=0;i<15;i++)ansPerm+=(ansSeed[i]<<i);
	__int64 all=0,u=0;
	std::map<int,__int64>::iterator it;
	for(int i=1;i<=500;i++){	
		for(it=memo[i].begin();it!=memo[i].end();it++){
			all+=(*it).second*3;
			if(ansPerm==(*it).first){
				if(isPrime[i]==true){
 					u+=(*it).second;
				}else{
					u+=(*it).second*2;
				}
			}
 		}
	}
	__int64 d=gcd(u,all);
	std::cout<<"\n"<<d<<" ans="<<u/d<<"/"<<all/d<<"\n";
}