ピタゴラス数

2chの宿題スレにて、3辺の和が1000以下の直角三角形の総数とそれらの中で面積が最大になる辺を出力という問題があった。
高々1000程度なら3重ループの総当たりでもスペックに任せて答えが出せるけど、10000以上とかになるともうだめ。
高速化するために、原始ピタゴラス数を求める方法を利用する。
m,nが互いに素であり、一方が偶数のとき、
a = k * (|m*m - n*n|)
b = k * 2 * m * n
c = k * (m*m + n*n)
で得られる。k=1のとき、a,b,cは原始ピタゴラス数と呼ぶ。
全てのピタゴラス数は原始ピタゴラス数に正整数kを掛けることで得ることが出来る。


重複を避けるために、n>mとする。
最大公約数gcdを使って素か否かを判定。
mとnでループを廻して、互いに素かつ一方が偶数のときにk倍してピタゴラス数を生成すればおk。
うっかり「一方が偶数のとき」というのを忘れると、互いに素なm,nとm',n'でk倍したときに重複解が出てくる。
これでめっちゃ速くなるお

#include <stdio.h>
#include <stdlib.h>
#define MAX 1000
int gcd(int x,int y)
{
	return y?gcd(y,x%y):x;
}
int main()
{
	int m,n,cnt=0;
	int a,b,c;
	double maxArea=0.0;

	for(m=1;m*m+2*m<MAX;m++){
		for(n=m+1;m*m+n*n+2*m*n+abs(m*m-n*n)<=MAX;n++){
			if(gcd(m,n)==1){
				int aa = abs(m*m-n*n);
				int bb = 2*m*n;
				int cc = m*m+n*n;
				int k;
				if(m%2 && n%2)continue;
				for(k=1;k*(aa+bb+cc)<=MAX;k++){
					cnt++;
					if(maxArea<(double)k*k*aa*bb){
						maxArea = (double)k*k*aa*bb;
						a=k*aa,b=k*bb,c=k*cc;
					}
					//printf("%d,%d,%d\n",aa*k,bb*k,cc*k);
				}
			}
		}
	}
	printf("cnt:%d\narea:%lf (%d,%d,%d)\n",cnt,maxArea/2.0,a,b,c);	
	return 0;
}