Brute force search for squares that are the difference of a K'th power and a cube

cubesquaresearch.py

import heapq
from fractions import gcd
import math
import sys
import time

def initialValue( maxA, b, power ):
    y = b ** power
    # Better too low than too high and miss one
    a = min( int( math.ceil( math.pow( y, 1.0 / 3.0 ) ) ),
             maxA )
    return ( - a ** 3 + b ** power, a, b )

# Values of -a^3 + b^{power} in sorted order
def powersMinusCubes( maxA, bRange, power ):
    # Want to visit -a^3 term from smallest to largest values
    # But we can prune those below zero early.
    h  = [ initialValue( maxA, b, power ) for b in bRange ]
    heapq.heapify( h )

    arraySize = sys.getsizeof( h )
    tupleSize = sys.getsizeof( h[-1] ) + sys.getsizeof( h[-1][0] ) + sys.getsizeof( h[-1][1] ) + sys.getsizeof( h[-1][2] )
    totalSize = ( arraySize + len( h ) * tupleSize ) / ( 1024 ** 2 )
    print "Heap size estimate:", totalSize, "MB"

    start = time.clock()
    
    while len( h ) >  0:
        smallest = h[0]
        yield smallest
        (y, a, b) = smallest
        if a > -maxA:
            replacement = ( y + a ** 3 - (a-1) ** 3, a-1, b )
            heapq.heapreplace( h, replacement )
        else:
            heapq.heappop( h )
            
    end = time.clock()
    print "Completed in", end-start, "seconds."

def isqrt(n):
    x = n
    y = (x + 1) // 2
    while y < x:
        x = y
        y = (x + n // x) // 2
    return x

def issquare( n ):
    if n <= 0:
        return ( False, 0 )
    
    r = isqrt( n )
    if r * r == n:
        return ( True, r )
    else:
        return ( False, 0 )

def search( aMax=10000, bChunk = 100, power=7 ):
    bStart = 1
    while True:
        bRange = xrange( bStart, bStart + bChunk )
        print "Searching b=", bStart, "...", bStart + bChunk, "a=", -aMax, "...", aMax
        for ( y, a, b ) in powersMinusCubes( maxA = aMax,
                                             bRange = bRange, 
                                             power = power ):
            ( square, r ) = issquare( y )
            if square:
                if gcd( a, b ) == 1:
                    if a < 0:
                        print "{}^2 = {}^{} - ({})^3".format( r, b, power, a )
                    else:   
                        print "{}^2 = {}^{} - {}^3".format( r, b, power, a )
        bStart += bChunk