Brent's variant of Pollard's rho in Python
Pollard's rho with Brent's cycle detection, which needs only one sequence, and with m differences multiplied before each gcd. Heuristically O(p^(1/2)) steps, about a quarter faster than Floyd's.
- File
brent_pollard_rho.py, 107 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/brent-pollard-rho/. The authors are listed in prov/ES28TA.jsonl.
from math import gcd
def brent_pollard_rho(n: int, c: int = 1, s: int = 2, m: int = 128) -> int:
"""Brent's variant of Pollard's rho with f(x) = x^2 + c mod n.
It compares x_{2r-2} with x_j for 3r - 1 <= j <= 4r - 2, r = 1, 2, 4, ...
(Brent's cycle detection), and multiplies m differences before each gcd.
Return a divisor g > 1 of n; g == n means failure for this c and s.
Precondition: n >= 2, 0 <= s < n, m >= 1. Raises ValueError otherwise.
Cost: heuristically O(p^(1/2)) multiplications modulo n, p the smallest prime
factor of n, with one gcd per m steps. Never more than O(n) steps.
"""
if n < 2 or not 0 <= s < n or m < 1:
raise ValueError("need n >= 2, 0 <= s < n and m >= 1")
y, r, q, g = s, 1, 1, 1 # @line L1
x = ys = s
while g == 1: # @line L2
x = y # @line L3
for _ in range(r):
y = (y * y + c) % n
k = 0 # @line L4
while k < r and g == 1: # @line L5
ys = y # @line L6
for _ in range(min(m, r - k)): # @line L7
y = (y * y + c) % n
q = q * abs(x - y) % n
g = gcd(q, n) # @line L8
k += m
r *= 2 # @line L9
if g == n: # @line L10
while True: # @line L11
ys = (ys * ys + c) % n
g = gcd(abs(x - ys), n)
if g > 1:
break
return g # @line L12
_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 _split(m: int) -> int:
"""A proper factor of an odd composite m, by Brent's rho with c = 1, 2, ...
(skipping c = m - 2). Trial division is a last resort."""
for c in range(1, 100):
if c == m - 2:
continue
g = brent_pollard_rho(m, c, 2 % m)
if g < m:
return g
d = 3
while m % d:
d += 2
return d
def factorize(n: int) -> list[int]:
"""Prime factors of n >= 1 in non-decreasing order, with multiplicity.
The factors 2 are removed first; every odd composite part is split by Brent's
rho. A loop with an explicit stack replaces recursion."""
if n < 1:
raise ValueError("n must be a positive integer")
out: list[int] = []
while n % 2 == 0:
out.append(2)
n //= 2
stack = [n]
while stack:
v = stack.pop()
if v == 1:
continue
if is_probable_prime(v):
out.append(v)
continue
g = _split(v)
stack.extend((g, v // g))
return sorted(out)
test_brent_pollard_rho.py 151 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/brent-pollard-rho/. The authors are listed in prov/ES28TA.jsonl.
import math
import random
import pytest
from brent_pollard_rho import brent_pollard_rho, factorize
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, 2199733160881]
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 next_prime(x: int) -> int:
x += 1
while not is_prime(x):
x += 1
return x
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 brent_pollard_rho(8051, 1, 2, 1) == 97
assert brent_pollard_rho(8051, 1, 2) == 97
assert brent_pollard_rho(1387, -1, 2) == 19
assert factorize(8051) == [83, 97]
def test_backtrack() -> None:
# With large batches the product can reach 0 mod n; the backtrack must
# still return a divisor, and usually a proper one.
proper = 0
for n in range(15, 3000, 2):
g = brent_pollard_rho(n, 1, 2, 1000)
assert 1 < g <= n and n % g == 0
proper += g < n
assert proper > 0
def test_divisor_always() -> None:
rng = random.Random(1975)
for _ in range(500):
n = rng.randrange(2, 10**5)
g = brent_pollard_rho(n, rng.randrange(1, 50), rng.randrange(0, n), rng.randrange(1, 20))
assert g > 1 and n % g == 0
def test_edges() -> None:
for args in ((1,), (0,), (15, 1, 15), (15, 1, -1), (15, 1, 2, 0)):
with pytest.raises(ValueError):
brent_pollard_rho(*args)
with pytest.raises(ValueError):
factorize(0)
with pytest.raises(ValueError):
factorize(-6)
assert factorize(1) == []
assert factorize(2) == [2]
assert factorize(3) == [3]
assert factorize(4) == [2, 2]
assert factorize(9) == [3, 3]
assert factorize(2**61 - 1) == [2**61 - 1]
assert factorize(2**64) == [2] * 64
assert factorize(3**40) == [3] * 40
assert factorize(1000003**2) == [1000003, 1000003]
assert factorize(10007**4) == [10007] * 4
assert factorize(2**64 + 1) == [274177, 67280421310721]
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(8051)
for bits in range(4, 33, 2):
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:
rng = random.Random(1387)
for bits in range(8, 33, 4):
p = random_prime(rng, bits)
q = next_prime(p)
assert check(p * q) == [p, q]
def test_unbalanced() -> None:
rng = random.Random(97)
for bits in range(8, 25, 4):
p, q = random_prime(rng, bits), random_prime(rng, 80 - bits)
assert check(p * q) == sorted([p, q])
def test_random_smooth() -> None:
rng = random.Random(1981)
primes = [p for p in range(2, 5000) if is_prime(p)]
for k in range(1, 30):
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