ピタゴラス数の計算

ちょっとおもしろいアルゴリズムがあったので試してみる。求めるピタゴラス数の計算でのcの上限がつまり65535にすると一秒弱であっという間に計算できてしまうのがすごい。

/* ピタゴラス数 高速版*/

/*
	http://club.pep.ne.jp/~asuzui/page9.html
*/

#include <windows.h>

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* 求めるピタゴラス数 a^2+b^2=c^2 の c の上限 */
#define MAXNUMBER 0xffff

/* 変数を交換する */
#define SWAP(type,a,b) {type tmp=a;a=b;b=tmp;}

/* 最大公約数を計算する関数 */
int gcm(int x, int y)
{
	int z;
	
	/* x を y で割った余りを z に代入し z が 0 以外のあいだループする */
	while((z = x % y) != 0)  {
		x = y; /* xにyを代入し */
		y = z; /* yに余りzを代入する */
	}
	
	return y;
}

int main(int argc, char *argv[])
{
	int count = 0;
	int u = 0, v = 0, ulimit = 0;
	int a = 0, b = 0, c = 0;
	
	const int limit = MAXNUMBER;
	
	DWORD time_before = 0, time_after = 0;
	double time = 0;
	TCHAR buf[80];
	
	FILE *fp = NULL;
	
	remove("fast.txt");
	
	MessageBox(NULL, TEXT("ピタゴラス数の計算"), TEXT("これから計算を始めます"), MB_OK);
	time_before = GetTickCount();

	ulimit = sqrt(limit);
	
	fp = fopen("fast.txt", "a");
	
	for (u = 2; u <= ulimit; u++) {
		for (v = u - 1; v > 0; v -= 2) {
			
			c = u * u + v * v;
			
			if(c > limit) {
				continue;
			}
			
			if (gcm(u, v) != 1) {
				continue;
			}
			
			count++;
			a = 2 * u * v;
			b = u * u - v * v;
			
			if (a > b) {
				SWAP(int, a, b);
			}
			
			
			sprintf(buf, "%d:\ta = %d, b = %d, c = %d\n", count, a, b, c);
			fputs(buf, fp);
		}
	}
	
	time_after = GetTickCount();
	
	time = (double)(time_after - time_before)/1000.0;
	
	sprintf(buf, TEXT("処理時間は %.3lf 秒です"), time);
	
	fclose(fp);
	
	MessageBox(GetDesktopWindow(), buf, TEXT("ピタゴラス数の計算"), MB_OK);
	
	exit(EXIT_SUCCESS);
}