Williams' p+1 algorithm in Python
Computes g = gcd(V_M(A) - 2, N) for a Lucas sequence V and a smooth exponent M. It finds a prime p of N when p + 1 or p - 1, chosen by the seed A, has only small prime power factors.
- File
williams_p_plus_1.py, 153 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/williams-p-plus-1/. The authors are listed in prov/XKCQGC.jsonl.
from math import gcd, isqrt
def primes_up_to(b: int) -> list[int]:
"""The primes p <= b, by the sieve of Eratosthenes in O(b log log b) time."""
if b < 2:
return []
sieve = bytearray([1]) * (b + 1)
sieve[0] = sieve[1] = 0
for i in range(2, isqrt(b) + 1):
if sieve[i]:
sieve[i * i :: i] = bytearray(len(range(i * i, b + 1, i)))
return [i for i in range(b + 1) if sieve[i]]
def lucas_v(v: int, m: int, n: int) -> int:
"""V_m(v) mod n for the Lucas sequence V_0 = 2, V_1 = v, V_(k+1) = v V_k - V_(k-1).
Precondition: m >= 1, n >= 2. It keeps the pair (V_j, V_(j+1)) and reads the
binary digits of m from the top (a binary Lucas chain).
Cost: 2 multiplications modulo n per binary digit of m."""
x, y = v % n, (v * v - 2) % n # @line K1
for b in bin(m)[3:]: # @line K2
if b == "1":
x, y = (x * y - v) % n, (y * y - 2) % n # @line K3
else:
x, y = (x * x - 2) % n, (x * y - v) % n # @line K4
return x # @line K5
def williams_p_plus_1(n: int, bound: int, seed: int) -> int:
"""Stage 1 of Williams' p + 1 method. Return g = gcd(V_M(seed) - 2, n), where
M is the product of the largest powers q^e <= bound of the primes q <= bound.
Precondition: n >= 3, bound >= 1, 0 <= seed < n. The result g divides n. If a
prime p | n does not divide seed^2 - 4 and p - (D/p) divides M, with D =
seed^2 - 4 and (D/p) the Legendre symbol, then p divides g.
Cost: about 2.9 * bound multiplications modulo n, plus a sieve up to bound.
"""
if n < 3 or bound < 1 or not 0 <= seed < n:
raise ValueError("need n >= 3, bound >= 1 and 0 <= seed < n")
v = seed # @line L1
for q in primes_up_to(bound): # @line L2
qe = q # @line L3
while qe * q <= bound:
qe *= q
v = lucas_v(v, qe, n) # @line L4
g = gcd(v - 2, n) # @line L5
return g # @line L6
_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
_PRIMES: list[int] = primes_up_to(1 << 12)
_PRIME_CAP = 1 << 22
def _primes_to(b: int) -> list[int]:
"""A cached list of primes that covers every prime up to min(b, 2^22)."""
global _PRIMES
b = min(b, _PRIME_CAP)
if _PRIMES[-1] < b:
_PRIMES = primes_up_to(max(b, 2 * _PRIMES[-1]))
return _PRIMES
def _split(m: int) -> int:
"""A proper factor of a composite m with no prime factor below 100.
Stage 1 runs with the exponents q, one prime at a time (q^k <= m), and a gcd
after each step, over all primes q <= sqrt(m). If p is the smallest prime
factor of m and p does not divide seed^2 - 4, then p - (D/p) divides the
exponent by the end, so the gcd becomes p or m. When it jumps from 1 to m in
one step, the next seed is tried. Trial division is the last resort."""
limit = isqrt(m)
primes = _primes_to(limit)
for seed in range(3, 40):
g = gcd(seed * seed - 4, m)
if 1 < g < m:
return g
v = seed
for q in primes:
if q > limit:
break
qk = q
while qk <= m:
v = lucas_v(v, q, m)
g = gcd(v - 2, m)
if 1 < g < m:
return g
if g == m:
break
qk *= q
if g == m:
break
d = 101
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 100, then each composite part is split
by the p + 1 method. A loop with an explicit stack replaces recursion."""
if n < 1:
raise ValueError("n must be a positive integer")
out: list[int] = []
for d in range(2, 100):
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_williams_p_plus_1.py 181 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/williams-p-plus-1/. The authors are listed in prov/XKCQGC.jsonl.
import math
import random
import pytest
from williams_p_plus_1 import factorize, lucas_v, primes_up_to, williams_p_plus_1
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 plus_smooth_prime(rng: random.Random, bits: int, b: int) -> int:
"""A random prime p of about the given size with p + 1 = 2 * (distinct odd
primes <= b), so every prime power in p + 1 is at most b."""
small = primes_up_to(b)[1:]
while True:
m, used = 2, set()
while m.bit_length() < bits:
r = rng.choice(small)
if r not in used:
used.add(r)
m *= r
if is_prime(m - 1):
return m - 1
def legendre(a: int, p: int) -> int:
a %= p
if a == 0:
return 0
return 1 if pow(a, (p - 1) // 2, p) == 1 else -1
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_lucas_v() -> None:
for v in range(-3, 8):
seq = [2, v]
for _ in range(60):
seq.append(v * seq[-1] - seq[-2])
for m in range(1, 60):
assert lucas_v(v, m, 1000003) == seq[m] % 1000003
# The composition rule V_m(V_k(v)) = V_(mk)(v).
assert lucas_v(lucas_v(6, 4, 5959), 15, 5959) == lucas_v(6, 60, 5959)
def test_example() -> None:
assert lucas_v(6, 4, 5959) == 1154
assert lucas_v(1154, 3, 5959) == 497
assert lucas_v(497, 5, 5959) == 5607
assert williams_p_plus_1(5959, 5, 6) == 59
assert williams_p_plus_1(5959, 5, 3) == 1
assert williams_p_plus_1(5959, 5, 5) == 101
assert williams_p_plus_1(5959, 5, 9) == 5959
assert math.gcd(pow(2, 60, 5959) - 1, 5959) == 1
assert factorize(5959) == [59, 101]
def test_divisibility_theorem() -> None:
rng = random.Random(1982)
for _ in range(400):
n = rng.randrange(3, 10**6)
b = rng.randrange(1, 40)
seed = rng.randrange(0, n)
g = williams_p_plus_1(n, b, seed)
assert n % g == 0
m = 1
for q in primes_up_to(b):
qe = q
while qe * q <= b:
qe *= q
m *= qe
d = seed * seed - 4
for p in primes_up_to(1000)[1:]:
if n % p == 0 and d % p and m % (p - legendre(d, p)) == 0:
assert g % p == 0
def test_edges() -> None:
for args in ((2, 5, 0), (15, 0, 3), (15, 5, -1), (15, 5, 15)):
with pytest.raises(ValueError):
williams_p_plus_1(*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(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:
check(n)
def test_random_semiprimes() -> None:
rng = random.Random(1403)
for bits in range(4, 19):
for _ in range(3):
p, q = random_prime(rng, bits), random_prime(rng, bits)
assert check(p * q) == sorted([p, q])
def test_smooth_p_plus_1_grow() -> None:
rng = random.Random(299)
for bits in range(16, 49, 4):
p = plus_smooth_prime(rng, bits, 500)
q = random_prime(rng, bits)
hits = 0
for seed in range(3, 23):
g = williams_p_plus_1(p * q, 500, seed)
assert g in (1, p, q, p * q)
if legendre(seed * seed - 4, p) == -1:
assert g % p == 0
hits += g == p
assert hits > 0
if bits <= 32:
assert check(p * q) == sorted([p, q])
def test_random_smooth() -> None:
rng = random.Random(2047)
primes = primes_up_to(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