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

https://projecteuler.net/problem=292
周長120までのピタゴラス多角形の種類数をこたえる問題。
数分の計算時間で答えが出るところまで到達したので満足。
多分エレガントな漸化式で解くのだろうけど、力技の漸化式でだましだまし解。
P(LIMIT)に設定しているのでコード中のLIMIT=30をLIMIT=120に変更すると答えを計算する。

#include<stdio.h>
#include<set>
#include<map>
#include<iostream>
int gcd ( int a, int b ){
 	int c;
 	while ( a != 0 ) {
 	   c = a; a = b%a;  b = c;
 	}
  	return b;
}

struct E{
	int dx,dy,len;
	bool operator<(const E& e1)const{
		if(dx!=e1.dx)return dx<e1.dx;
		return dy<e1.dy;
	}
};
 
struct E2{
	int x,y,len,dx,dy,sx,sy;
	bool operator<(const E2& e2)const{
		if(x!=e2.x)return x<e2.x;
		if(y!=e2.y)return y<e2.y;
		if(dx!=e2.dx)return dx<e2.dx;
		if(dy!=e2.dy)return dy<e2.dy;
		if(sx!=e2.sx)return sx<e2.sx;
		if(sy!=e2.sy)return sy<e2.sy;
		return len<e2.len; 
	}
};
 
bool ok(int x1,int y1,int x2,int y2){ 
	return x1*y2-x2*y1>0;
}
 
std::set<E> memo;
const int LIMIT=30;//とりあえずP(30)に設定ここを１２０に変えると数分の計算で問題の答えが出る
 
int main(){
 
	E e1;
	e1.len=1;
	e1.dx=-1;
	e1.dy=0;
	memo.insert(e1);
	e1.dx=1;
	e1.dy=0;
	memo.insert(e1);
	e1.dx=0;
	e1.dy=1;
	memo.insert(e1);
	e1.dx=0;
	e1.dy=-1;
	memo.insert(e1);
	for(int m=2;m*m+1<=LIMIT;m++){
		for(int n=1;n<m;n++){
 
			int c=m*m+n*n;
			if(c>LIMIT)break;
			if(gcd(m,n)!=1)continue;
			if((m-n)%2!=1)continue;
 
			e1.len=c;
 			e1.dx=m*m-n*n;
			e1.dy=2*m*n;
			memo.insert(e1);
			e1.dx=-e1.dx;
			memo.insert(e1);
			e1.dx=-e1.dx;
			e1.dy=-e1.dy;
			memo.insert(e1);
			e1.dx=n*n-m*m;
			e1.dy=-2*m*n;
			memo.insert(e1);
 
			e1.dy=m*m-n*n;
		e1.dx=2*m*n;
			memo.insert(e1);
			e1.dx=-e1.dx;
			memo.insert(e1);
			e1.dx=-e1.dx;
			e1.dy=-e1.dy;
			memo.insert(e1);
			e1.dy=n*n-m*m;
			e1.dx=-2*m*n;
			memo.insert(e1);
		}
	}
	printf("%d\n",memo.size());
	std::set<E>::iterator it;
	E2 e2;
	std::map<E2,long long int> ms;
	for(it=memo.begin();it!=memo.end();it++){
		e1=(*it);
		if(e1.dx>0&&e1.dy>=0){
 
			e2.dx=e1.dx;
			e2.dy=e1.dy;
			e2.sx=e1.dx;
			e2.sy=e1.dy;
			for(int i=1;i*e1.len<LIMIT/2;i++){
				e2.len=e1.len*i;
				e2.x=e1.dx*i;
				e2.y=e1.dy*i;
				ms[e2]++;
			}
		}
	}
	long long int ans=0;
	for(int j=2;j<=LIMIT;j++){
		std::map<E2,long long int> next;
		std::map<E2,long long int>::iterator mIt;
		long long int add=0;
		for(mIt=ms.begin();mIt!=ms.end();mIt++){
			e2=(*mIt).first;
			long long int c=(*mIt).second;
 			for(it=memo.begin();it!=memo.end();it++){
 
				e1=(*it);
 
				E2 e3; 
				for(int i=1;e2.len+i*e1.len<=LIMIT;i++){
					e3.x=e2.x+e1.dx*i;
					e3.y=e2.y+e1.dy*i;
					e3.len=e2.len+e1.len*i;
					e3.dx=e1.dx;
 					e3.dy=e1.dy;
					e3.sx=e2.sx;
					e3.sy=e2.sy;
					if(ok(e1.dx,e1.dy,e2.dx,e2.dy)==false)continue;
					if(e3.x==0&&e3.y==0){
 						add+=c;
					}else{
						if(ok(e3.x,e3.y,e2.x,e2.y)==false)continue;
						if(ok(e3.x,e3.y,e2.sx,e2.sy)==false)continue;
						int len2=e3.x*e3.x+e3.y*e3.y;
						int len3=(LIMIT-e3.len)*(LIMIT-e3.len);
 						if((len2>len3)||(len2>(LIMIT/2)*(LIMIT/2)))continue;
						next[e3]+=c;
					}
 				}
			}
		}
		std::cout<<add<<" "<<j<<" "<<ans<<" "<<next.size()<<"\n";
		add=add*4/j;
 		ans+=add;
		ms.clear();
		ms.insert(next.begin(),next.end());
 	}
	std::cout<<"\nans="<<ans<<"\n";
}