pythagorean.py

#! /usr/bin/env python

import heapq

def gcd (m, n):
    while 0 < n:
        (m, n) = (n, m % n)
    return m

def pythagorean_triples ():
    triples = [(5, 1, 2, 1)]
    while True:
        (c, m, n, k) = heapq.heappop(triples)
        if (1 == gcd(m, n)):
            yield (k*n*n - k*m*m, 2*k*m*n, c)
            heapq.heappush(triples, ((k+1)*n*n + (k+1)*m*m, m, n, k+1))
        if k == 1:
            if m+1 == n:
                heapq.heappush(triples, ((n+1)*(n+1) + n*n, m+1, n+1, 1))
            heapq.heappush(triples, ((n+2)*(n+2) + m*m, m, n+2, 1))

for t in pythagorean_triples():
    if 1000 < t[2]:
        break
    print(t)