Finding Primes Using Complex Numbers

Pythagoras In LegoWith complex numbers, I always feel as if I’m getting a glimpse of something truly awesome that lies hidden within mathematics. The first time I understood how they worked, I thought it was some form of magic.

I get the same feeling with prime numbers. Like many, I’ve looked at them from all angles – prime gaps, large primes, prime densities, prime sieves – and they continue to fascinate. A few months ago I was thumbing through Henry Warren’s programmers cookbook Hacker’s Delight (A) and discovered a whole chapter on the various formulas for (some of) them. Mind-bending stuff.

Then a few weeks ago, after reading Simon Singh’s fascinating book Fermat’s Last Theorem (B), and having a go at solving it (as you do), I went back to Pythagorean triplets, wondering how many there were and what the largest one every found was, i.e. the largest Pythagorean triplet that is not just a multiple of a smaller one.

And then I remembered something I’d learned about primes…and something clicked. OK, I’d had a little coffee, so I’m not sure what the thought sequence was, but here is an attempt to piece it all together with a series of assertions.

Pythagoras’ Theorem

1. Pythagoras tells us how to calculate the hypotenuse of a right-angle triangle: the sum of the squares of the two shorter sides is equal to the square of the hypotenuse.

Pythagoras triangles

2. When all the sides are whole numbers, we have a Pythagorean triplet. The smallest of these is a 3-4-5 triangle. (see top image).

 

3. But there are many sums of 2 squares that don’t total to another square, but instead to a integer which is prime, e.g:

1 & 16 (=17), 1 & 36 (=37), 1 & 100 (=101), 1 & 196 (=197), etc…,

4 & 9 (=13), 4 & 25 (=29), 4 & 49 (=53), 4 & 169 (=173), etc…,

9 & 64 (=73), 9 & 100 (=109), 9 & 400 (=409), etc…

 

4.  Since each of the numbers in the above pairs are squares, their square roots will be whole numbers. But since the hypotenuses are prime, their square roots will be irrational, e.g:

(1, 16, 17) => (1, 4, √17), (4, 9, 13) => (2, 3, √13), (9, 64, 73) => (3, 8, √73).

Complex Number Multiplication

5.  The result of multiplying 2 complex numbers together, where each component is whole, is always another complex number with whole number components, e.g:

(2 + 3j)(3 + j) = 6 + 2j + 11j – 3 = 3 + 11j

 

6.  If we plot a complex number, superimposing both its Cartesian and polar forms, we obtain a right-angle triangle whose hypotenuse is the polar length of the complex number:

 

Argand diagram showing Pythagoras plotted as a complex number

 

7.  The equivalent in polar terms of multiplication is to add the angles and multiply the lengths. e.g:

2 ∠45° x 3 ∠25° = 6 ∠70°

 

8.  If we square a complex number, we therefore double the angle and square the length:

Plotting the square of a complex number

 

9.  We can therefore plot any of the square-root-triplets listed above in (4) in the complex plane, e.g:

a. The (2, 3, √13) triangle becomes the complex number 2 + 3j.
and
b. The (3, 8, √73) triangle becomes the complex number 3 + 8j.

 

10. From Assertion (3), if we square any of the complex numbers found this way, the result will have a hypotenuse whose length is a prime number, regardless of its new angle or the lengths of its sides.

 

11.  From Assertion (5), the real and imaginary components will also be whole numbers, meaning that the three corners of the squared complex number will also fall exactly onto the complex cartesian grid.

 

12.  Since the hypotenuse will now be a prime number, the lengths of the three sides will be irreducible to a smaller Pythagorean triplet (in the way that 6-8-10 is reducible to a 3-4-5 triplet). Rather, they will now form a new, much larger right-angle triangle whose sides form a perfect Pythagorean triplet.

 

By this method, Pythagorean triplets can easily be found by simply checking for sums of squares that are prime, and by applying some elementary arithmetic on complex numbers.

Worked Examples

Example 1:

4 squared (16), plus 1 squared (1) = 17, which is prime.

The complex number (4 + j) has a hypotenuse of √17, therefore the square of this complex number has a hypotenuse that is 17 long.

Squaring it, we get (4 + j)(4 + j) = 16 + 8j -1 = 15 + 8j. We have therefore discovered (8, 15, 17) as a perfect Pythagorean triplet, the elements of which form the lengths of a right-angle triangle.

 

Example 2:

2 squared (4), plus 5 squared (25) = 29, which is prime.

The complex number (5 + 2j) has a hypotenuse of √29, therefore the square of this complex number has a hypotenuse that is 29 long.

Squaring it, we get (5 + 2j) (5 + 2j) = 25 + 10j – 4 = 21 + 10j. We therefore know that (10, 21 ,29) form a perfect Pythagorean triplet, the elements of which form the lengths of a right-angle triangle.

 

Example 3:

2 squared (4), plus 7 squared (49) = 53, which is prime.

The complex number (7 + 2j) has a hypotenuse of √53, therefore the square of this complex number has a hypotenuse that is 53 long.

Squaring it, we get (7 + 2j) (7 + 2j) = 49 + 28j – 4 = 45 + 28j. We therefore know that (28, 45, 53) forms a perfect Pythagorean triplet, the elements of which form the lengths of a right-angle triangle.

 

So What Is Happening?

What I seem to have done here is stumble upon a different method to that discovered by Greek mathematician Euclid of Alexandria (active c. 300 BC) for calculating new Pythagorean triples. By treating Pythagorean triples as complex numbers, and carefully choosing their sizes, new Pythagorean triples are revealed. (I’m pretty sure it’s been done before – this is just about me having fun and working it out for myself!)

Historically, Norwegian-Danish surveyor & mathematician Caspar Wessel (1745–1818) was the first to propose in 1799 a method of treating complex numbers as vectors on a plane. Plotting a complex number on a Cartesian plane with real and imaginary axes was later later called an Argand diagram, after another complex number pioneer, the amateur Swiss mathematician Jean-Robert Argand (1768–1822).

Regardless of who came up with this method first, we can be fairly certain it never occurred to Euclid, since the Greeks knew nothing of complex numbers, as far as we know. But what he did do was describe a method for calculating an infinite number of Pythagorean triples by starting with any two numbers, m, n, where one is odd and one is even. He asserted that if they are chosen so that they are co-prime (i.e. sharing no factors and neither being a multiple of the other, so that they cannot be reduced to a smaller proportional pair) then a primitive (irreducible) Pythagorean triple (x, y, z) exists such that:

x = m² – n²

y = 2mn

z = m² + n²

Now, if we take Wessel’s idea of treating a complex number C = (a + bj) as a vector, its square C² is calculated as follows:

C² = (a + bj)²

= (a + bj)(a + bj)

= a² – b² + 2abj

= (a² – b²) + j(2ab)

Which just happens to correspond to Euclid’s assertion. Which I think is pretty damned cool.

It would be nice if Euclid’s hypotenuse z was always prime, which would give us a great way to calculate massive primes, but it’s not (e.g. m=8, n=1, give a (16, 63, 65) triangle). In fact, as Euclid’s a and b tend to infinity, only about 40% of his hypotenuses are prime. I know, because I tested it in Python up to a=b=500,000, with this code:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Finds all of Euclid's Pythagorean triplets that result in a prime hyptotenuse,
up to triangles of sides 500 units long. Yields primes up to 500,000.
 
Created on Sat Sep  9 07:18:09 2017
 
author: matta_idlecoder at protonmail dot com
"""
 
from math import sqrt
import sys
import timeit
 
def is_prime(Number, return_first_factor=False):
    """Returns True if Number is prime, else False
 
    Its job is not to find all the factors, only to see if a factor exists
    """
    if type(Number) != int:
        print ("Error in function is_prime() - it was passed a non-int.")
        sys.exit(0)
 
    elif Number < 2:
        print ("Error in function is_prime() - it was passed a number less than 2.")
        sys.exit(0)
 
    elif (Number == 2) or (Number == 3):
        return (True, Number) if return_first_factor else True
 
    elif Number % 2 == 0:
        return (False, 2) if return_first_factor else False
 
    elif Number % 3 == 0:
        return (False, 3) if return_first_factor else False
 
    # if you're here, Number >= 5, an int, odd, and not a multiple of 3:
    go_no_further = int(round(sqrt(Number), 0))
 
    for poss_factor in range(3, go_no_further+2, 2):
        if Number % poss_factor == 0:
            # finds first factor of Number, don't care which:
            return (False, poss_factor) if return_first_factor else False
 
    return (True, 1) if return_first_factor else True
 
 
def get_prime_list_to_half(number):
    """returns a set of the primes up the square root of number (if whole)
    """
    prime_list = []
    if number > 2:
        prime_list.append(2)
 
    for poss_prime in range(3, int(number/2)+1, 2):
        if is_prime(poss_prime):
            prime_list.append(poss_prime)
 
    return prime_list
 
 
def co_prime(x, y, primes):
    """Returns True if x and y share any factors
 
    Algoirthm is first to check if one is not a factor of the other, then
    to see that they share no factors from the list of primes.
    """
    if x == y:
        print ("Error in function co_prime() - has been passed the same 2 numbers.")
        sys.exit(0)
    elif x % y == 0 or y % x == 0:
        return False
 
    for prime in primes:
        if (x % prime == 0) and (y % prime == 0):
            return False
    return True
 
 
 
def main():
    """Main function
    """
    Max_N = 500
 
    prime_list = get_prime_list_to_half(Max_N)
 
    SumSquares = []
    for i in range (2, Max_N+1, 1):
        for j in range (i, Max_N+1, 1):
 
            if i == j:
                continue
 
            if (i % 2 == 1) and (j % 2 == 1):
                # if i and j are both odd:
                continue
 
            if (i % 2 == 0) and (j % 2 == 0):
                # if i and j are both even:
                continue
 
            if not co_prime(i, j, prime_list):
                continue
 
            SumSquares.append([i, j, (i**2 + j**2)])
 
    SumSquares_copy = SumSquares[:]  
    for sq_set in SumSquares_copy:
        sum_sq = sq_set[2]
        if sum_sq % 5 == 0:
            SumSquares.remove(sq_set)
        elif (sum_sq - int(sqrt(sum_sq))**2) == 0:
            # remove any perfect squares:
            SumSquares.remove(sq_set)
 
    MaxPrime = 1
    non_prime_count  = 0
    SumSquares_copy = SumSquares[:] 
    for index, sq_set in enumerate(SumSquares_copy):
        if not is_prime(sq_set[2]):
            non_prime_count += 1
            SumSquares[index].append([is_prime(sq_set[2], return_first_factor=True)[1]])
        else:
            MaxPrime = sq_set[2] if sq_set[2] > MaxPrime else MaxPrime
 
    square_total = len(SumSquares)
    prime_total = square_total - non_prime_count
    print ("\nCounting up to {}, {:,d} out of the {:,d} Euclid's co-prime square-sums are prime, or about {:.0%}.".
           format(Max_N, prime_total, square_total, prime_total/square_total))
 
    print ("\nThe largest prime found was {:,d}".format(MaxPrime))
    return
 
 
if __name__ == "__main__":
    start_time = timeit.default_timer()
    main()
    stop_time = timeit.default_timer()
    prime_time = stop_time - start_time
    print ('\nThat took', round(prime_time, 1), 'seconds.')

 

One problem with this method is that it doesn’t find all the primes (e.g. 11 isn’t the sum of any coprime squares). Nor is it a very fast way of calculating primes. To illustrate, using the above python code on my laptop, here’s a comparison of this method for finding primes up to 500,000, against an Eratosthenes sieve method. First, the timing using Euclid’s coprime method:

Counting up to 500, 13,654 out of the 33,658 Euclid's co-prime square-sums are prime, or about 41%.

The largest prime found was 495,017

That took 18.1 seconds.

And here’s the timing for running an Eratosthenes sieve in Python, finding all primes up to 500,000:

Enter the number you want to find all the primes up to: 500000

Finding all primes up to 500,000................................................
Completed the search in 0 mins 0.11 seconds.

Saving prime data set to ./KnownPrimes.db..........

prime_data['MaxPrimeFound']: 499,979
prime_data['MaxPrimeUsedToSieve']: 701
prime_data['NumberOfPrimes']: 41,538
prime_data['MaxPrimeGapFound']: 114

Which makes the Eratosthenes sieve method about 165 times faster, and that was finding all the primes.

Whether or not I’m on to something here, I’m afraid my maths isn’t good enough (yet) to know where to go with this. Anyone interested in looking further into this should probably first check out Wikipedia’s page on the Pythagorean triples, which tells us that “Primitive Pythagorean triples have been used in cryptography as random sequences and for the generation of keys”.(C)

On the other hand, that was so enjoyable, I think I might take the rest of the week off.

 

NOTES
(A) Henry Warren, Hacker’s Delight, Addison-Wesley, Boston, 2013

(B) Simon Singh, Fermat’s Last Theorem, Fourth Estate, London, 1997

(C) Kak, S. and Prabhu, M. Cryptographic applications of primitive Pythagorean triples. Cryptologia, 38:215–222, 2014