Gilbreath's conjecture in Python

Question

I've written some code for testing Gilbreath's conjecture in python (2.x). The conjecture was shown (in the 90s) to be true up to n=10¹². My code uses the primesieve module to quickly generate primes and, by ignoring blocks of 0s and 2s, only calculates values in the array when necessary. However, it's quite slow; it makes heavy use of dictionary lookup, addition and deletion, and I expect would benefit from modification of this aspect.

def gilbreath(nMax):
    import primesieve
    def fill(D, j, k):
        if (j, k-1) not in D:
            fill(D, j, k-1)
        D[(j, k)] = abs(D[(j+1, k-1)]-D[(j, k-1)])
        # remove unneeded entries
        if (j-1, k) in D:
            del D[(j, k-1)]
        if (j+1, k) in D:
            del D[(j+1, k-1)]
    primes = primesieve.Iterator()
    D = {}
    depthLeft, depth, maxDepth, nDepth = -1, -1, -1, -1
    for n in xrange(1, nMax+1):
        D[(n, 0)] = int(primes.next_prime())
        j, k = n, 0
        while (D[(j, k)] > 2) or (k <= depthLeft):
            if D[(j, k)] > 2:
                depth = k
            j -= 1
            k += 1
            fill(D, j, k)
            if (j == 1) and D[(j, k)] > 1:
                print "conjecture false at n = %d" %n
        depthLeft = depth
        if depth > maxDepth:
            maxDepth, nDepth = depth, n
    print "max depth %d at n = %d of %d" %(maxDepth, nDepth, nMax)

Answer 1

12% speed-up by simplification

Your code as is takes (with input 10**5):

max depth 96 at n = 92717 of 100000
        4027699 function calls (3764972 primitive calls) in 9.303 seconds

I simplified it by removing the:

    if (j-1, k) in D:
        del D[(j, k-1)]
    if (j+1, k) in D:
        del D[(j+1, k-1)]

And it runs in:

max depth 96 at n = 92717 of 100000
         4027699 function calls (3764972 primitive calls) in 8.129 seconds

So it both runs faster and is less code to maintain, a win-win.

Answer 2

You should change the names of your variables according to the style guide:

Use the function naming rules: lowercase with words separated by underscores as necessary to improve readability.

For example, nMax should be n_max. re-writing this we get:

def gilbreath(n_max):
    import primesieve
    def fill(D, j, k):
        if (j, k-1) not in D:
            fill(D, j, k-1)
        D[(j, k)] = abs(D[(j+1, k-1)]-D[(j, k-1)])
    primes = primesieve.Iterator()
    D = {}
    depth_left, depth, max_depth, n_depth = -1, -1, -1, -1
    for n in xrange(1, n_max+1):
        D[(n, 0)] = int(primes.next_prime())
        j, k = n, 0
        while (D[(j, k)] > 2) or (k <= depth_left):
            if D[(j, k)] > 2:
                depth = k
            j -= 1
            k += 1
            fill(D, j, k)
            if (j == 1) and D[(j, k)] > 1:
                print "conjecture false at n = %d" % n
        depth_left = depth
        if depth > max_depth:
            max_depth, n_depth = depth, n
    print "max depth %d at n = %d of %d" %(max_depth, n_depth, n_max)

gilbreath(100000)

(Well, with @Caridorc's optimization).

Using cProfile I get the following:

❯ python -m cProfile gilbreath.py
max depth 96 at n = 92717 of 100000
         4027466 function calls (3764739 primitive calls) in 3.534 seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.001    0.001    0.001    0.001 __init__.py:1(<module>)
        1    0.304    0.304    3.534    3.534 gilbreath.py:1(<module>)
        1    1.192    1.192    3.230    3.230 gilbreath.py:1(gilbreath)
1963731/1701004    1.912    0.000    2.027    0.000 gilbreath.py:3(fill)
  1963731    0.115    0.000    0.115    0.000 {abs}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
   100000    0.010    0.000    0.010    0.000 {method 'next_prime' of 'primesieve._primesieve.Iterator' objects}

Now, as for optimizations, I would consider using something like Cython.

Note: From here on I deal with Cython, if you are not interested in this, then all I have is naming advice.

For ease of use I use runcython to make build development easier. To optimize this code for Cython, declare some of the types for to be integer types.

# Change file name to gilbreath.pyx
def gilbreath(int n_max):
    import primesieve
    def fill(D, int j, int k):
        if (j, k-1) not in D:
            fill(D, j, k-1)
        D[(j, k)] = abs(D[(j+1, k-1)]-D[(j, k-1)])
    primes = primesieve.Iterator()
    D = {}
    cdef int depth_left, depth, max_depth, n_depth
    depth_left, depth, max_depth, n_depth = -1, -1, -1, -1
    for n in xrange(1, n_max+1):
        D[(n, 0)] = int(primes.next_prime())
        j, k = n, 0
        while (D[(j, k)] > 2) or (k <= depth_left):
            if D[(j, k)] > 2:
                depth = k
            j -= 1
            k += 1
            fill(D, j, k)
            if (j == 1) and D[(j, k)] > 1:
                print "conjecture false at n = %d" % n
        depth_left = depth
        if depth > max_depth:
            max_depth, n_depth = depth, n
    print "max depth %d at n = %d of %d" % (max_depth, n_depth, n_max)

Secondly, for testing I make a file called test_gil.py:

import gilbreath
gilbreath.gilbreath(100000)

Finally I run:

makecython gilbreath.pyx

And:

❯ python -m cProfile test_gil.py  
max depth 96 at n = 92717 of 100000
         4 function calls in 2.185 seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.001    0.001    0.001    0.001 __init__.py:1(<module>)
        1    0.000    0.000    2.185    2.185 test_gil.py:1(<module>)
        1    2.184    2.184    2.185    2.185 {gilbreath.gilbreath}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

Shaves off some more time. About a second or so.