#!/usr/bin/env python # # Find reduced Pythagorean triples upto some bound. # Author: Yotam Medini yotam.medini@gmail.com -- 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 """ % 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)