foo.py

def ispower(a,b): #a==b^i ?
    while a>1:
        if a%b==0:a/=b
        else:return False
    return True

def beal(max_base, max_power):
    bases, powers, table, pow = initial_data(max_base, max_power)
    for x in bases:         
        flag2=flag3=False # for 2^m+3^n=z^r  x=2^i or x=3^i ?
        if ispower(x,2):flag2=True
        elif ispower(x,3):flag3=True        
        powx = pow[x]
        if x % 10 == 0: print 'x =', x,time.asctime()# add time        
        for y in bases:
            if y > x or gcd(x,y) > 1: continue
            if flag2:
                if ispower(y,3):continue # for 2^m+3^n=z^r
            elif flag3 and ispower(y,2):continue            
            powy = pow[y]
            for m in powers:
                if x==max_base and m==max_power:continue # no need check
                xm = powx[m]
                for n in powers:
                    sum = xm + powy[n]
                    try:
                        r = table[sum]
                        report(x, m, y, n, nth_root(sum, r), r)
                    except KeyError:
                        pass
                
def initial_data(max_base, max_power):
    bases = range(2, max_base+1) # no need 1 
    powers = range(3, max_power+1)
    table = {}
    pow = [None] * (max_base+1)
    for z in bases:
        pow[z] = [None] * (max_power+1)
        for r in powers:
            zr = long(z) ** r
            pow[z][r] = zr
            table[zr] = r
    print 'initialized %d table elements' % len(table)
    return bases, powers, table, pow  

def report(x, m, y, n, z, r):
    x, y, z = map(long, (x, y, z))
    assert min(x, y, z) > 0 and min(m, n, r) > 2 
    assert x ** m + y ** n == z ** r
    print '%d ^ %d + %d ^ %d = %d ^ %d = %s' % \
          ( x,   m,   y,   n,   z,   r,  z**r)
    if gcd(x,y) == gcd(x,z) == gcd(y, z) == 1: 
        raise 'a ruckus: SOLUTION!!'

def gcd(x, y):
    while x:
        x, y = y % x, x
    return y

def nth_root(base, n): 
    return long(round(base ** (1.0/n)))

import time
def timer(b, p):
    start = time.clock()
    beal(b, p)
    secs = time.clock() - start
    print  {'secs': secs, 'mins': secs/60, 'hrs': secs/60/60}

import cProfile, pstats, StringIO
pr = cProfile.Profile()
pr.enable()
timer(200, 20)
pr.disable()
s = StringIO.StringIO()
sortby = 'cumulative'
ps = pstats.Stats(pr, stream=s).sort_stats(sortby)
ps.print_stats()
print s.getvalue()