Euler #46

#!/usr/bin/env python

from libeuler import is_prime, primes
from math import sqrt

if __name__ == '__main__':
    P = primes(100000)
    n = 1
    while True:
        found = False
        for p in P[1:]:
            a2 = n - 0.5*(p - 1)
            if a2 < 0:
                break
            
            if sqrt(a2) == float(int(sqrt(a2))):
                print "%d -> %d + 2*%d^2" % (2*n + 1, p, sqrt(a2))
                found = True
                n += 1
                break

        if not found and not (2*n+1) in P:
            print 2*n + 1
            break