fermfactor.py

from math import sqrt, ceil
from gmpy import mpz, is_square, is_prime

fermat = 5959
rsa64 = 16748810522526493651
n = rsa64

report_batch = 1000000

def fermat(N):
    cnt = 0
    a = ceil(sqrt(N))
    b2 = mpz(a * a - N)
    while not is_square(b2):
        a += 1
        b2 = mpz(a * a - N)
        cnt += 1

        if cnt % report_batch == 0:
            print a, sqrt(b2), b2

    print a, b2
    return a - sqrt(b2), a + sqrt(b2)


p1, p2 = fermat(n)

print int(p1), int(p2), int(p1 * p2), int(p1 * p2) == n, is_prime(mpz(p1)), is_prime(mpz(p2))