The Miller-Rabin probabilistic primality test

Rabin’s 1980 randomization of Miller’s test declares \(n\) “probably prime” after \(k\) random witnesses; for a composite \(n\) each witness is fooled with probability at most \(1/4\), so the error shrinks like \(4^{-k}\). This script shows the test exposing a Carmichael number that fools Fermat’s test, measures the error rate against the number of rounds for two hard composites, and certifies a 39-digit prime far beyond trial division’s reach.

import matplotlib.pyplot as plt
import numpy as np

from mathematicskit.number_theory import fast_mod_pow, is_prime_miller_rabin, prime_factorization

A Carmichael number fools Fermat but not Miller-Rabin

n = 3215031751
print(f"{n} = {' * '.join(map(str, prime_factorization(n)))}")
print("Fermat test a^(n-1) = 1 for a = 2, 3, 5, 7:", all(fast_mod_pow(a, n - 1, n) == 1 for a in (2, 3, 5, 7)))
print("Miller-Rabin (40 random witnesses) says prime:", is_prime_miller_rabin(n))

3215031751 = 151 * 751 * 28351
Fermat test a^(n-1) = 1 for a = 2, 3, 5, 7: True
Miller-Rabin (40 random witnesses) says prime: False

The error probability shrinks like 4^-k

Run the test with k rounds on many random seeds and record how often a composite slips through.

trials = 4000
ks = np.arange(1, 7)
fig, ax = plt.subplots()
for composite in (3215031751, 9080191):
    rates = [np.mean([is_prime_miller_rabin(composite, k=int(k), seed=s) for s in range(trials)]) for k in ks]
    print(f"n = {composite}: false 'prime' rate by rounds = {np.round(rates, 4).tolist()}")
    seen = np.array(rates) > 0  # rates below 1/trials are unresolved
    ax.semilogy(ks[seen], np.array(rates)[seen], "o-", label=f"n = {composite}")
ax.semilogy(ks, 4.0 ** (-ks), "k--", label="Rabin's bound 4^-k")
ax.set_xlabel("rounds k (random witnesses)")
ax.set_ylabel("probability a composite passes")
ax.set_title("Miller-Rabin error rate vs. number of witnesses")
ax.legend()
plt.show()

n = 3215031751: false 'prime' rate by rounds = [0.2458, 0.0595, 0.013, 0.0038, 0.0, 0.0]
n = 9080191: false 'prime' rate by rounds = [0.2502, 0.0635, 0.0145, 0.0032, 0.001, 0.0002]

Certifying a large prime

Trial division would need about \(10^{19}\) divisions for \(2^{127} - 1\); Miller-Rabin needs 40 modular exponentiations.

print("2^127 - 1 is prime:", is_prime_miller_rabin(2**127 - 1))
print("2^127 + 1 is prime:", is_prime_miller_rabin(2**127 + 1))

2^127 - 1 is prime: True
2^127 + 1 is prime: False