素因数分解プログラム

/*		factor.c		*/
#include <stdio.h>
#include <math.h>

#define MAXCOUNT	500
#define MAXFRACT	200
#define MAXFRACT1	201		/* MAXFRACT + 1 */

long GCD(long m, long n);
void fract(long f[], long n);

int main(void)
{
	long f[MAXFRACT1];
	long n, *p;

	n = 691891200L;
	fract(f, n);
	p = f;
	printf("%ld = %ld ", n, *p);
	while(*(++p) != 0)	printf("x %ld ", *p);
	putchar('\n');
	return 1;
}

void fract(long f[], long n)
{
	long a0, p, pa, pb, b, bs, c, m, na;
	int count,i = 0,j;

	f[0] = 0L;
	if(n <= 0L)	return;
	do
	{
		a0 = (long)sqrt((double)n);
		pa = 0L;
		pb = 1L;
		m = a0;
		b = 1L;
		c = 0L;
		count = 0;

		while(count++ <= MAXCOUNT)
		{
			for(j = 1, p = pa; j <= m; j++)
			{
				p += pb;
				p %= n;
			}
			c = (m * b - c);
			b = (n - c * c) / b;
			m = (a0 + c) / b;
			bs = (long)sqrt((double)b);
			if((b == bs * bs) && (p - bs > 1) && (p - bs != n)){
			
				na = GCD(p - bs,n);
				if(na != 1L)
				{
					f[i++] = na;
					if(i == MAXFRACT)
					{
						f[i] = 0L;
						return;
					}
					n /= na;
					break;
				}
			}
			pa = pb;
			pb = p;
		}
		if(count > MAXCOUNT)
		{
			f[i++] = n;
			break;
		}
	}while(n > 1);
	f[i] = 0L;
	return;
}

long GCD(long m, long n)
{
	long w;

	if(m <= 0L || n <= 0L)	return 0L;
	while(n > 0L)
	{
		w = m % n;
		m = n;
		n = w;
	}
	return m;
}

691891200 = 26208 x 160 x 11 x 3 x 5