Shanks' square forms factorization in Python
Walks the continued fraction of sqrt(kN) until a denominator at an even step is a square, then walks a second cycle to a symmetric point that shows a factor. Heuristically O(N^(1/4)) steps.
- File
shanks_squfof.py, 139 lines- Level
- I0 not checked yet
- Toolchain
- Python 3.14, pytest 9
- Test framework
- pytest
- Run the tests
python3 -m pytest
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/shanks-squfof/. The authors are listed in prov/QBJ183.jsonl.
from math import gcd, isqrt
MULTIPLIERS = (1, 3, 5, 7, 11, 15, 21, 33, 35, 55, 77, 105, 165, 231, 385, 1155)
def squfof(n: int, k: int = 1, limit: int = 0) -> int:
"""Shanks' square forms factorization of n with the multiplier k.
Precondition: n >= 3 odd, not a perfect square, k >= 1 squarefree, and k * n
not a perfect square. The forward and backward cycles each make at most
`limit` steps; the default is 4 * ceil(sqrt(2 * sqrt(k * n))) + 8.
Returns a proper factor of n, or 1 if this multiplier fails.
Cost: O(limit) steps, each one division and a few products of numbers below
2 sqrt(k n); heuristically O(n^(1/4)) steps.
"""
if n < 3 or n % 2 == 0 or k < 1:
raise ValueError("need odd n >= 3 and k >= 1")
d = k * n # @line L1
a0 = isqrt(d)
if a0 * a0 == d:
raise ValueError("k * n must not be a perfect square")
if limit <= 0:
limit = 4 * (isqrt(2 * a0) + 1) + 8
p, q_prev, q = a0, 1, d - a0 * a0 # @line L2
i = 1
while i < limit:
b = (a0 + p) // q # @line L3
p_new = b * q - p
q_prev, q = q, q_prev + b * (p - p_new)
p = p_new
i += 1
if i % 2 == 1: # @line L4
continue
r = isqrt(q)
if r * r != q:
continue
b = (a0 - p) // r # @line L5
s_p = b * r + p
s_prev, s = r, (d - s_p * s_p) // r
for _ in range(limit): # @line L6
b = (a0 + s_p) // s # @line L7
s_new = b * s - s_p
if s_new == s_p: # @line L8
break
s_prev, s = s, s_prev + b * (s_p - s_new)
s_p = s_new
else:
continue
g = gcd(n, s_p) # @line L9
if 1 < g < n:
return g # @line L10
return 1 # @line L11
_BASES = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41)
def is_probable_prime(n: int) -> bool:
"""Strong probable-prime test to the 13 prime bases up to 41.
It is a proof of primality for n < 3317044064679887385961981 (about 3.3 * 10^24)."""
if n < 2:
return False
for p in _BASES:
if n % p == 0:
return n == p
d, s = n - 1, 0
while d % 2 == 0:
d, s = d // 2, s + 1
for a in _BASES:
x = pow(a, d, n)
if x in (1, n - 1):
continue
for _ in range(s - 1):
x = x * x % n
if x == n - 1:
break
else:
return False
return True
_SMALL_PRIMES = [d for d in range(2, 1 << 16) if all(d % e for e in range(2, isqrt(d) + 1))]
def _split(m: int) -> int:
"""A proper factor of a composite m with no prime factor below 2^16.
A perfect square is split by its square root. Otherwise SQUFOF runs with
the multipliers 1, 3, 5, ..., 1155 in turn, skipping those that share a
factor with m or make k * m a square. Trial division is the last resort."""
r = isqrt(m)
if r * r == m:
return r
for k in MULTIPLIERS:
if gcd(k, m) > 1:
continue
s = isqrt(k * m)
if s * s == k * m:
continue
g = squfof(m, k)
if g > 1:
return g
d = 65537
while m % d:
d += 2
return d
def factorize(n: int) -> list[int]:
"""Prime factors of n >= 1 in non-decreasing order, with multiplicity.
Trial division removes the primes below 2^16, then each composite part is
split by SQUFOF. A loop with an explicit stack replaces recursion. SQUFOF
costs about m^(1/4) steps for a part m, so this is meant for n up to about
2^80 once the small factors are gone."""
if n < 1:
raise ValueError("n must be a positive integer")
out: list[int] = []
for d in _SMALL_PRIMES:
if d * d > n:
break
while n % d == 0:
out.append(d)
n //= d
stack = [n]
while stack:
m = stack.pop()
if m == 1:
continue
if is_probable_prime(m):
out.append(m)
continue
g = _split(m)
stack.extend((g, m // g))
return sorted(out)
test_shanks_squfof.py 159 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/shanks-squfof/. The authors are listed in prov/QBJ183.jsonl.
import math
import random
import pytest
from shanks_squfof import MULTIPLIERS, factorize, squfof
CARMICHAEL = [561, 1105, 1729, 2465, 2821, 6601, 8911, 10585, 15841, 29341, 41041, 46657,
52633, 62745, 63973, 75361, 101101, 115921, 126217, 162401, 172081, 188461,
252601, 278545, 294409, 314821, 334153, 340561, 399001, 410041, 449065,
488881, 512461, 3215031751]
def is_prime(n: int) -> bool:
"""Deterministic Miller-Rabin for n < 3.3 * 10^24, written here independently."""
assert n < 3317044064679887385961981
if n < 2:
return False
small = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41)
for p in small:
if n % p == 0:
return n == p
d, s = n - 1, 0
while d % 2 == 0:
d, s = d // 2, s + 1
for a in small:
x = pow(a, d, n)
if x in (1, n - 1):
continue
for _ in range(s - 1):
x = x * x % n
if x == n - 1:
break
else:
return False
return True
def random_prime(rng: random.Random, bits: int) -> int:
while True:
p = rng.getrandbits(bits) | (1 << (bits - 1)) | 1
if is_prime(p):
return p
def check(n: int) -> list[int]:
f = factorize(n)
assert math.prod(f) == n
assert all(is_prime(p) for p in f)
assert f == sorted(f)
return f
def small_primes(b: int) -> list[int]:
return [p for p in range(2, b + 1) if is_prime(p)]
def test_examples() -> None:
assert squfof(11111) == 41
assert squfof(13231) == 101
assert factorize(11111) == [41, 271]
assert factorize(13231) == [101, 131]
def test_returns_divisor() -> None:
rng = random.Random(1975)
for _ in range(2000):
n = rng.randrange(3, 10**7) | 1
if math.isqrt(n) ** 2 == n:
continue
for k in (1, 3, 5, 7, 11):
if math.isqrt(k * n) ** 2 == k * n:
continue
g = squfof(n, k)
assert g == 1 or (1 < g < n and n % g == 0)
if is_prime(n):
assert g == 1
def test_squfof_all_small() -> None:
# Every odd composite non-square n < 30000 is split by some multiplier.
for n in range(9, 30001, 2):
if is_prime(n) or math.isqrt(n) ** 2 == n:
continue
for k in MULTIPLIERS:
if math.gcd(k, n) > 1 or math.isqrt(k * n) ** 2 == k * n:
continue
g = squfof(n, k)
if g > 1:
assert n % g == 0 and g < n
break
else:
raise AssertionError(n)
def test_edges() -> None:
for args in ((2,), (1,), (10,), (9,), (15, 0), (7, 7)):
with pytest.raises(ValueError):
squfof(*args)
with pytest.raises(ValueError):
factorize(0)
assert factorize(1) == []
assert factorize(2) == [2]
assert factorize(3) == [3]
assert factorize(2**61 - 1) == [2**61 - 1]
assert factorize(2**40) == [2] * 40
assert factorize(101**2) == [101, 101]
assert factorize(101**3) == [101] * 3
assert factorize(1000003**2) == [1000003, 1000003]
assert factorize(10007**3) == [10007] * 3
assert factorize(1000003 * 1000033) == [1000003, 1000033]
def test_all_small() -> None:
for n in range(1, 100001):
f = factorize(n)
assert math.prod(f) == n
assert all(is_prime(p) for p in f)
def test_carmichael() -> None:
for n in CARMICHAEL + [2199733160881]:
check(n)
def test_random_semiprimes() -> None:
rng = random.Random(1403)
for bits in range(4, 27):
for _ in range(3):
p, q = random_prime(rng, bits), random_prime(rng, bits)
assert check(p * q) == sorted([p, q])
def test_close_primes() -> None:
p = 1000003
for _ in range(10):
q = p + 2
while not is_prime(q):
q += 2
assert check(p * q) == [p, q]
p = q
def test_products_of_large_primes() -> None:
rng = random.Random(1155)
for k in range(1, 5):
f = sorted(random_prime(rng, 17) for _ in range(k))
assert check(math.prod(f)) == f
def test_random_smooth() -> None:
rng = random.Random(2047)
primes = small_primes(3000)
for k in range(1, 25):
f = sorted(rng.choice(primes) for _ in range(k))
assert check(math.prod(f)) == f
What the next level needs
- the tests did not run in this build