Generating Pythagorean triples using Gaussian (complex) integers

Question

I have only just recently found out about a way of generating Pythagorean triples through this video explaining it, involving the use of Gaussian (complex) integers. I have so far managed to write a function returning a list of Pythagorean triples generated by every Gaussian integer where the imaginary part is less than the real part.

def pyt(max_real):
    t = []
    real = 2
    imag = 1
    while real <= max_real:
        z = complex(real, imag)**2
        t.append((z.real, z.imag, abs(z)))
        if imag + 1 == real:
            real += 1
            imag = 1
        else:
            imag += 1
    return t

The problem with this is that some triplets (such as {9, 12, 15}) are not generated through the initial step in the video that the function has been based on, and I'm unsure of how to generate those.

>>> for i in pyt(4):
        print(i)


(3.0, 4.0, 5.0)
(8.0, 6.0, 10.0)
(5.0, 12.0, 13.0)
(15.0, 8.0, 17.0)
(12.0, 16.0, 20.0)
(7.0, 24.0, 25.0)
>>> # missing: (9, 12, 15), possibly others

How would I go about generating every possible triplet, somehow using the ones I already have or otherwise?

Answer 1

Edit: I realized this could actually miss some triples, see the notice on the last line.

This answer is based on the link you provided. We will use the following information:

• We can use the method of Gaussian integers to find all triple generators

• Any triple is a multiple of one of the above generators

• To find a triple, we never need to scale a generator by less than 1/2, giving an upper bound to the biggest generator we will need.

So here is some pseudocode of how we can proceed. Some details on possible implementation will follow.

def pyt(max_real):
    max_generator = 2 * max_real
    generators = set()

    # Find every generator inside our upper bound
    for x in [Gaussian integers if abs(x) < max_generator and x.imag < x.real]:
        y = x**2
        triple = (y.real, y.imag, abs(y))
        generators.add(triple)


    # Scale up
    scaled_up_generators = set()

    for a, b, c in generators:

        for i in range(max_real / c):
            scaled_up_generators.add((i * a, i * b, i * c))

    # Scale down
    triples = set()

    for a, b, c in generators:

        common_factors = find_common_factors(a, b, c)
        for factor in common_factors:
            triples.add((a / factor, b / factor, c / factor))

    triples = set()

    # Before returning we filter out the triples that are too big.
    return filter(lambda triple: triple[2] <= max_real, triples)

The above code recovers all triple generators up to twice the provided bound. Then by scaling them up and down, we recover all triples inside our boundary.

You might want to have a look at an efficient way to find common factors to implement find_common_factors, here is a start.

Again, this implementation is based solely on the information provided in your link. Also it will catch way more triple, but might not catch triples that need to be scaled by finer fractions. There might be a more efficient way to proceed, but for that I recommend turning to MathExchange.

Answer 2

Consider this paragraph¹ from the Wikipedia page about Pythagorean triples:

Despite generating all primitive triples, Euclid's formula does not produce all triples—for example, (9, 12, 15) cannot be generated using integer m and n. This can be remedied by inserting an additional parameter k to the formula. The following will generate all Pythagorean triples uniquely:

a = k ⋅ (m² − n²), b = k ⋅ ( 2 ⋅ m ⋅ n ) , c = k ⋅ (m² + n² )

where m, n, and k are positive integers with m > n, and with m and n coprime and not both odd.

While the OP used complex, in the following implementation I avoid floating-point types.

from math import gcd as gcd

def pythagorean_triples(max_c):
    triples = []
    m = 2
    while m * m < max_c:
        # "... not both odd..."
        for n in range(1 + m % 2, m, 2):

            # "... with m and n coprime ..."
            if gcd(m, n) != 1:
                continue
            
            c = m * m + n * n
            if c > max_c:
                break

            a = m * m - n * n
            b = 2 * m * n
            if b < a:
                a, b = b, a            
            
            k = 1
            while c * k <= max_c:
                triples.append((a * k, b * k, c * k))
                k += 1
        m += 1
    return triples

Note that I also used the maximum value of the hypotenuse as limit, not the maximum real part of the Gaussian integers:

>>> for a, b, c in sorted(pythagorean_triples(25)):
       print(a, b, c)

3 4 5
5 12 13
6 8 10
7 24 25
8 15 17
9 12 15
12 16 20
15 20 25

^{1) https://en.wikipedia.org/wiki/Pythagorean_triple#Generating_a_triple}