Skip to content

Instantly share code, notes, and snippets.

@inaz2
Last active January 5, 2023 12:42
Show Gist options
  • Save inaz2/481b2b415716b763fb25d79abea95211 to your computer and use it in GitHub Desktop.
Save inaz2/481b2b415716b763fb25d79abea95211 to your computer and use it in GitHub Desktop.
SQUFOF algorithm (Shanks's square forms factorization)
$ time python3 squfof.py 60766145992321225002169406923
60766145992321225002169406923 = 242950340194949 * 250117558771727
real 0m23.864s
user 0m23.850s
sys 0m0.008s
# 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))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment