Goldbach's Conjecture

I am doing an exercise on finding the nearest pair of prime numbers p₁ and p₂ such that:

• p₁ ≤ p₂
• p₁ + p₂ = n (for 4 ≤ n ≤ 10⁷)
• p₁ and p₂ are primes
• p₂ - p₁ is minimal

I think it should work perfectly except that it got a TLE. Can someone point out how can I avoid it?

import math
import random

def _try_composite(a, d, n, s):
if pow(a, d, n) == 1:
return False
for i in range(s):
if pow(a, 2 ** i * d, n) == n - 1:
return False
return True  # n  is definitely composite

def is_prime(n, _precision_for_huge_n=16):
if n in _known_primes:
return True
if n in (0, 1) or any((n % p) == 0 for p in _known_primes):
return False
d, s = n - 1, 0
while not d % 2:
d, s = d >> 1, s + 1
# Returns exact according to http://primes.utm.edu/prove/prove2_3.html
if n < 1373653:
return not any(_try_composite(a, d, n, s) for a in (2, 3))
if n < 25326001:
return not any(_try_composite(a, d, n, s) for a in (2, 3, 5))

_known_primes = [2, 3]
_known_primes += [x for x in range(5, 1000, 2) if is_prime(x)]

def main():
n = int(input())
if n > 4:
for smaller in range(n // 2, -1, -1):
if n - smaller >= smaller:
if is_prime(n - smaller) and is_prime(smaller):
print (smaller, n - smaller)
flag = True
break
else:
print ('2 2')

main()

• Since I need to minimize p_2 - p_1, it should be from n/2 to 3 and thus doesn't change? Mar 13, 2018 at 2:29
• Also using a sieve is a good idea... except sieving from 1 to 10^7 takes around 19.3 second of my code. Mar 13, 2018 at 2:29
• Mar 13, 2018 at 9:57
• @TobySpeight Sorry, I did not copy it. I typed it out, so it is a typo. Mar 15, 2018 at 4:43
• Typing it yourself is the form of copying that I was thinking of - thanks for correcting the description. It's good to know it's not a bug in the code! Mar 15, 2018 at 8:42

My solution consists of a single change:

• Changed the is_prime function into a prime sieve (from this answer) + lookup

My final code (Also handled odd numbers just for fun, definitely not because I misread the question):

def main():
import itertools
import math
izip = itertools.zip_longest
chain = itertools.chain.from_iterable
compress = itertools.compress
def rwh_primes2_python3(n):
""" Input n>=6, Returns a list of primes, 2 <= p < n """
zero = bytearray([False])
size = n//3 + (n % 6 == 2)
sieve = bytearray([True]) * size
sieve[0] = False
for i in range(int(n**0.5)//3+1):
if sieve[i]:
k=3*i+1|1
start = (k*k+4*k-2*k*(i&1))//3
sieve[(k*k)//3::2*k]=zero*((size - (k*k)//3 - 1) // (2 * k) + 1)
sieve[  start ::2*k]=zero*((size -   start  - 1) // (2 * k) + 1)
ans = [2,3]
poss = chain(izip(*[range(i, n, 6) for i in (1,5)]))
ans.extend(compress(poss, sieve))
return ans

string2 = "Impossible"
n = int(input())
sieve = [t for t in rwh_primes2_python3(n) if t <= math.floor((n // 2) / 2) * 2 + 1][::-1]
another = [t for t in rwh_primes2_python3(n) if t >= math.floor((n // 2) / 2) * 2 + 1]

if n > 5 and n % 2 == 0:
for smaller in sieve:
if n - smaller in another:
print (smaller, n - smaller)
break
elif n % 2 == 1 and n != 5:
if n - 2 in another or n-2 in sieve: print (2, n-2)
else: print ('Impossible')
elif n == 4:
print ('2 2')
else:
print ('2 3')

main()

• One advantage of this approach is that if you want to test multiple n, you just sieve for the largest one, and that cost is amortized over all the inputs. I should point out that you invoke rwh_primes2_python3(n) twice - I think you can call it just once, and then split the result into upper and lower halves (actually just one half, for the for smaller in another loop - assuming primes get sparser, that's a better choice - and we can keep sieve complete for the if n-smaller in sieve test). Mar 14, 2018 at 8:39

I opted to not create a sieve or a list of primes. Instead I check each prime only when forced to check. As a result, for numbers larger than 2000 many fewer primes are found. This increases the speed exponentially the larger the number becomes. 18 and 20 digit numbers require fewer than 100 prime checks usually. This code will find 18 digit even numbers pair of primes in about 20 seconds. 9 digits takes less than a second.

# 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
else:
mid = int(number//2) + (1 - (number//2%2) )
for p in range(mid, 0, -2): # just odds
if IsPrime(p) and IsPrime(number-p):
print(f"\n{p:,} + {number-p:,} = {number:,}\n")
return
elif IsPrime(p+2) and IsPrime(number-(p+2)):
print(f"\n{p+2:,} + {number-(p+2):,} = {number:,}\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)')


Output:

python3 goldbach_near.py
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.

Please enter an even number > 3 to check with Goldbach's Conjecture> 12345678

6,172,799 + 6,172,879 = 12,345,678

$python3 goldbach_near.py 12348 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. 6,151 + 6,197 = 12,348$ python3 goldbach_near.py 44000000000000000
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.

21,999,999,999,999,979 + 22,000,000,000,000,021 = 44,000,000,000,000,000

30 function calls in 14.669 seconds