Lehman's factorization method in Python
Trial division to N^(1/3), then Fermat's method on 4kN for every multiplier k up to N^(1/3), each over a short range of x. It factors N or proves it prime in O(N^(1/3)) operations.
- File
lehman_factorization.py, 106 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/lehman-factorization/. The authors are listed in prov/DT65XZ.jsonl.
from math import gcd, isqrt
def icbrt(n: int) -> int:
"""The largest integer c with c^3 <= n, for n >= 0, by Newton's method."""
if n < 8:
return 1 if n > 0 else 0
c = 1 << ((n.bit_length() + 2) // 3)
while True:
d = (2 * c + n // (c * c)) // 3
if d >= c:
break
c = d
while c * c * c > n:
c -= 1
while (c + 1) ** 3 <= n:
c += 1
return c
def lehman_factorization(n: int) -> int:
"""Lehman's method. Return a divisor d of n with 1 < d < n, or n if n is prime.
Precondition: n >= 22. Raises ValueError otherwise.
Cost: O(n^(1/3)) arithmetic operations on numbers of the size of n: trial
division to n^(1/3), then about n^(1/3) perfect-square tests.
"""
if n < 22:
raise ValueError("n must be at least 22")
c = icbrt(n) # @line L1
for d in range(2, c + 1):
if n % d == 0:
return d
for k in range(1, c + 2): # @line L2
a = isqrt(4 * k * n) # @line L3
if a * a < 4 * k * n:
a += 1
while True:
r = a * a - 4 * k * n # @line L4
if r * r * r >= n * n:
break
b = isqrt(r) # @line L5
if b * b == r:
return gcd(a + b, n)
a += 1 # @line L6
return n # @line L7
_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
def factorize(n: int) -> list[int]:
"""Prime factors of n >= 1 in non-decreasing order, with multiplicity.
Parts below 22 are factored by trial division. Larger parts are tested with
Miller-Rabin and, when composite, split by Lehman's method. A loop with an
explicit stack replaces recursion."""
if n < 1:
raise ValueError("n must be a positive integer")
out: list[int] = []
stack = [n]
while stack:
m = stack.pop()
if m == 1:
continue
if m < 22:
d = 2
while m % d:
d += 1
out.append(d)
stack.append(m // d)
continue
if is_probable_prime(m):
out.append(m)
continue
d = lehman_factorization(m)
stack.extend((d, m // d))
return sorted(out)
test_lehman_factorization.py 137 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/lehman-factorization/. The authors are listed in prov/DT65XZ.jsonl.
import math
import random
import pytest
from lehman_factorization import factorize, icbrt, lehman_factorization
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 test_example() -> None:
assert lehman_factorization(2581) == 89
assert factorize(2581) == [29, 89]
def test_icbrt() -> None:
for x in range(0, 3000):
c = icbrt(x)
assert c**3 <= x < (c + 1) ** 3
for e in range(1, 200, 7):
x = 3**e
c = icbrt(x)
assert c**3 <= x < (c + 1) ** 3
def test_split_or_prime() -> None:
for n in range(22, 30000):
d = lehman_factorization(n)
if is_prime(n):
assert d == n
else:
assert 1 < d < n and n % d == 0
def test_edges() -> None:
for bad in (-1, 0, 1, 21):
with pytest.raises(ValueError):
lehman_factorization(bad)
with pytest.raises(ValueError):
factorize(0)
assert factorize(1) == []
assert factorize(2) == [2]
assert factorize(3) == [3]
assert factorize(4) == [2, 2]
assert factorize(21) == [3, 7]
assert lehman_factorization(23) == 23
assert lehman_factorization(25) == 5
assert factorize(2**61 - 1) == [2**61 - 1]
assert factorize(2**50) == [2] * 50
assert factorize(1000003**2) == [1000003] * 2
assert factorize(10007**3) == [10007] * 3
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:
check(n)
def test_random_semiprimes_grow() -> None:
rng = random.Random(1974)
for bits in range(4, 21):
for _ in range(3):
p, q = random_prime(rng, bits), random_prime(rng, bits)
assert check(p * q) == sorted([p, q])
def test_unbalanced_and_close() -> None:
rng = random.Random(637)
for bits in range(10, 17, 2):
p, q = random_prime(rng, bits), random_prime(rng, 2 * bits)
assert check(p * q) == sorted([p, q])
p = random_prime(rng, 2 * bits)
q = p + 2
while not is_prime(q):
q += 2
assert check(p * q) == [p, q]
def test_random_smooth() -> None:
rng = random.Random(1643)
primes = [p for p in range(2, 2000) if is_prime(p)]
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