jemima
jemima

Reputation: 29

Weak Goldbach Conjecture in python

I've tried to write a code for the weak Goldbach Conjecture, which states that every odd number greater than 5 can be expressed as the sum of three prime numbers. However, the code only returns (0, 0, 0). I only need one triple that works rather than a list of triples. Any ideas where I'm going wrong? Also I know that the code isn't the most efficient as it is, especially using the eratosthenes function to generate my primes, but this is the form that I've been asked to code in.

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 weak_goldbach(N):
    x, y, z = 0, 0, 0
    result = 0
    if not N % 2:
        prime = eratosthenes(N)
        while result != N:
            for i in range(len(prime)):
                x = prime[i]
                if result == N: 
                    break
                for j in range(i, len(prime)):
                    y = prime[j]
                    result = x + y
                    if result == N: 
                        break 
                    for k in range (j, len(prime)):
                        z = prime[k]
                        result = x + y + z
                        if result == N:break
    return x, y, z

Upvotes: 1

Views: 880

Answers (2)

Andy Richter
Andy Richter

Reputation: 189

2^100+3 in less than a second


2^1000 takes about ten seconds

You can always assume the first prime number is 3. Then subtract three from the number and then solve for Goldbach's Strong Conjecture with the even number you have left. This code works quickly for fairly large numbers I like to test with 2**100+3

# Goldbach's Conjecture tester using Miller-Rabin primality test. 
import sys,secrets 
#from gmpy2 import is_prime as IsPrime
# or use home grown IsPrime

def IsPrime( n ):
    from secrets import randbelow
    k = 6
    if (n <= 3):
        return n > 1
    if n%6 not in [1,5]: return False
    d = n - 1;
    while (d % 2 == 0):
        d //= 2;
    for i in range(k):
        x = pow(2 + randbelow(n - 4), d, n)
        if (x == 1 or x == n - 1): continue
        while (d != n - 1):
           x = (x * x) % n
           d *= 2;
           if (x == 1): return False
           if (x == n - 1):
              break
        else: return False
    return True        

def goldbach(N:int):
    number = N
    gold = []
    if number%2 and number > 5: 
        gold.extend([3])
        number -= 3 # Change to even 
    if number%2 or number < 4: 
        return [0,0,0]
    # Now use Goldbach's Strong Conjecture on this even number > 3 
    if number == 4:
       gold.extend([2,2])
       return gold
    elif IsPrime(number-3):
        gold.extend([3,number-3]) 
        return gold
    else:
       for p,q in ((i,i+2) for i in range(5,number//2+1,6)): 
            if IsPrime(number-p) and IsPrime(p): 
                gold.extend([p,number-p ])
                return gold
            if IsPrime(number-q) and IsPrime(q):
                gold.extend([q,number-q ])
                return gold
    raise Exception("Looks like %d is the exception to the weak Goldbach conjecture!" % N)
    return

if __name__=="__main__":
    args = len(sys.argv)
    #g = 121005022304007026
    g = 2**100+3
    if args > 1:
        g = int(sys.argv[1])
    j = "" if (g%2==0) else "Weak " 
    print(f"Test Goldbach's {j}Conjecture.");  
    print(f"{g} = ",end =''); 
    print(*goldbach(g),sep=" + ") 

Output:

 $ time python goldbach_miller.py 1267650600228229401496703205379
Test Goldbach's Weak Conjecture.
1267650600228229401496703205379 = 3 + 479 + 1267650600228229401496703204897

real    0m0.036s
user    0m0.027s
sys     0m0.009s

 $ time python goldbach_miller.py 12345678901234567
Test Goldbach's Weak Conjecture.
12345678901234567 = 3 + 191 + 12345678901234373

real    0m0.034s
user    0m0.029s
sys     0m0.005s

$ time python goldbach_miller.py 123456789012345678
Test Goldbach's Conjecture.
123456789012345678 = 7 + 123456789012345671

real    0m0.034s
user    0m0.028s
sys     0m0.006s

2^1000:

$ time python goldbach_miller.py 
Test Goldbach's Conjecture.
10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069376 = 16607 + 10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668052769

real    0m10.849s
user    0m10.831s
sys     0m0.013s

Using a Sieve. Slower and more limited but:

This code uses a bitarray for the sieve. So the sieve can be larger. it is also faster than a sieve of integers. 2^30+3 takes about 5 seconds. So slow and, even with a bitarray(), it cannot do 2**50+3 at all.

from math import isqrt
from bitarray import bitarray
import sys

def sieve_of_bits(n):
    """Create a sieve of Eratosthenes up to the given size."""
    p = bitarray(n)
    p.setall(True)
    p[0:2] = False   # Clear zero and one
    p[4::2] = False  # Clear multiples of 2
    p[9::3] = False  # Clear multiples of 3
    for k6m,k6p in ((k6,k6+2) for k6 in range(5,isqrt(n)+1,6)): # Process only numbers of the form 6k-1 and 6k+1
        if p[k6m]:  # If 6k-1 is prime
            p[k6m*k6m::2*k6m] = False  # Clear multiples of 6k-1
        if p[k6p]:  # If 6k+1 is prime
            p[k6p*k6p::2*k6p] = False  # Clear multiples of 6k+1
    return p
    # for array of ints: return p.search(bitarray('1')) 

    
def goldbach(N:int):
    number = N
    prime = sieve_of_bits(N)
    gold = []
    if number%2 and number > 5: 
        gold.extend([3])
        number -= 3 # Change to even 
    if number%2 or number < 4: 
        return [0,0,0]
    # Now use Goldbach's Strong Conjecture on this even number > 3 
    if number == 4:
       gold.extend([2,2])
       return gold
    elif prime[number-3]:
        gold.extend([3,number-3]) 
        return gold
    else:
       for p,q in ((i,i+2) for i in range(5,number//2+1,6)): 
            if prime[number-p] and prime[p]: 
                gold.extend([p,number-p ])
                return gold
            if prime[number-q] and prime[q]:
                gold.extend([q,number-q ])
                return gold
    raise Exception("Looks like %d is the exception to the weak Goldbach conjecture!" % N)
    return

if __name__=="__main__":
    args = len(sys.argv)
    #g = 121005022304007026
    g = 2**30+3
    if args > 1:
        g = int(sys.argv[1])
    j = "" if (g%2==0) else "Weak " 
    print(f"Test Goldbach's {j}Conjecture.");  
    print(f"{g} = ",end =''); 
    print(*goldbach(g),sep=" + ") 

Output:

$ time python goldbach_sieve.py 
Test Goldbach's Weak Conjecture.
1073741827 = 3 + 41 + 1073741783

real    0m4.627s
$ time python goldbach_sieve.py 1234567890
Test Goldbach's Conjecture.
1234567890 = 79 + 1234567811

real    0m5.437s

Upvotes: 0

Silas Coker
Silas Coker

Reputation: 500

Bugs

There are a few problems with your code, but the first, and the reason it fails, is that not N % 2 always evaluates to false for odd numbers, skipping your loop and returning the initial values you set x,y and z.

There are also logical problems with it; your code breaks in the most inner loop when x+y+z == N, then the outer loops break when result is correctly set, but only after changing x or y! This means that even if you fixed the first problem, your code would always return the wrong result.

Inefficiencies

Firstly, you don't need the complicated break logic at all! You can just return x, y, z when you first find it sums to N.

Secondly, you don't need the result=x+y code in the middle loop, since it has no relevance in the weak Goldbach conjecture and is never true anyway.

Thirdly, the outer while loop is completely useless. It does nothing at all, except create an infinite loop if the inner loops didn't find a result for some reason.

Other issues

It's best to reduce nesting; so, the condition making sure it is an odd number greater than 5 should be negative, and should return an exception if it doesn't hold. That way the body of the function isn't nested in a conditional, which makes it a bit more readable.

Corrected code that works

def weak_goldbach(N):
     x, y, z = 0, 0, 0
     result = 0
     if not N % 2 or N < 7:
         raise Exception("Bad input - must be odd number greater than 5.")
     prime = eratosthenes(N)
     for i in range(len(prime)):
         x = prime[i]
         for j in range(i, len(prime)):
             y = prime[j]
             for k in range (j, len(prime)):
                 z = prime[k]
                 if x+y+z == N:
                     return x, y, z
     raise Exception("Looks like %d is the exception to the weak Goldbach conjecture!" % N)

Upvotes: 3

Related Questions