Goldbach Conjecture in Python

Question

I have attempted to write a code that returns a single pair that satisfies the Goldbach Conjecture for a given N. The conjecture states that every even number greater than 4 can be expressed as the sum of two prime numbers. The function returns a pair that is just slightly off, for example, goldbach(34) returns (5, 31) rather than the correct answer (3, 31). Similarly, goldbach(38) returns (11, 31). Any ideas where I am going wrong here? I understand that this code is not massively efficient, however this is the way I have been asked to code for my assignment.

def eratosthenes(n):
    primes = list (range(2, n+1))
    for i in primes:
        j=2
        while i*j<= primes[-1]:
            if i*j in primes:
                primes.remove(i*j)
            j=j+1
    return primes

def odd_primes(N):
    oddprimes = eratosthenes(N)
    oddprimes.remove(2)
    return(oddprimes)

def goldbach(N):
    x, y = 0, 0
    result = 0
    if N % 2 == 0:
        prime = odd_primes(N)
        while result != N:
            for i in range(len(prime)):
                x = prime[i]
                if result == N: break
                for j in range(len(prime)):
                    y = prime[j]
                    result = x + y
                    if result == N: break 
    return x, y 

Answer 1

You are assigning x before breaking the loop once your condition is met. Simply invert your break lines in the first for loop:

def goldbach(N):
    x, y = 0, 0
    result = 0
    if N % 2 == 0:
        prime = odd_primes(N)
        while result != N:
            for i in range(len(prime)):
                if result == N: break  # this line first
                x = prime[i]   # this line after
                for j in range(len(prime)):
                    y = prime[j]
                    result = x + y
                    if result == N: break 
    return x, y 

Answer 2

def eratosthenes(n):
    primes = list (range(2, n+1))
    for i in primes:
        j=2
        while i*j<= primes[-1]:
            if i*j in primes:
                primes.remove(i*j)
            j=j+1
    return primes

def odd_primes(N):
    oddprimes = eratosthenes(N)
    oddprimes.remove(2)
    return(oddprimes)

def goldbach(N):
    x, y = 0, 0
    result = 0
    if N % 2 == 0:
        prime = odd_primes(N)
        while result != N:
            for i in range(len(prime)):
                if result == N: break 
                x = prime[i]
                for j in range(len(prime)):
                    y = prime[j]
                    result = x + y
                    if result == N: break 
    return x, y 

is the correct version. When you have found a pair, you set x to the next prime.

Answer 3

def is_prime(a):
    for div in range(2, a//2 + 1):
        if a % div == 0:
            return False
    return True


def goldbach(n):
    summ = 0
    for p1 in range(2, n):
        if is_prime(p1) and is_prime(n-p1):
            p2 = n-p1
            return str(n)+"="+str(p1)+"+"+str(p2)

Answer 4

    def isPrime(n):
        for i in range(2,n):
            if n%i == 0:
                return 0
        return 1

    no = int(input("Enter Number: "))

    for i in range(3,no):
        if isPrime(i) == 1:
            for l in range(i,no):
                if isPrime(l) == 1:
                    if no == (i+l):
                        print(i,"+",l,"=",no)

Answer 5

Goldbach Conjecture is observed only in even numbers greater than 2, therefore I improve on Punnie's answer to catch odd numbers and any number less done 2.

def isPrime(n):
    for i in range(2,n):
        if n%i == 0:
            return 0
    return 1

def goldbach(no):
  if no%2 !=0 :
    return "Error {} is not an even number ".format(no) 
  elif no <= 2:
    return "Error {} is not greater than 2, Goldbach Conjecture is observed only in even numbers greater than 2".format(no)

  else:
      for i in range(3,no):
          if isPrime(i) == 1:
              for l in range(i,no):
                  if isPrime(l) == 1:
                      if no == (i+l):
                          print(i,"+",l,"=",no)

no = int(input("Enter Even Number greater than 2: "))
goldbach(no)            ```

Answer 6

It is possible to use a sieve that separates the primes 1 (mod 6) with the ones -1 (mod 6) in this way it is possible to speed up using numpy without using the for loop:

import numpy as np

def sieve(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]
    return P5m6,P1m6


def goldbach(n):
    P5m6,P1m6=sieve(n)
    nmod6=n%6
    if nmod6==0:
        r=(6*np.nonzero(P5m6[1:n//6]&P1m6[n//6-1:0:-1])[0]+5)[0]
    elif nmod6==2:
        if P5m6[n//6]:
            r=3
        else:
            r=(6*np.nonzero(P1m6[1:(n-6)//12+(n//6-1)%2+1]&P1m6[n//6-1:(n-6)//12:-1])[0]+7)[0]
    elif nmod6==4:
        if P1m6[n//6]:
            r=3
        else:
            r=(6*np.nonzero(P5m6[1:n//12+(n//6)%2+1]&P5m6[n//6:n//12:-1])[0]+5)[0]
    else:
        r=0
    return r,n-r

Answer 7

The Goldbach decomposition of an even number into a pair of primes numbers is not (always) unique so you need a method to select one possible solution (as required from the question). A natural possibility would be using min/max function.

I did not use the function odd_primes (just a slice) and reimplemented goldbach function using a combinatorics approach.

import itertools as it

def eratosthenes(n):
    primes = list (range(2, n+1))
    for p in primes:
        j = 2
        while p*j <= primes[-1]:
            if p*j in primes:
                primes.remove(p*j)
            j += 1
    return primes

def goldbach(n):
    if n % 2 != 0 or n <= 2: raise Exception('Conjecture assumptions not satisfied.')
    ps = eratosthenes(n)[1:]
    # trivial cases to overcome limitations of the combinations
    if n == 4: return [(2, 2)]
    if n == 6: return [(3, 3)]

    return max(it.combinations(ps, 2), key=lambda p: p[1] if sum(p) == n else 0)

for n in range(6, 40, 2):
    print(n, goldbach(n))

Output

8 (3, 5)
10 (3, 7)
12 (5, 7)
14 (3, 11)
16 (3, 13)
18 (5, 13)
20 (3, 17)
22 (3, 19)
24 (5, 19)
26 (3, 23)
28 (5, 23)
30 (7, 23)
32 (3, 29)
34 (3, 31)
36 (5, 31)
38 (7, 31)

Answer 8

You want to get, for any even number N, a pair of prime number whose sum is N.

For this, we will simply go through each number between 2 and N. Let's call this number i. We will then check if i and N-i are both prime.

If that is the case, then we have a pair of prime number so that their sum is N.

Indeed, i + (N-i) = N.

def is_prime(number):
    for div in range(2, number//2 + 1):
        if number % div == 0:
            return False
    return True


def compute_goldbach_pair(number):
    for p1 in range(2, number):
        p2 = number - p1
        if is_prime(p1) and is_prime(p2):
            return [p1, p2]

In this code, we return as soon as we have found a pair. But multiple pairs can exists, in this case, you'll have to modify the code a bit to collect each existing pairs.

Answer 9

Don't check for prime unless forced.

This takes milliseconds for trillions. Using gmpy2 or home grown IsPrime. Just don't call for a Prime check unless you have to. For very large numbers you only call about 50 times.

# Goldbach's Conjecture tester. 
from gmpy2 import is_prime
import sys
import cProfile
# or use home grown IsPrime
def IsPrime(n): 
        if (n == 2 or n == 3): 
            return True 
        if (n <= 1 or n % 2 == 0 or n % 3 == 0): 
            return False 
        for i  in range(5, int(n**.5)+1,  6): 
            if (n % i == 0 or n % (i + 2) == 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 odds 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)')

Runs very fast with GMPY2.is_prime():

   python goldbach.py 4444444444444444
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.


107 + 4,444,444,444,444,337 = 4,444,444,444,444,444

         55 function calls in 0.001 seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    0.001    0.001 <string>:1(<module>)
        1    0.000    0.000    0.001    0.001 goldbach.py:16(goldbach)
        1    0.000    0.000    0.001    0.001 {built-in method builtins.exec}
        1    0.001    0.001    0.001    0.001 {built-in method builtins.print}
       50    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 Home grown IsPrime():

python goldbach.py 4444444444444444
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.


107 + 4,444,444,444,444,337 = 4,444,444,444,444,444

         55 function calls in 3.957 seconds

   Ordered by: standard name

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

10 trillion using all python(No GMPY2)

python goldbach.py 10000000000000
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.


29 + 9,999,999,999,971 = 10,000,000,000,000

         20 function calls in 0.160 seconds

   Ordered by: standard name

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