p94.py

import sys

# these are called isosceles Heronian triangles, and there's a generating function for them:
# https://en.wikipedia.org/wiki/Heronian_triangle

# considering that side lengths must be in the forms
# n, n, n+1 or 
# n, n, n-1

# these must be primitive Heronian triangles, 
# because for all n, gcd(n, n+1) = 1,
# and primitives are defined as gcd(a, b, c) = 1

# pairs of isosceles heronian triangles are generated by:
# a = 2(u**2 - v**2)
# b = u**2 + v**2
# c = u**2 + v**2

# a = 4uv
# b = u**2 + v**2
# c = u**2 + v**2

# for coprime u, v with u > v and u + v odd

def gcd(u, v):
    i = u + 0
    j = v + 0
    while i > 0 and j > 0:
        if i > j:
            i = i % j
        else:
            j = j % i
    return max(i, j)

print gcd(13, 1)
print gcd(12, 21)
print gcd(121, 22)
print gcd(9, 9)

def coprime(u, v):
    return gcd(u, v) == 1

print coprime(9, 8)
print coprime(6, 9)

total = 0

limit = 1000000000

i = 1

while True:
    # ensure i + j is odd
    if i % 2 == 0:
        j = 1
    else:
        j = 2
    # find heronian triangles with side-lengths n, n, n +/- 1
    while j < i:
        a = 2 * (i ** 2 - j ** 2)
        b = i ** 2 + j ** 2 
        c = i ** 2 + j ** 2
        a2 = 4 * i * j
        if 2 * i ** 2 + 4 * i > limit:
            print i
            sys.exit("limit reached")
        if coprime(i, j):
            if abs(a - b) == 1:
                print i, a, b, c
                total += a + b + c
                print total
            if abs(a2 - b) == 1:
                print i, a2, b, c
                total += a2 + b + c
                print total
        j += 2
    i += 1