Lucas primality test in Python
Proves n prime from the prime factors of n - 1. For each prime q dividing n - 1 it finds a base a with a^(n-1) ≡ 1 and a^((n-1)/q) ≢ 1 (mod n).
- File
lucas.py, 97 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/lucas-primality-test/. The authors are listed in prov/12XYK5.jsonl.
"""The Lucas primality test (n - 1 test), in the form of Brillhart, Lehmer and Selfridge.
lucas_primality_test(n, primes) decides whether n is prime, given the distinct
prime factors of n - 1. For each such prime q it looks for a base a with
a^(n-1) = 1 (mod n) and a^((n-1)/q) != 1 (mod n). The answer is always right:
true is a proof of primality, false comes with a Fermat witness.
Cost, with l binary digits in n and m distinct primes q: each base tried costs
two modular exponentiations, O(l) products modulo n, so O(l^3) bit operations
with schoolbook arithmetic. For a prime n the number of bases tried for q is
the least a that is not a q-th power residue, usually 2 or 3. For a composite
n each search stops at the latest at the smallest prime factor of n. Integers
are of arbitrary precision (Python int).
"""
def distinct_prime_factors(m: int) -> list[int]:
"""The distinct prime factors of m >= 1, by trial division (for small m only)."""
if m < 1:
raise ValueError("m must be at least 1")
primes = []
d = 2
while d * d <= m:
if m % d == 0:
primes.append(d)
while m % d == 0:
m //= d
d += 1 if d == 2 else 2
if m > 1:
primes.append(m)
return primes
def _check_factors(n: int, primes: list[int] | tuple[int, ...]) -> None:
"""Raise ValueError unless primes are distinct and n - 1 is a product of their powers."""
m = n - 1
if len(set(primes)) != len(primes):
raise ValueError("the primes must be distinct")
for q in primes:
if q < 2 or m % q != 0:
raise ValueError("each q must be at least 2 and divide n - 1")
while m % q == 0:
m //= q
if m != 1:
raise ValueError("the primes must include every prime factor of n - 1")
def lucas_witnesses(n: int, primes: list[int] | tuple[int, ...]) -> dict[int, int] | None:
"""For a prime n, a dict q -> a with a^(n-1) = 1 and a^((n-1)/q) != 1 (mod n).
Return None if n is composite. primes must be the distinct prime factors of
n - 1; this is checked, except that the primality of each q is the caller's
responsibility.
"""
if n < 3: # @line L1
return {} if n == 2 else None
_check_factors(n, primes)
found: dict[int, int] = {}
for q in primes: # @line L2
a = 2 # @line L3
while True: # @line L4
if pow(a, n - 1, n) != 1: # @line L5
return None
if pow(a, (n - 1) // q, n) != 1: # @line L6
found[q] = a
break
a += 1 # @line L7
return found # @line L8
def lucas_primality_test(n: int, primes: list[int] | tuple[int, ...]) -> bool:
"""True if n is prime, given the distinct prime factors of n - 1 (n >= 2)."""
return lucas_witnesses(n, primes) is not None
def is_prime_lucas(n: int) -> bool:
"""Primality of n by the Lucas test, factoring n - 1 by trial division (small n only)."""
if n < 3:
return n == 2
return lucas_primality_test(n, distinct_prime_factors(n - 1))
def verify_lucas_witnesses(n: int, primes: list[int] | tuple[int, ...], witnesses: dict[int, int]) -> bool:
"""Check a certificate from lucas_witnesses: one exponentiation pair per prime q."""
if n == 2:
return True
if n < 3:
return False
_check_factors(n, primes)
for q in primes:
a = witnesses.get(q)
if a is None or pow(a, n - 1, n) != 1 or pow(a, (n - 1) // q, n) == 1:
return False
return True
lucas_primality_test.py 20 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/lucas-primality-test/. The authors are listed in prov/12XYK5.jsonl.
"""Entry point for the Lucas primality test; the code is in lucas.py."""
from lucas import (
distinct_prime_factors,
is_prime_lucas,
lucas_primality_test,
lucas_witnesses,
verify_lucas_witnesses,
)
__all__ = [
"distinct_prime_factors",
"is_prime_lucas",
"lucas_primality_test",
"lucas_witnesses",
"verify_lucas_witnesses",
]
test_lucas_primality_test.py 153 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/lucas-primality-test/. The authors are listed in prov/12XYK5.jsonl.
import math
import random
import pytest
from lucas_primality_test import (
distinct_prime_factors,
is_prime_lucas,
lucas_primality_test,
lucas_witnesses,
verify_lucas_witnesses,
)
LIMIT = 100_000
CARMICHAEL = (561, 1105, 1729, 2465, 2821, 6601, 8911, 41041, 825265, 321197185, 5394826801)
# Smallest strong pseudoprimes to the first 1 to 7 prime bases.
PSI = (2047, 1373653, 25326001, 3215031751, 2152302898747, 3474749660383, 341550071728321)
SMALL_PRIMES = [p for p in range(2, 10_000) if all(p % d for d in range(2, math.isqrt(p) + 1))]
def sieve(limit: int) -> bytearray:
is_p = bytearray([1]) * (limit + 1)
is_p[0] = is_p[1] = 0
for i in range(2, math.isqrt(limit) + 1):
if is_p[i]:
is_p[i * i :: i] = bytearray(len(is_p[i * i :: i]))
return is_p
def trial_division_is_prime(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
def miller_rabin(n: int, rng: random.Random, rounds: int = 30) -> bool:
"""An independent probable prime test for large n."""
if n < 4:
return n >= 2
if n % 2 == 0:
return False
d, s = n - 1, 0
while d % 2 == 0:
d //= 2
s += 1
for _ in range(rounds):
x = pow(rng.randrange(2, n - 1), 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 test_example_1009_and_561() -> None:
assert distinct_prime_factors(1008) == [2, 3, 7]
assert lucas_witnesses(1009, [2, 3, 7]) == {2: 11, 3: 2, 7: 2}
assert pow(2, 504, 1009) == 1 and pow(11, 504, 1009) == 1008
assert pow(2, 336, 1009) == 374 and pow(2, 144, 1009) == 105
assert lucas_primality_test(1009, [7, 3, 2]) is True
# 561: the base 2 passes the Fermat test and 2^280 = 1, then 3 divides 561.
assert pow(2, 560, 561) == 1 and pow(2, 280, 561) == 1 and pow(3, 560, 561) == 375
assert lucas_witnesses(561, [2, 5, 7]) is None
def test_against_trial_division_up_to_limit() -> None:
is_p = sieve(LIMIT)
for n in range(LIMIT + 1):
assert is_prime_lucas(n) == bool(is_p[n]), n
for n in range(0, 3000):
assert is_prime_lucas(n) == trial_division_is_prime(n)
def test_certificates_verify() -> None:
for n in range(3, 5000):
primes = distinct_prime_factors(n - 1)
w = lucas_witnesses(n, primes)
if w is None:
assert not trial_division_is_prime(n)
else:
assert trial_division_is_prime(n)
assert verify_lucas_witnesses(n, primes, w)
assert not verify_lucas_witnesses(1009, [2, 3, 7], {2: 2, 3: 2, 7: 2})
assert not verify_lucas_witnesses(561, [2, 5, 7], {2: 5, 5: 2, 7: 2})
def test_hard_composites() -> None:
for n in CARMICHAEL + PSI:
assert lucas_primality_test(n, distinct_prime_factors(n - 1)) is False, n
for n in (4, 6, 9, 25, 3**10, 2**20, 7919**2, 65537 * 65539):
assert lucas_primality_test(n, distinct_prime_factors(n - 1)) is False
def test_large_primes_and_composites_with_factored_n_minus_1() -> None:
rng = random.Random(5)
for digits in (30, 100, 200, 300):
primes_found = composites_found = tries = 0
while tries < 3000 and (primes_found < 1 or composites_found < 3):
tries += 1
primes = sorted(set(rng.sample(SMALL_PRIMES[1:], 4)) | {2})
m = 2
while m < 10 ** (digits - 1):
m *= rng.choice(primes)
n = m + 1
qs = distinct_prime_factors_from(m, primes)
probable = miller_rabin(n, rng)
assert lucas_primality_test(n, qs) is probable
if probable:
w = lucas_witnesses(n, qs)
assert w is not None and verify_lucas_witnesses(n, qs, w)
primes_found += 1
else:
composites_found += 1
assert primes_found >= 1 and composites_found >= 3, digits
# Fermat primes and a few Proth-style primes with n - 1 a power of 2 times a prime.
for k in range(5):
assert lucas_primality_test(2 ** (2**k) + 1, [2]) is True
for k in range(5, 10):
assert lucas_primality_test(2 ** (2**k) + 1, [2]) is False
assert lucas_primality_test(3 * 2**189 + 1, [2, 3]) is True
assert lucas_primality_test(3 * 2**190 + 1, [2, 3]) is miller_rabin(3 * 2**190 + 1, rng)
def distinct_prime_factors_from(m: int, primes: list[int]) -> list[int]:
return [q for q in primes if m % q == 0]
def test_bad_factor_lists() -> None:
with pytest.raises(ValueError):
lucas_primality_test(1009, [2, 3]) # 7 is missing
with pytest.raises(ValueError):
lucas_primality_test(1009, [2, 3, 7, 5]) # 5 does not divide 1008
with pytest.raises(ValueError):
lucas_primality_test(1009, [2, 2, 3, 7])
with pytest.raises(ValueError):
distinct_prime_factors(0)
assert lucas_primality_test(2, []) is True
assert lucas_primality_test(1, []) is False
assert lucas_primality_test(0, []) is False
What the next level needs
- the tests did not run in this build