// File		pythagoras.C
// Description	Program to compute Pythagorean triples
// Author	Eric Conrad
// Usage	pythagoras pmax
//		where pmax is a (small) positive integer which parameterizes
//		the triples.
// Output	x y z [p,q]
//		where x^2 + y^2 = z^2 and p and q are the parameters which
//		determine x, y and z.
// Corrections
//	March 8, 2000 -- corrected an error in the swap code for the case
//		m<n in the gcd algorithm.  It doesn't affect the program
//		since the generated triples are subject to the constraint
//		p>q, but I'd hate to see someone use the gcd code in another
//		program...  (Thanks to my student Ted Pavlic for noting
//		this error.)  Also main now returns a 0 for the benefit
//		of those who have compilers that scream about these things.

#include <iostream.h>
#include <stdlib.h>

#define INT	unsigned long

INT gcd(INT, INT);		// Euclid's gcd algorithm (see below)

int main(int argc, char** argv)
	{
	if( argc != 2 )
		{
		cerr << "Usage: " << argv[0] << " pmax\n";
		exit(1);
		}

	INT pmax = atol( argv[1] );
	for( INT p=2 ; p <= pmax ; p++ )
	    for( INT q=1 ; q < p ; q++ )
		{
			/* We want p and q relatively prime and of opposite
			   parity so that we only obtain primitive triples */
		if( (p+q)%2 == 0 ) continue;
		if( gcd(p,q) > 1 ) continue;

		INT x = p*p - q*q;
		INT y = 2 * p * q;
		INT z = p*p + q*q;

		if( x < y )
			cout << x << " " << y;
		else
			cout << y << " " << x;
		cout << " " << z << "   [" << p << "," << q << "]\n";
		}

	return 0;
	}

	// Euclid's gcd algorithm  (very elegant and very ancient!)

INT gcd( INT m, INT n )
	{
		// To start everything must be non-negative and m>=n
	if( m<n )
		{	// There are clever machine-dependent ways to do swaps.
			// This is not one of them.
		INT swap;
		swap = m;		// ** corrected March 9, 2000 **
		m = n;
		n = swap;
		}

	while( n!=0 )
		{
		INT r = m % n;		// division algorithm
		m = n;			// set up next iteration
		n = r;
		}

	return m;
	}