pythagorean3s.py

#!/usr/bin/env python
#
#  Find reduced Pythagorean triples upto some bound.
#  Author:  Yotam Medini  [email protected] -- Created: 2006/January/29
#

import sys

def gcd(m, n):
    if m < n:
        t = m;  m = n;  n = t; # swap
    # Now we have:  m >= n >= 0
    while n > 0: # Euclid
        remainder = m % n
        m = n
        n = remainder
    return m


def iroot(a, n):
    m = 0
    M = a + 1
    error = M - m
    progress = True
    while progress:
        x = (m + M)/2
        xPower = x**n
        if xPower == a:
            progress = False
        else:
            if xPower < a:
                m = x
            else:
                M = x
            new_error = M - m
            progress = (new_error < error)
            error = new_error
    return x


def hypotenuse(a, b):
    c2 = a*a + b*b
    c = iroot(c2, 2)
    if c*c != c2:
        c = -1
    return c


def usage_exit(a0, rc):
    sys.stderr.write("""
Find Pythagorean triples upto some bound.
Usage:
   %s <Bound>
""" % a0)
    sys.exit(rc)


if len(sys.argv) != 2:
    usage_exit(sys.argv[0], 1)

bound = int(sys.argv[1])
a = 1
n_triples = 0
while a <= bound:
    b = a
    while b <= bound:
        c = hypotenuse(a, b)
        if c != -1 and gcd(a, b) == 1:
            n_triples += 1
            sys.stdout.write("%5d %5d %5d\n" % (a, b, c))
        b += 1
    a += 1
sys.stdout.write("Found %d reduced Pythagorean triples\n" % n_triples)
sys.exit(0)