Pratt certificate in Python
A proof that p is prime from a primitive root of p and the factorization of p - 1, repeated for each odd prime factor. It has O((log p)^2) size and needs fewer than 2 lg^2 p products to check.
- File
pratt_certificate.py, 180 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/pratt-certificate/. The authors are listed in prov/QY13TB.jsonl.
"""Pratt certificates of primality.
A certificate for n is a list of entries (p, a, factors), where p is an odd
prime, factors maps each prime q of p - 1 to its exponent, and a satisfies
a^((p-1)/2) = -1 (mod p), a^((p-1)/(2q)) != -1 (mod p) for each odd q.
Each odd q of an entry must be the p of an earlier entry, so the list is the
Lucas tree in bottom-up order, with shared primes listed once.
verify_pratt(n, cert) checks a certificate: deterministic, with
O(log n) entries and O((log n)^2) modular multiplications, so O((log n)^4)
bit operations with schoolbook arithmetic on integers of bit length log n.
build_pratt(n) builds one; its cost is dominated by factoring each p - 1.
All loops are iterative; there is no recursion.
"""
import math
import random
Entry = tuple[int, int, dict[int, int]]
def verify_pratt(n: int, cert: list[Entry]) -> bool:
if n == 2: # @line L1
return True
if n < 3 or n % 2 == 0: # @line L2
return False
proven = {2} # @line L3
for p, a, factors in cert: # @line L4
prod = 1
for q, e in factors.items():
prod *= q**e if e >= 1 else 0
if p < 3 or p % 2 == 0 or prod != p - 1: # @line L5
return False
if any(q not in proven for q in factors): # @line L6
return False
if pow(a, (p - 1) // 2, p) != p - 1: # @line L7
return False
for q in factors: # @line L8
if q != 2 and pow(a, (p - 1) // (2 * q), p) == p - 1: # @line L9
return False
proven.add(p) # @line L10
return n in proven # @line L11
def multiplications(e: int) -> int:
"""Modular multiplications of the left-to-right binary method for x^e."""
return 0 if e <= 1 else e.bit_length() - 1 + bin(e).count("1") - 1
def verify_cost(cert: list[Entry]) -> int:
"""M(p): the multiplications that verify_pratt does with the binary method."""
total = 0
for p, _, factors in cert:
total += multiplications((p - 1) // 2)
total += sum(multiplications((p - 1) // (2 * q)) for q in factors if q != 2)
return total
# Helpers for the builder. They may err (a probable prime that is composite);
# the builder then fails to find a base, so it never returns a false certificate.
def _strong_probable(n: int, a: int) -> bool:
d, s = n - 1, 0
while d % 2 == 0:
d //= 2
s += 1
x = pow(a, d, n)
if x in (1, n - 1):
return True
for _ in range(s - 1):
x = x * x % n
if x == n - 1:
return True
return False
def _probable_prime(n: int) -> bool:
if n < 2:
return False
for p in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37):
if n % p == 0:
return n == p
return all(_strong_probable(n, a) for a in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41) if a % n)
def _rho(n: int, rng: random.Random) -> int:
"""A nontrivial factor of the odd composite n (Pollard rho, Brent's cycle)."""
while True:
y, c, m = rng.randrange(1, n), rng.randrange(1, n), 128
g, r, q = 1, 1, 1
x = ys = y
while g == 1:
x = y
for _ in range(r):
y = (y * y + c) % n
k = 0
while k < r and g == 1:
ys = y
for _ in range(min(m, r - k)):
y = (y * y + c) % n
q = q * abs(x - y) % n
g = math.gcd(q, n)
k += m
r *= 2
if g == n:
g = 1
while g == 1:
ys = (ys * ys + c) % n
g = math.gcd(abs(x - ys), n)
if g != n:
return g
def factorize(m: int, rng: random.Random, small: int = 10_000) -> dict[int, int]:
"""Prime factorization of m >= 1: trial division, then Pollard rho."""
f: dict[int, int] = {}
for d in (2, 3):
while m % d == 0:
f[d] = f.get(d, 0) + 1
m //= d
d = 5
while d <= small and d * d <= m:
for dd in (d, d + 2):
while m % dd == 0:
f[dd] = f.get(dd, 0) + 1
m //= dd
d += 6
stack = [m] if m > 1 else []
while stack:
x = stack.pop()
if _probable_prime(x):
f[x] = f.get(x, 0) + 1
continue
r = math.isqrt(x)
if r * r == x:
stack += [r, r]
continue
g = _rho(x, rng)
stack += [g, x // g]
return f
def build_pratt(n: int, rng: random.Random | None = None, tries: int = 10_000) -> list[Entry] | None:
"""A Pratt certificate for n, or None if n is not prime (or no base is found).
The result lists every odd prime of the Lucas tree once, children first,
with the smallest base that works for each prime."""
rng = rng if rng is not None else random.Random()
if n == 2:
return []
if n < 3 or n % 2 == 0 or not _probable_prime(n):
return None # a composite has no certificate; this check only saves time
order: list[int] = []
facts: dict[int, dict[int, int]] = {}
stack = [n]
while stack: # depth-first, iterative; each prime is factored once
p = stack.pop()
if p in facts:
continue
facts[p] = factorize(p - 1, rng)
order.append(p)
stack += [q for q in facts[p] if q != 2 and q not in facts]
cert: list[Entry] = []
for p in sorted(order): # every q of p - 1 is smaller than p
f = facts[p]
base = 0
for a in range(2, min(p, 2 + tries)):
if pow(a, (p - 1) // 2, p) == p - 1 and all(
pow(a, (p - 1) // (2 * q), p) != p - 1 for q in f if q != 2
):
base = a
break
if base == 0:
return None
cert.append((p, base, f))
return cert if verify_pratt(n, cert) else None
test_pratt_certificate.py 172 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/pratt-certificate/. The authors are listed in prov/QY13TB.jsonl.
import math
import random
from pratt_certificate import _probable_prime as probably_prime
from pratt_certificate import build_pratt, factorize, multiplications, verify_cost, verify_pratt
LIMIT = 100_000
def sieve(n: int) -> list[bool]:
is_p = [True] * (n + 1)
is_p[0] = is_p[1] = False
for i in range(2, math.isqrt(n) + 1):
if is_p[i]:
is_p[i * i :: i] = [False] * len(is_p[i * i :: i])
return is_p
def is_prime_td(n: int) -> bool:
if n < 2:
return False
if n % 2 == 0:
return n == 2
d = 3
while d * d <= n:
if n % d == 0:
return False
d += 2
return True
CARMICHAEL = [561, 1105, 1729, 2465, 2821, 6601, 8911, 10585, 15841, 29341, 41041, 46657, 52633, 62745, 63973, 75361]
STRONG_PSP = [2047, 1373653, 25326001, 3215031751, 2152302898747, 3474749660383, 341550071728321, 3825123056546413051]
def test_example_1279() -> None:
cert = build_pratt(1279)
assert cert == [
(3, 2, {2: 1}),
(5, 2, {2: 2}),
(7, 3, {2: 1, 3: 1}),
(71, 7, {2: 1, 5: 1, 7: 1}),
(1279, 3, {2: 1, 3: 2, 71: 1}),
]
assert verify_pratt(1279, cert)
assert pow(3, 639, 1279) == 1278 and pow(3, 213, 1279) == 775 and pow(3, 9, 1279) == 498
assert pow(7, 35, 71) == 70 and pow(7, 7, 71) == 14 and pow(7, 5, 71) == 51
assert verify_cost(cert) == 48 # Crandall and Pomerance: M(1279) = 48
assert [multiplications(e) for e in (639, 213, 9, 35, 7, 5, 2, 3, 1)] == [16, 11, 4, 7, 4, 3, 1, 2, 0]
def test_small_inputs() -> None:
assert verify_pratt(2, []) and build_pratt(2) == []
for n in (-7, 0, 1, 4, 6, 100):
assert not verify_pratt(n, []) and build_pratt(n) is None
assert verify_pratt(3, [(3, 2, {2: 1})])
assert not verify_pratt(3, [])
def test_all_n_up_to_limit() -> None:
is_p = sieve(LIMIT)
rng = random.Random(1)
for n in range(LIMIT + 1):
cert = build_pratt(n, rng)
if is_p[n]:
assert cert is not None and verify_pratt(n, cert), n
# The certificate is short: fewer than lg n entries, fewer than 2 lg^2 n products.
if n > 2:
assert len(cert) < math.log2(n) and verify_cost(cert) < 2 * math.log2(n) ** 2, n
else:
assert cert is None, n
def subcertificate(f: dict[int, int], rng: random.Random) -> list:
"""Certificates of the odd primes of f, merged, each prime once."""
out: list = []
for q in sorted(f):
if q != 2:
out += [e for e in build_pratt(q, rng) if e[0] not in {x[0] for x in out}]
return out
def test_composites_are_never_certified() -> None:
"""For an odd composite p no base satisfies the conditions, so an entry
for p always fails. Try every base for all odd composites below 1000."""
rng = random.Random(0)
for p in range(9, 1000, 2):
if is_prime_td(p):
continue
f = factorize(p - 1, rng)
sub = subcertificate(f, rng)
assert verify_pratt(sub[-1][0] if sub else 2, sub)
for a in range(p):
assert not verify_pratt(p, [*sub, (p, a, f)]), (p, a)
def test_carmichael_and_strong_pseudoprimes() -> None:
rng = random.Random(2)
for n in CARMICHAEL + STRONG_PSP:
assert build_pratt(n, rng) is None, n
f = factorize(n - 1, rng)
sub = subcertificate(f, rng)
for a in (2, 3, 5, 7, n - 1):
assert not verify_pratt(n, [*sub, (n, a, f)])
def test_tampered_certificates_fail() -> None:
rng = random.Random(3)
n = 1279
good = build_pratt(n)
assert not verify_pratt(n, good[:-1]) # missing root
assert not verify_pratt(n, good[1:]) # missing child 3
assert not verify_pratt(n, [good[-1], *good[:-1]]) # child after parent
assert not verify_pratt(n, [*good[:-1], (1279, 2, good[-1][2])]) # 2 is not a primitive root
assert not verify_pratt(n, [*good[:-1], (1279, 3, {2: 1, 639: 1})]) # 639 not proven
assert not verify_pratt(n, [*good[:-1], (1279, 3, {2: 1, 3: 1, 71: 1})]) # wrong product
assert not verify_pratt(n, [*good[:-1], (1279, 3, {2: 1, 3: 2, 71: 1, 5: 0})]) # zero exponent
assert not verify_pratt(1281, [*good[:-1], (1281, 3, {2: 5, 5: 1, 8: 1})])
# A certificate for one prime does not certify another.
assert not verify_pratt(1283, good)
for _ in range(200):
bad = [list(e) for e in good]
i = rng.randrange(len(bad))
bad[i][1] = rng.randrange(2, bad[i][0])
bad_cert = [(p, a, f) for p, a, f in bad]
ok = all(
pow(a, (p - 1) // 2, p) == p - 1 and all(pow(a, (p - 1) // (2 * q), p) != p - 1 for q in f if q != 2)
for p, a, f in bad_cert
)
assert verify_pratt(n, bad_cert) == ok
def test_large_primes() -> None:
rng = random.Random(4)
for k in (61, 89, 107, 127):
n = 2**k - 1
cert = build_pratt(n, rng)
assert cert is not None and verify_pratt(n, cert) and verify_cost(cert) < 2 * k * k
n = 3 * 2**534 + 1 # 162 digits
cert = build_pratt(n, rng)
assert cert is not None and verify_pratt(n, cert) and cert[-1][0] == n
# Primes 2 * (product of primes below 10^4) + 1 with several hundred digits.
small = [p for p in range(3, 10_000) if is_prime_td(p)]
found = 0
while found < 3:
m = 2
while m.bit_length() < 900:
m *= rng.choice(small)
n = m + 1
cert = build_pratt(n, rng)
assert (cert is not None) == probably_prime(n)
if cert is not None:
assert verify_pratt(n, cert) and verify_cost(cert) < 2 * math.log2(n) ** 2
found += 1
def test_large_composites_rejected() -> None:
rng = random.Random(5)
m521, m607 = 2**521 - 1, 2**607 - 1
assert build_pratt(m521 * m607, rng) is None
assert build_pratt(m521**2, rng) is None
assert build_pratt(2**512 + 1, rng) is None
def test_random_inputs_against_trial_division() -> None:
rng = random.Random(6)
for _ in range(300):
n = rng.randrange(10**9, 10**12)
assert (build_pratt(n, rng) is not None) == is_prime_td(n)
What the next level needs
- the tests did not run in this build