Goldbach's conjecture algorithm

Question

I did a Goldbach's Conjecture exercise and got it to work. It's pretty slow though and I was wondering how I could optimize it.

number = int(input("Enter your number >> "))
print("\nCalculating...")
if number % 2 == 0: #Only even numbers
    primes = primenums(number) #returns all prime numbers <= input number
    addend1 = primes[0]
    addend2 = primes[0]

    while addend1 + addend2 != number:
        if primes[-1] == addend2:
            addend2 = primes[primes.index(addend1) + 1]
            addend1 = primes[primes.index(addend1) + 1]
        else:
            addend2 = primes[primes.index(addend2) + 1]

Right now, up to 10,000 the algorithm is pretty fast, but at 100,000 it takes about 3 seconds to finish. Is that just how it is or could I make it faster?

Answer 1

One thing that makes your code slow is your repeated calls to primes.index. Each of these calls needs to linearly scan primes until it finds the value you are looking for, making this \$\mathcal{O}(n)\$. This is especially bad in the first if case since you do this twice with exactly the same argument. Instead just keep the index around in addition to the number:

def goldbach_keep_index(number):
    if number % 2 == 1: #Only even numbers
        raise ValueError("Goldbach conjecture is only defined for even numbers")
    primes = primenums(number) #returns all prime numbers <= input number
    i, j = 0, 0
    addend1 = primes[i]
    addend2 = primes[j]

    while addend1 + addend2 != number:
        if addend2 == primes[-1]:
            i = j = i + 1
            addend2 = primes[i]
            addend1 = primes[i]
        else:
            j += 1
            addend2 = primes[j]
    return addend1, addend2

Note that I made this a function, so you can reuse it.

Your main calling code could look like this:

if __name__ == "__main__":
    number = int(input("Enter your number >> "))
    p1, p2 = goldbach_keep_index(number)
    print(f"{p1} + {p2} = {number}")

---

However, what you are really doing is taking all combinations of prime numbers, without repeating yourself. So just use itertools.combinations_with_replacement:

from itertools import combinations_with_replacement

def goldbach_combinations(number):
    if number % 2 == 0:
        raise ValueError("Goldbach conjecture is only defined for even numbers")
    primes = primenums(number)
    for addend1, addend2 in combinations_with_replacement(primes, 2):
        if addend1 + addend2 == number:
            return addend1, addend2
    raise Exception(f"Found a counter-example to the Goldbach conjecture: {number}")

In this case primes does not even need to be a list, it can just be a generator.

---

If you instead make primes a set, you can easily implement the idea suggested by @Josiah in their answer by just checking if number - p is in primes:

def goldbach_set(number):
    if number % 2 == 0: #Only even numbers
        raise ValueError("Goldbach conjecture is only defined for even numbers")
    primes = set(primenums(number)) #returns all prime numbers <= input number
    for p in primes:
        k = number - p
        if k in primes:
            return p, k
    raise Exception(f"Found a counter-example to the Goldbach conjecture: {number}")

---

And now some timing comparisons, where goldbach is your code put into a function:

def goldbach(number):
    if number % 2 == 0: #Only even numbers
        raise ValueError("Goldbach conjecture is only defined for even numbers")
    primes = primenums(number) #returns all prime numbers <= input number
    addend1 = primes[0]
    addend2 = primes[0]

    while addend1 + addend2 != number:
        if primes[-1] == addend2:
            addend2 = primes[primes.index(addend1) + 1]
            addend1 = primes[primes.index(addend1) + 1]
        else:
            addend2 = primes[primes.index(addend2) + 1]
    return addend1, addend2

And primenums is a simple sieve:

def prime_sieve(limit):
    prime = [True] * limit
    prime[0] = prime[1] = False

    for i, is_prime in enumerate(prime):
        if is_prime:
            yield i
            for n in range(i * i, limit, i):
                prime[n] = False

def primenums(limit):
    return list(prime_sieve(limit))

And when pulling the generation of the primes outside of the function (and calculating them up to the largest number in the plot):

Of course the first three functions are vastly slower than the set because they need to go through more combinations now. However, all functions are then constant time, so the increase in time comes solely from the fact that you have to consider more primes.

Answer 2

As suggested by Graipher in the comment, you're probably spending a lot of time generating a list of primes. That would be worth profiling, but I can't speculate as to how it might be improved without at least seeing the code.

There are two things that stand out to me as very clear performance sinks.

The first, and easiest to change, is index. If you call index on a list, python has to check every element in the list one at a time to see whether it's the thing you're after. If you have millions of primes, every one of those index calls is secretly a hidden loop that runs millions of times. The good news is you can resolve this easily enough: just keep another variable which stores where in the list of primes you're up to. Then you don't have to search for your place again every time you go round the loop, because you have a bookmark to it.

Second, and slightly more subtle, is the way that this is looping. You are looking for two numbers that add up to a target. If your target is 100, you might first say "what if it's 3?" and then try all primes with 3. Then you say "what if it's 5?" and try all the primes with 5. And then you do the same with 7, and so on. This means your trying all the primes with all the primes.

Instead, if you are checking whether the first prime is 3, you know that the second prime has to be 97. So all you have to do is check whether 97 is prime. If so, it's a match, if not, you can move on to checking 5 and 95, and then 7 and 93, and so on. In this way your algorithm checks each prime against one other number instead of each prime against all the others, and should be massively faster.

Answer 3

Let's take n = 1,000,000,000 for example and see what your code does.

It calculates all primes from 1 to 1 billion. That takes a while, but it gives you an array of all primes in sorted order.

It then calculates 2 + 2, 2 + 3, 2 + 5, 2 + 7, 2 + 999,999,xxx to check if one of these numbers equals 1,000,000,000. But obviously when addend1 = 2, addend2 has to be 999,999,998 to add to one billion, so you are checking tens of millions of sums unnecessarily.

It then calculates 3 + 3, 3 + 5, 3 + 7 etc., and again, addend2 would have to be 999,999,997 million. And then again for addend1 = 5, 7, 11 etc. So you are checking a huge number of sums needlessly.

Start with addend1 = first prime, addend2 = last prime in your array. If their sum is too small, replace addend1 with the next prime to make the sum larger by the smallest possible amound. If their sum is too large, replace addend2 with the previous prime. If the sum is right, you have found a solution. And once you reached addend1 > addend2, you know there is no solution. This will be very quick, since usually you don't need to check too many values for addend1.

That makes the search a lot quicker, but doesn't help with the sieve trying to find all the primes. You usually don't need all the primes, only a very small number. In the example with n = 1,000,000,000 you probably find two primes adding up to a billion with addend1 ≤ 1000 and addend2 ≥ 999,999,000.

So here's what you do: To find all primes say in the range 999,999,000 to 1,000,000,000, you need to check if these numbers are divisible by any prime up to the square root of 1,000,000,000 which is a bit less than 32,000. So first you use a sieve to find all numbers up to 32,000. Then you create a sieve that finds all primes from 999,999,000 to 1,000,000,000. You run your search algorithm until it tries to examine primes that are not in this range. This likely doesn't happen, but if it happens, you create another sieve for the numbers 999,998,000 to 999,999,000 and so on. So instead of maybe 50 million primes, you look only for maybe 50 or 100. Again, that makes it a lot faster.

Answer 4

This will be much faster. (less than a second to check 44 trillion) It does fewer prime checks, usually less than 100, and does not check numbers that are even for prime. Building a sieve of 700k primes for goldbach(12,000,000) takes forever. This code makes 7 IsPrime() calls for 12 million.
With native Python Prime code it takes about .001 second to find answer for 1 Billion making 42 calls to IsPrime().

Three major differences from other IsPrime(n-p) programs:

Checks Goldbach prime 3 then uses 6k ± 1 for the rest. So 1/3 of the Goldbach checks vs every odd number > 3.

More efficient IsPrime() tests using every 6th odd number using 6k ± 1 as well.

Can use GMPY2 to check for prime. (takes less than .001 second for 25 digit number.

Solutions using a sieve get n*k slower as the number size (k) grows. 10,000,000 takes 2.5 seconds 100,000,000 takes 25 seconds and so on until memory is exhausted building the sieve. The call prime_check only as needed grows linearly (k) as k grows: 1,000,000,000,000 takes 0.037 seconds while 100,000,000,000,000 takes 0.450 seconds.

With GMPY2:

$ python3 goldbach.py 12345678912345678
This is a test of Goldbach's Conjecture that for all even integers
greater than 2 there are two primes that add up to that even number.


61 + 12,345,678,912,345,617 = 12,345,678,912,345,678

         42 function calls in 0.000 seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    0.000    0.000 <string>:1(<module>)
        1    0.000    0.000    0.000    0.000 goldbach.py:16(goldbach)
        1    0.000    0.000    0.000    0.000 {built-in method builtins.exec}
        1    0.000    0.000    0.000    0.000 {built-in method builtins.print}
       37    0.000    0.000    0.000    0.000 {built-in method gmpy2.gmpy2.is_prime}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}


$ python3 goldbach.py 4444444444444444444444
This is a test of Goldbach's Conjecture that for all even integers
greater than 2 there are two primes that add up to that even number.


257 + 4,444,444,444,444,444,444,187 = 4,444,444,444,444,444,444,444

         144 function calls in 0.000 seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    0.000    0.000 <string>:1(<module>)
        1    0.000    0.000    0.000    0.000 goldbach.py:16(goldbach)
        1    0.000    0.000    0.000    0.000 {built-in method builtins.exec}
        1    0.000    0.000    0.000    0.000 {built-in method builtins.print}
      139    0.000    0.000    0.000    0.000 {built-in method gmpy2.gmpy2.is_prime}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

With my Home Grown IsPrime:

    $ python3 goldbach.py 12345678912345678
This is a test of Goldbach's Conjecture that for all even integers
greater than 2 there are two primes that add up to that even number.


61 + 12,345,678,912,345,617 = 12,345,678,912,345,678

         42 function calls in 3.869 seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    3.869    3.869 <string>:1(<module>)
        1    0.000    0.000    3.869    3.869 goldbach.py:16(goldbach)
       37    3.869    0.105    3.869    0.105 goldbach.py:6(IsPrime)
        1    0.000    0.000    3.869    3.869 {built-in method builtins.exec}
        1    0.000    0.000    0.000    0.000 {built-in method builtins.print}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

$ python goldbach.py 1000000000000000
This is a test of Goldbach's Conjecture that for all even integers
greater than 2 there are two primes that add up to that even number.


11 + 999,999,999,999,989 = 1,000,000,000,000,000

         14 function calls in 1.101 seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    1.101    1.101 <string>:1(<module>)
        1    0.000    0.000    1.101    1.101 goldbach.py:20(goldbach)
        3    0.000    0.000    0.000    0.000 goldbach.py:28(<genexpr>)
        6    1.101    0.184    1.101    0.184 goldbach.py:7(IsPrime)
        1    0.000    0.000    1.101    1.101 {built-in method builtins.exec}
        1    0.000    0.000    0.000    0.000 {built-in method builtins.print}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

Listing:

# Goldbach's Conjecture tester. 
from gmpy2 import is_prime
import sys
import cProfile

# or use home grown IsPrime
def IsPrime(n):  
   if (n <= 3): return n > 1      
   if not (n%6 == 1 or n%6 == 5): return False      
   for p,q in ((i,i+2) for i in range(5,int(n**.5)+1,6)): 
      if (n%p == 0 or n%q == 0):  return False          
   return True  


def goldbach(number):
    if number == 4:
        print("\n2 + 2 = 4\n")
        return
    elif IsPrime(number - 3):
        print(f"\n3 + {number-3:,} = {number}\n")
        return  
    else:
        for p in range(5, number, 6): # just 6k±1 after 3 
            if IsPrime(p) and IsPrime(number-p): 
                print(f"\n{p:,} + {number-p:,} = {N:,}\n")
                return
            elif IsPrime(p+2) and IsPrime(number-(p+2)):
                print(f"\n{p+2:,} + {number-(p+2):,} = {N:,}\n")
                return
        raise Exception(f"Found a counter-example to the Goldbach conjecture: {number}")

if __name__=="__main__":
    N = 1
    args = len(sys.argv)
    if args > 1:
        N = int(sys.argv[1])
    print("This is a test of Goldbach's Conjecture that for all even integers")
    print("greater than 2 there are two primes that add up to that even number.\n")
    while (N < 3 or N%2):
        N = int(input("Please enter an even number > 3 to check with Goldbach's Conjecture> "))
    cProfile.run('goldbach(N)')