fermat_factorial.py

import gmpy


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


def fermat(n, verbose=True):
    a = isqrt(n)
    b2 = a * a - n
    b = isqrt(n)
    count = 0
    while b * b != b2:
        if verbose:
            pass
            #   print('Trying: a=%s b2=%s b=%s' % (a, b2, b))
        a = a + 1
        b2 = a * a - n
        b = isqrt(b2)  # int(b2**0.5)
        count += 1
    p = a + b
    q = a - b
    assert n == p * q
    print('a=', a)
    print('b=', b)
    print('p=', p)
    print('q=', q)
    print('pq=', p * q)
    return p, q


def fermat_factor(n, n1, n2):
    a = gmpy.sqrt(n) + 1
    b = a * a - n
    bsq = gmpy.sqrt(b)
    while bsq * bsq != b or a + bsq == n1 or a + bsq == n2:
        a += 1
        b = a * a - n
        bsq = gmpy.sqrt(b)
    return a + bsq
mersenne = [2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127,
            521, 607, 1279, 2203, 2281, 3217, 4253, 4423, 9689]
n = int("798321817573328185527646107613495929846147444322791353283989998016" +
        "2788028361090036128124997317580506991621017956050649707513252490" +
        "2086881120372213626641879468491936860976686933630869673826972619" +
        "9383219515991467448076533010760265779495796183315027763039834855" +
        "6604648543103954170846714140826022009859276124501067859234750189" +
        "4176269580510459729633673468068467144199744563731826362102608811" +
        "0334008878137547802826280994434901700160878386069980174904566013" +
        "1580244856777241162382628174724566095424541378151979429533619755" +
        "5688543537992197142258053220453757666537840276416475602759374950" +
        "715283890232230741542737319569819793988431443")
for n1 in mersenne:
    for n2 in mersenne:
        m1 = (2 ** n1) - 1
        m2 = (2 ** n2) - 1
        if m1 * m2 == n:
            print "Match! ", m1, m2
print n