Coding Lemoine's Conjecture in Python

Question

I am trying to code a function that when given N, will return a pair of numbers that comply with Lemoine's Conjecture (every odd number greater than 5 can be expressed as the sum of a prime and the double of a prime). I have created this code based on a previous function relating to Goldbach's Conjecture (this function works fine) and have used a different function to generate a list of prime numbers up to N, however my new code isn't giving me correct results and I can't see why- any ideas? Thanks

def eratosthenes2(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 lemoine(N):
    x, y = 0, 0
    result = 0
    if N % 2:
        prime = eratosthenes2(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 = 2*x + y
                    if result == N: 
                        break 
    return x, y 

Answer 1

Here. This works fairly quickly. The problem with getting all the pairs is that for large numbers like 2**45+1 it takes a while to collect thousands and thousands of pairs.

'''
Lemoine's Conjecture (every odd number greater than 5 can be expressed as the sum of a prime and the double of a prime). 
'''
from gmpy2 import is_prime 
# or use a home grown one
#def is_prime(n):
#  if n<=3: return n>1
#  if n%6 not in [1,5]: 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 gen_primes():
    """Generator that yields the next prime number."""
    yield 2
    yield 3
    num = 5
    while True:
        if is_prime(num):
            yield num
        if is_prime(num+2):
            yield num+2
        num += 6

def lemoine(n,once=True):
   pairs = []
   if n < 7 or not n%2: return pairs
   primes = gen_primes()
   p = next(primes)
   while n-p*2 > 2:
      if is_prime(n-(p*2)):
          pairs+=[(n-(p*2),p)]
          if once: return pairs
      p=next(primes)
   return pairs
   
start = 2**64+1
end = start+200
for i in range(start,end+1):
   pairs = lemoine(i)
   #pairs = lemoine(i,once=False) # get all pairs 
   for a,b in pairs: 
      print(f'{i} = {a} + 2 * {b}  or {i} = {a} + {2*b}') 

Output:

$ python lemoine.py 
18446744073709551617 = 18446744073709551163 + 2 * 227  or 18446744073709551617 = 18446744073709551163 + 454
18446744073709551619 = 18446744073709551557 + 2 * 31  or 18446744073709551619 = 18446744073709551557 + 62
18446744073709551621 = 18446744073709551427 + 2 * 97  or 18446744073709551621 = 18446744073709551427 + 194
18446744073709551623 = 18446744073709550341 + 2 * 641  or 18446744073709551623 = 18446744073709550341 + 1282
18446744073709551625 = 18446744073709551263 + 2 * 181  or 18446744073709551625 = 18446744073709551263 + 362
18446744073709551627 = 18446744073709551533 + 2 * 47  or 18446744073709551627 = 18446744073709551533 + 94
18446744073709551629 = 18446744073709551427 + 2 * 101  or 18446744073709551629 = 18446744073709551427 + 202
18446744073709551631 = 18446744073709551557 + 2 * 37  or 18446744073709551631 = 18446744073709551557 + 74
18446744073709551633 = 18446744073709551629 + 2 * 2  or 18446744073709551633 = 18446744073709551629 + 4
18446744073709551635 = 18446744073709551629 + 2 * 3  or 18446744073709551635 = 18446744073709551629 + 6
18446744073709551637 = 18446744073709551359 + 2 * 139  or 18446744073709551637 = 18446744073709551359 + 278
18446744073709551639 = 18446744073709551629 + 2 * 5  or 18446744073709551639 = 18446744073709551629 + 10
18446744073709551641 = 18446744073709551427 + 2 * 107  or 18446744073709551641 = 18446744073709551427 + 214

Answer 2

If you have to use the sieve method, then first:

Change your sieve to use n instead of n + 1:

primes = list (range(2, n))

Next, change your lemoine function to have:

if result == N:                             
    return x, y 

instead of the current way, where you break when result == N. The way you have it now, you're exiting the function after x has been incremented one more time, leading to incorrect results. (e.g. 2 instead of 3 in the below n = 47 example.

Here's a working implementation to compare with:

def isPrime (n):
  if n < 2:
    return False

  for i in range(2, (int(n ** (1/2)) + 1)):
    if n % i == 0:
      return False
  return True

def lemoine(n):
  pairs = {}

  # n = p + (2 * q)
  for q in range(1, int(n / 2)):
    p = n - 2 * q

    # Are p and q prime?
    if isPrime(p) and isPrime(q):
      pairs[p] = q

  return pairs

n = 47
pairs = lemoine(n)
for key in pairs:
  print('{} is {} + 2 * {}'.format(n, key, pairs[key]))

This gives the output:

47 is 43 + 2 * 2
47 is 41 + 2 * 3
47 is 37 + 2 * 5
47 is 13 + 2 * 17

This example is from the Wikipedia page.