Goldbach's Conjecture - Find the number of ways an even number can be written as the sum of two primes

Question

I want to know the number of ways that a given positive even number can be written as the sum of two primes.

At the moment, I have this code:

n = int(input(">"))
def genprimes(n):#generate primes from 2 to n
    primes = []
    for i in range(2,n):
        prime = True
        for a in range(2,i):
            if i % a == 0:
                prime = False
        if prime == True:
            primes.append(i)
    return primes

pairs = []
primes = genprimes(n)

for prime1 in primes:
    for prime2 in primes:
        if prime1 + prime2 == n and [prime1,prime2] not in pairs and [prime2,prime1] not in pairs:
            pairs.append([prime1,prime2])
print(len(pairs))

It works, but it becomes a little slow when n goes above 10,000 ish. Can anyone think of a more efficient way to find this value, so that it yields quick results for high values?

Answer 1

You have two slow algorithms.

To get your list of primes, you separately check each number for prime-hood. The fastest way get a list of primes is to use a pre-built list of primes--that is what I do for such problems. The list of primes up to 2**16 (65,536) takes only 6542 integers and I store that list in a module. You may prefer to use the Sieve of Eratosthenes which is generally recognized as the fastest way to get such a list from scratch.

You then try each pair of primes to see if they add to your target number. It would be much faster to, for each prime p, see if n-p is also prime. You could check n-p with a binary search, as @Shubham suggests, but it would probably be faster to have two indices, one regularly going up and yielding p and the other going down as needed checking for n-p. It might be quickest to copy your prime list to a set and use that to check if n-p is in the list. These last two techniques are of order n, but for numbers only up to 10,000 the overhead may affect which technique is best. And yes, for such problems 10,000 is not very large. Finally, if memory is not a concern, @BoulderKeith points out that you could leave the prime list as a list of Booleans (True/False) and very quickly check if n-p is prime--this is the fastest, if most memory-consuming, technique of all.

Answer 2

This can be done efficiently by optimizing both part of your program.

First, we can optimize the genprimes() by generating primes using Sieve of Eratosthenes which can generate prime upto n in O(nlog(log(n)).

After that, once we have the list of primes, we can do the following:

• Iterate through the prime list and select a prime, say p1.
• Using Binary Search check whether n-p1 also exist in the prime list.

If p = number of primes till n

then the complexity of this part will be O(plog(p)).

As @PaulHankin suggested, from @SrujanKumarGulla's answer, we can build up a O(p) solution once we have the prime list.

Implementation of second part (the first part is standard Sieve of Eratosthenes):

# prime list is assumed to be sorted. It will be if it's generated
# using Sieve of Eratosthenes

def solutions(n, prime_list):
    low = 0, high = len(prime_list)-1
    sol = []
    while low < high:
        temp = prime_list[low] + prime_list[high]
        if temp == n:
            sol.append((prime_list[low], prime_list[high], ))
        elif temp < n:
            low += 1
        else:
            high -= 1

    return len(sol)

Answer 3

If memory is not an issue, you can optimize speed even more by:

1) Building an array from 2..n with array[i]=true if i is prime. (Sieve of Eratosthenes or something fancier can be used.)

2) For all 2 <= i < n/2 see if both array[i] and array[n-i] are both true. If so, then i and n-i are a pair. (Setting the upper limit of i to n/2 means you don't have to remove duplicates either.)

The entire algorithm is O(n log n). (O(n log n) to build the primes; O(n) to find pairs.)

I did this in C# rather than Python, but was able to find all pairs with n = 100,000,000 in 10 seconds.

Answer 4

    def prim(n):
        primes = []
        for num in range(1, n + 1):     # current number check
            limit = int(num ** 0.5) + 1 # search limit
            for div in range(2, limit): # searching for divider
                if (num % div) == 0:    # divider found
                    break               # exit the inner loop
            else:                       # out of the inner loop
                primes.append(num)      # no divider found hence prime
        return(primes)

    num = int(input('Input an even number: '))
    primes = prim(num)                  # create list of primes

    for i in range(num//2 + 1):
        if num % 2 != 0:
            print(num, 'not even'); break
        if i in primes and (num - i) in primes:
            print(i, '+', num - i)
"""
Examples:
=========
Input an even number: 18
1 + 17
5 + 13
7 + 11

Input an even number: 88
5 + 83
17 + 71
29 + 59
41 + 47

Input an even number: 17
17 not even
"""

Answer 5

This is the actual solution from Rory's suggestions, which uses the Sieve of Eratosthenes to efficiently generate the primes and then sees if n-prime is also prime.

The count of solutions for which this is true then needed to be halved and rounded up to prevent duplicates.

n = int(input(">"))
def genprimes(n):
    limitn = n+1
    not_prime = set()
    primes = []

    for i in range(2, limitn):
        if i in not_prime:
            continue
        for f in range(i*2, limitn, i):
            not_prime.add(f)
        primes.append(i)
    return primes

nos=0
primes = genprimes(n)

for prime in primes:
    if n-prime in primes:
        nos += 1
print(str(int((nos/2)+0.5)))

Answer 6

With this version of the Sieve of Eratosthenes it is possible to quickly find the number of ways in which a even number can be written as the sum of two primes:

import numpy as np
def goldbach(n):
    P5m6 =np.ones((n//6+1), dtype=bool)
    P1m6 =np.ones((n//6+1), dtype=bool)
    for i in range(1,int((n**0.5+1)/6)+1):
        if P5m6[i]:
            P5m6[6*i*i::6*i-1]=[False]
            P1m6[6*i*i-2*i::6*i-1]=[False]
        if P1m6[i]:
            P5m6[6*i*i::6*i+1]=[False]
            P1m6[6*i*i+2*i::6*i+1]=[False]

    nmod6=n%6
    if nmod6==0:
        return np.count_nonzero(P5m6[1:n//6]&P1m6[n//6-1:0:-1])
    elif nmod6==2:
        return np.count_nonzero(P1m6[1:(n-6)//12+(n//6-1)%2+1]&P1m6[n//6-1:(n-6)//12:-1])+np.count_nonzero(P5m6[n//6])
    elif nmod6==4:
        return np.count_nonzero(P5m6[1:n//12+(n//6)%2+1]&P5m6[n//6:n//12:-1])+np.count_nonzero(P1m6[n//6])
    else:
        return -1
print(goldbach(100000000))