Lemoine's Conjecture Verified to 10^10

The mathematician who first thought of this idea: Emile Lemoine

Lemoine's conjecture, sometimes referred to as Levy's conjecture, states that all odd numbers can be broken into the sum of an odd prime and an even semi-prime.

Here are three examples of "lemoine partitions", i.e. where the odd numbers are broken down into their constituent elements based on the definition above.
$$7 = 3 + 2 * 2$$
$$9 = 3 + 2 * 3$$
$$11 = 5 + 2 * 3$$
In 1999, Dann Corbit was able to verify the conjecture to 1,000,000,000 numbers utilizing the following code in C (https://groups.google.com/forum/?hl=en#!msg/sci.math/HoCpH8gPDHM/gtU-fuSWnFMJ):

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#define TEST( f, x ) ( *( f+( x )/16 )&( 1<<( ( ( x )%16L )/2 ) ) )
#define SET( f, x ) *( f+( x )/16 )|=1<<( ( ( x )%16L )/2 )

unsigned        ip(unsigned long j, unsigned char *d)
{
    if (j == 1)
        return 0;
    if ((j % 2) == 0) {
        if (j == 2)
            return 1;
        else
            return 0;
    } else
        return (!(*(d + (j) / 16) & (1 << (((j) % 16L) / 2))) && j != 1);
}

int             main(int argc, char *argv[])
{
    unsigned char  *feld = NULL,
    *zzz = NULL;
    unsigned long   teste = 1,
    max,
    mom,
    hits = 1,
    alloc;
    time_t          begin;
    unsigned char   quiet = 0;
    if (argc > 1) {
        max = atol(argv[1]) + 10000;
        if (argc > 2)
            quiet = 1;
    } else
        max = 1000000000L;
    
    zzz = feld = malloc(alloc = (((max -= 10000L) >> 4) + 1L));
    if (feld) {
        char            found;
        memset(zzz, 0, alloc);
        printf("Searching prime numbers to : %lu\n", max);
        begin = time(NULL);
        while ((teste += 2) < max)
            if (!TEST(feld, teste)) {
                ++hits;
                for (mom = 3L * teste; mom < max; mom += teste << 1)
                    SET(feld, mom);
            }
        printf(" %ld prime numbers found in %.2f secs.\n\n",
               hits, difftime(time(NULL), begin));
        {
            long            o,
            p,
            q,
            j;
            for (o = 7; o <= max; o += 2) {
                found = 0;
                for (j = 2; j < o;) {
                    p = j;
                    q = o - 2 * p;
                    if (ip(p, feld) && ip(q, feld)) {
                        if (!quiet)
                            printf("%lu = 2*%lu + %lu\n", o, p, q);
                        found = 1;
                        break;
                    }
                    if (j == 2)
                        j++;
                    else
                        j += 2;
                }
                if (!found)
                    printf("*** Problem!! %lu\n", o);
            }
            
            free(feld);
            printf(" Finished calculations in %.2f secs.\n\n",
                   difftime(time(NULL), begin));
        }
    } else {
        puts("Memory allocation failure.");
        exit(EXIT_FAILURE);
    }
    return 0;
}

A new verifier was implemented in Python to try and either match/exceed this number.

"""
Description: Traditional Implementation of Lemoine's Conjecture
"""

import numpy as np
import time

def primesfrom2to(n):
    # https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Input n>=6, Returns a array of primes, 2 <= p < n """
    sieve = np.ones(n/3 + (n%6==2), dtype=np.bool)
    sieve[0] = False
    for i in xrange(int(n**0.5)/3+1):
        if sieve[i]:
            k=3*i+1|1
            sieve[      ((k*k)/3)      ::2*k] = False
            sieve[(k*k+4*k-2*k*(i&1))/3::2*k] = False
    return np.r_[2,3,((3*np.nonzero(sieve)[0]+1)|1)]

def main():
    start_time = time.time()
    n = 7
    x = 0
    l = list(primesfrom2to(1000000000))
    sl = set(l)
    print("Done")
    print("--- %s seconds ---" % (time.time() - start_time))            

    while (n < 1000000000):
        if(n-(2*l[x]) in sl):
            n = n+2
            x = 0
        x = x + 1
    print("Done 2")  
    print("--- %s seconds ---" % (time.time() - start_time))            
    
main()

This version relied on a pre-generated list of prime numbers in both list and set form to generate odd numbers sequentially. It rearranges the formal definition of Odd # = 2q + p to Odd # - 2q = p to take advantage of python's native set data type and enhanced searching capabilities. Replacing these with just lists, numpy arrays, or iterators was considered but not successfully implemented in a time-conserving manner.

You can track the progress of the generator by inserting this code snippet within the if statement within the while loop within the main method.

if n % 1000 == 1:
                print("\r{}/{}".format(n, 100000))

Here is a table with the speed results from this implementation. The program reached 2*10^9 before stressing the hardware due to lack of available memory.

Odd #’s up to ___ checked	Time
10,000	--- 0.0129 seconds ---
100,000	--- 0.129 seconds ---
1,000,000	--- 1 seconds ---
10,000,000	--- 11 seconds ---
100,000,000	--- 138 seconds ---
1,000,000,000	--- 1593 seconds ---
2,000,000,000	--- 3,487 seconds ---

Another approach would be to generate the results from all of the combinations of p + 2q and match those up to a list of the odd numbers greater than or equal to seven. This was implemented with the following program.

"""
Description: New strategy to tackle verifying Lemoine's conjecture
"""

import numpy as np
import time
from operator import add

def primesfrom2to(n):
    # https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Input n>=6, Returns a array of primes, 2 <= p < n """
    sieve = np.ones(n/3 + (n%6==2), dtype=np.bool)
    sieve[0] = False
    for i in xrange(int(n**0.5)/3+1):
        if sieve[i]:
            k=3*i+1|1
            sieve[      ((k*k)/3)      ::2*k] = False
            sieve[(k*k+4*k-2*k*(i&1))/3::2*k] = False
    return np.r_[2,3,((3*np.nonzero(sieve)[0]+1)|1)]

def cartesian_add(arr1, arr2):
    arr1_e = np.repeat(np.expand_dims(arr1, 1), arr2.size, axis=1)
    arr2_e = np.repeat(np.expand_dims(arr2,  0), arr1.size, axis=0)
    return arr1_e + arr2_e

def main():
    start_time = time.time()
    list1 = primesfrom2to(1000000)
    list2 = 2 * list1
    list4 = np.arange(7,3000000,2)
    list4 = np.setdiff1d(list4,np.unique(cartesian_add(list1, list2)))
    print(list4)    
    print("--- %s seconds ---" % (time.time() - start_time))
  
main()

While the original versions of this implementation utilized nested for loops (later changed to for loop containing a map(add,primes,semiprimes) function with rolling numpy primes loop), Ernie Parke contributed the cartesian_add method which shaved off several seconds (it runs about 4x faster without the use of the np.unique function than the original implementation).

The cartesian_add function essentially generates all of the combinations of addition while the set.diff1d function figures out which numbers are only a part of the set of odd numbers (originally represented as list4) and therefore haven't been verified.

Schematic of the principle utilized in the "new strategy." Source

Nevertheless, this version was able to only check 100,000 numbers in approx. five seconds which was much slower when compared to the original implementation (see the table above).

Jacob G. of Futuresight Technologies also developed a C++ code to to verify Lemoine's conjecture for a certain number of digits. Please note that it requires -fopenmp flag if you want to utilize the parallelism within the program to decrease run-time.

#include <iostream>
#include <vector>
#include <set>

typedef unsigned long long ull;

using namespace std;

const ull maxNumber = 1000000000;

bool numberWorks(ull n, vector<ull>& primes, set<ull>& primeSet) {
    ull factor;
    ull remainder;
    ull primeSize = primes.size();
    ull semiPrime;
    for (ull i = 0; i < primeSize; i++) {
        factor = primes[i];
        if (factor > n) {
            return false;
        }
        remainder = n - factor;
        
        if (remainder % 2 == 0 && primeSet.find(remainder / 2) != primeSet.end()) {
            return true;
        }
    }
    return false;
}

int main() {
    bool* primes = new bool[maxNumber];
    bool* semiPrimes = new bool[maxNumber];
    
    vector<ull> primeList;
    set<ull> primeSet;
    
    for (int i = 0; i < maxNumber; i++) {
        primes[i] = true;
    }
    
    for (int i = 2; i < maxNumber; i++) {
        if (primes[i]) {
            for (int j = i; j < maxNumber; j += i) {
                primes[j] = false;
            }
            primeList.push_back(i);
            primeSet.insert(i);
        }
    }
    cout << "Computed primes!" << endl;
    
    #pragma omp parallel for num_threads(8)
    for (ull numberToCheck = 7; numberToCheck < maxNumber; numberToCheck += 2) {
        if (!numberWorks(numberToCheck, primeList, primeSet)) {
            cout << "Number failed: " << numberToCheck << endl;
        } else if ((numberToCheck - 1) % 1000000 == 0) {
            cout << (numberToCheck / (double)maxNumber) * 100 << "%" << endl;
        }
    }
    
    delete[] primes;
    delete[] semiPrimes;
}

He contributed the following speed readings:

Odd #’s up to ___ checked	Time
1,000,000	--- .6 seconds ---
10,000,000	--- 15 seconds ---
100,000,000	--- 301 seconds ---

Out of curiosity, the original C code was run again to see how efficient it was. Due to C being a compiled language, this test ended up being more successful than the python programs at and allowed the conjecture to be verified to 10^10 in a reasonable time period. Here are some times for those runs:

Odd #’s up to ___ checked	Time
1,000,000,000	--- 108 seconds ---
10,000,000,000	--- 1009 seconds ---

Therefore, it is concluded that the conjecture is now verified up to 10^10 through the original C code, and up to 2*10^9 through the new Python implementation.

Thank you to Ernie Parke and Jacob G. answering all of my questions throughout this process. 

Sources:
https://en.wikipedia.org/wiki/Lemoine%27s_conjecture
http://mathworld.wolfram.com/LevysConjecture.html
https://planetmath.org/levysconjecture
https://primes.utm.edu/curios/page.php?number_id=82&submitter=Capelle
https://www.mathsteacher.com.au/year7/ch11_ratios/04_prob/unit.htm
https://oeis.org/A046927