SQUFOF algorithm (Shanks's square forms factorization)

gistfile1.txt

$ time python3 squfof.py 60766145992321225002169406923
60766145992321225002169406923 = 242950340194949 * 250117558771727

real    0m23.864s
user    0m23.850s
sys     0m0.008s

squfof.py

# squfof.py

import sys
import gmpy2
from math import gcd, prod

# https://en.wikipedia.org/wiki/Shanks%27s_square_forms_factorization
def squfof(n):
    if n % 2 == 0:
        return 2
    if gmpy2.is_square(n):
        raise Exception('n is perfect square')

    L = 2*gmpy2.isqrt(2*gmpy2.isqrt(n))
    B = 3*L
    ks = [1, 3, 5, 7, 11, 3*5, 3*7, 3*11, 5*7, 5*11, 7*11, 3*5*7, 3*5*11, 3*7*11, 5*7*11, 3*5*7*11]
    for k in ks:
        Pinit = gmpy2.isqrt(k*n)
        P0 = Pinit
        Q0 = 1
        Q1 = k*n - P0*P0
        for i in range(2, B):
            b = (Pinit + P0)//Q1
            P1 = b*Q1 - P0
            Q2 = Q0 + b*(P0 - P1)
            if i % 2 == 0 and gmpy2.is_square(Q2):
                break
            P0, Q0, Q1 = P1, Q1, Q2
        else:
            continue

        q = gmpy2.isqrt(Q2)
        b0 = (Pinit - P1)//q
        P0 = b0*q + P1
        Q0 = q
        Q1 = (k*n - P0*P0)//Q0
        while True:
            b = (Pinit + P0)//Q1
            P1 = b*Q1 - P0
            Q2 = Q0 + b*(P0 - P1)
            if P0 == P1:
                break
            P0, Q0, Q1 = P1, Q1, Q2

        d = gcd(n, P1)
        if 1 < d < n:
            return d


if __name__ == '__main__':
    n = int(sys.argv[1])
    p = squfof(n)
    if p:
        print('{} = {} * {}'.format(n, p, n//p)) 
    else:
        print('{} is prime'.format(n))