In [14]:

    

import math
from itertools import count
from time import clock

sqrt = math.sqrt

def pe44(prove=False, printit=0):
    """ D, a, b, S are pentagon numbers
        a + b  = S 
        |b- a| = D 
        phase 1: search for a solution (0.2 s) -> min_D
                 given D, b -> test a = b - D and S = 2 * b - D
                 a < b < S, D < b < S in the method below
                 a may be larger or smaller than D, but is positive
                 As a < b, we can test if a in pentagonals.
                 S is larger than b and we have to use is_pent
        If prove is True then phase 2 proves that there is no solution with
          a D lower than min_D:
            phase 2: try to find a solution with a lower D
            if D is smaller than b - b_prev it can be removed from Ds
            search stops if D < b - b_prev for all D in Ds (Ds is empty)
            total time 250 sec. In C this should be only a few secs.
    """
    t0 = clock()
    def pent(n): return n * (3 * n - 1) // 2
    def is_pent(N): return ((sqrt(24 * N + 1) + 1) / 6.0).is_integer()
    ps = set() # keeps pentagonal numbers smaller than b
    Ds = []    # keeps the list of possible D
    min_D = None
    b = 0
    # for b_index in xrange(1, 2000000): 
    for b_index in count(1): 
        b_prev = b
        b = pent(b_index)
        delta = b - b_prev 
        # any D smaller than delta can be removed, as a=b-D cannot be pentagonal
        if Ds and Ds[0] < delta: 
            Ds.pop(0)
        if printit and b_index % printit == 0:
            print b_index, clock() - t0, len(Ds)
        # if min_D is found and no Ds we can stop
        if min_D and not (prove and Ds):
            break
        for D in Ds:
            if b - D in ps and is_pent(2 * b - D): 
                min_D = D
                while Ds[-1] >= min_D: # from now on only check D < min_D
                    Ds.pop()
                break
        ps.add(b)
        if not min_D: 
            Ds.append(b)
    return min_D


    

In [15]:

    

pe44(prove=False, printit=100000)


    

    
Out[15]:





5482660

In [16]:

    

pe44(prove=True, printit=100000)


    

    




100000 42.87 1464
200000 93.46 1279
300000 143.68 1137
400000 190.73 1017
500000 233.38 911
600000 273.39 816
700000 311.4 728
800000 342.1 646
900000 369.92 570
1000000 394.94 497
1100000 417.16 428
1200000 436.61 362
1300000 453.32 299
1400000 467.48 238
1500000 477.11 179
1600000 484.21 122
1700000 488.82 67
1800000 490.94 14






    
Out[16]:





5482660