fermat.py

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

import sys

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 usage_exit(a0, rc):
    sys.stderr.write("""
Find Pythagorean/Fermat triples upto some bound.
Usage:
   %s <Bound> <Power>
""" % a0)
    sys.exit(rc)


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


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

bound = int(sys.argv[1])
power = int(sys.argv[2])
a = 1
n_triples = 0
while a <= bound:
    b = a
    while b <= bound:
        c = hypotenuse_power(a, b, power)
        if c != -1:
            n_triples += 1
            sys.stdout.write("%5d %5d %5d\n" % (a, b, c))
        b += 1
    a += 1
triple_type_name = "Pythagorean"
if power > 2:
    triple_type_name = "Fermat"
sys.stdout.write("Found %d %s triples\n" % (n_triples, triple_type_name))
sys.exit(0)