Miller-Rabin, in depth: strong probable primes, witnesses, deterministic base sets and the bugs that break primality tests

The problem: proving compositeness cheaply

A primality test answers one question about an integer n. The naive answer, trial division up to the square root of n, costs about the square root of n operations, which is exponential in the number of bits. For a 64-bit number that is up to four billion divisions; for cryptographic sizes it is hopeless.

Miller-Rabin flips the problem. Rather than searching for a factor, it checks a property that every prime must have. If n fails the property for some base a, n is certainly composite and a is called a witness. If n passes, n is a strong probable prime to base a. Composites that pass for a base a are called strong pseudoprimes to that base, and a is then a strong liar. The whole design rests on two facts: the property is cheap to check, taking O(log n) multiplications, and no composite has many liars.

This one-sided error is the key to using the test well. A composite answer comes with a certificate, the witness, which anyone can re-check. A prime answer is a statement of probability unless you use a base set that is proven complete for your input range.

From the Fermat test to the strong test

Fermat's little theorem says that if n is prime and a is not a multiple of n, then a^n-1 = 1 (mod n). The Fermat test checks this for a few bases. It fails badly: Carmichael numbers such as 561 = 3 x 11 x 17 satisfy the congruence for every base coprime to them, so the Fermat test calls them prime no matter how many coprime bases you try.

Miller's refinement uses a second fact about primes. Modulo a prime, the equation x² = 1 has exactly two solutions, 1 and n - 1, because a prime modulus admits no zero divisors: (x - 1)(x + 1) = 0 forces one factor to be zero. Modulo a composite there can be extra square roots of 1. So instead of looking only at a^n-1, the test watches how the value gets there.

Write n - 1 = 2^s d with d odd. Compute x = a^d mod n, then square it s - 1 times. The final square would be a^n-1. If n is prime, this chain must either start at 1, or reach n - 1 at some step, after which every later square is 1. Any other pattern means either Fermat failed or the chain passed through a square root of 1 that is neither 1 nor n - 1, and either way n is composite.

Worked examples, traced by hand

Three traces, each easy to reproduce with Python's three-argument pow.

n = 97, a = 5. 96 = 2⁵ x 3, so s = 5 and d = 3. The chain is 5³ = 125 = 28, then 28² = 784 = 8, then 64, then 64² = 4096 = 22, then 22² = 484 = 96. The value 96 is n - 1, so 97 passes base 5, as a prime must.

n = 561, a = 2. 560 = 2⁴ x 35. The chain is 263, 166, 67, 1. It reaches 1 without passing through 560, so 2 is a witness and 561 is composite, even though 2⁵⁶⁰ = 1 mod 561 and the Fermat test is fooled. Better still, 67 is a nontrivial square root of 1, so gcd(67 - 1, 561) = 33 and gcd(67 + 1, 561) = 17 are factors. A witness of this kind leaks a factorisation for free.

n = 221 = 13 x 17. 220 = 2² x 55. With a = 174 the chain is 47, 220: it hits n - 1, so 174 is a strong liar. With a = 137 the chain is 188, 205: no n - 1, so 137 is a witness. Among the bases coprime to 221, the only strong liars are 1, 21, 47, 174, 200 and 220, six out of 192.

A correct implementation

The Python version below is complete and deterministic for every n below 2⁶⁴ when used with its default bases, and still correct, though no longer proven complete, beyond that. Three lines carry most of the correctness: the small-prime filter, which also handles n up to 37 where some bases would equal n; the reduction a %= n; and the skip when the reduced base is zero.

SMALL_PRIMES = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37)
# Seven bases that make the test deterministic for every n < 2^64.
BASES_64 = (2, 325, 9375, 28178, 450775, 9780504, 1795265022)

def is_probable_prime(n: int, bases=BASES_64) -> bool:
    if n < 2:
        return False
    for p in SMALL_PRIMES:              # cheap filter, also handles tiny n
        if n % p == 0:
            return n == p
    d, s = n - 1, 0
    while d % 2 == 0:                   # n - 1 = 2^s * d with d odd
        d //= 2
        s += 1
    for a in bases:
        a %= n                          # bases may exceed n: reduce them
        if a == 0:
            continue                    # a multiple of n tells us nothing
        x = pow(a, d, n)
        if x == 1 or x == n - 1:
            continue
        for _ in range(s - 1):
            x = x * x % n
            if x == n - 1:
                break
        else:
            return False                # a is a witness: definitely composite
    return True

Python integers are arbitrary precision, so x * x % n never overflows. In C, C++, Java long or Rust u64 it does: two 64-bit residues multiply to a 128-bit product. Use a 128-bit intermediate, Montgomery multiplication, or a big-integer library.

#include <stdint.h>
typedef unsigned __int128 u128;

static uint64_t mulmod(uint64_t a, uint64_t b, uint64_t m) {
    return (uint64_t)((u128)a * b % m);     /* a*b overflows 64 bits otherwise */
}
static uint64_t powmod(uint64_t a, uint64_t e, uint64_t m) {
    uint64_t r = 1; a %= m;
    while (e) { if (e & 1) r = mulmod(r, a, m); a = mulmod(a, a, m); e >>= 1; }
    return r;
}

Exponentiation by squaring makes each round cost about log₂ n squarings plus at most as many multiplications. With schoolbook multiplication on b-bit numbers, one round is O(b³) bit operations, and k rounds cost O(k b³), comfortably fast for thousands of bits. Compare that with the exponential cost of trial division using the growth-rate vocabulary: the gap is the reason probabilistic tests won.

How wrong can a probable prime be?

Monier and Rabin proved independently that for odd composite n greater than 9, at most a quarter of the bases in [1, n - 1] are strong liars. Choose k bases independently and uniformly at random, and the probability that a composite survives all of them is at most 4^-k. Forty rounds give at most 2^-80. That bound is per input and holds for every composite, including ones an adversary picked, provided the bases are truly random and secret from the adversary.

For random inputs, such as candidates during key generation, the real error is far below the bound because most composites have very few liars; the 221 example has 6 liars out of 192 coprime bases, not 48. Standards for key generation account for this and specify round counts as a function of key size. Do not hard-code a number from memory; take it from the standard you must comply with, such as FIPS 186-5 for RSA.

Deterministic base sets

For bounded inputs, exhaustive searches have found small base sets with no strong pseudoprimes below a threshold. These turn the test into a proof for that range. The base sets below are published results, checked here against a sieve for every n below 10⁶ and against known pseudoprimes.

Two theoretical results frame these tables. Assuming the generalized Riemann hypothesis, testing every base up to 2(ln n)² is deterministic; Miller's 1976 test rested on GRH, and Bach later proved this explicit bound. Without any hypothesis, the AKS algorithm proves primality in polynomial time, but it is far slower in practice and is not used in production.

Beyond 64 bits the common production choice is Baillie-PSW: one Miller-Rabin round to base 2 followed by a strong Lucas probable-prime test. No counterexample is known, and it has been verified to have none below 2⁶⁴. Its Lucas half selects a parameter with the Jacobi symbol, the same tool covered in the article on Legendre and Jacobi symbols.

Failure modes

• Unreduced bases. With the seven-base set, a small n such as 5 gives a = 325 = 0 (mod 5). Without reduction and the zero skip, pow(325, d, 5) is 0 and the test calls a prime composite. Unit tests that start at large n never see it.
• Edge inputs. 0, 1, 2, 3, even numbers and negative values each need explicit handling.
• Overflow. A 64-bit x * x % n silently wraps for n above 2³², producing a test that is right for small inputs and randomly wrong for large ones.
• Fixed bases on adversarial input. A deterministic base list is a public target. Albrecht, Massimo, Paterson and Somorovsky showed in Prime and Prejudice (2018) that several cryptographic libraries could be fed crafted composites that passed their fixed-base or few-round checks. When an attacker supplies the number, for example Diffie-Hellman parameters received from a peer, use random bases or Baillie-PSW.
• Weak randomness. Bases from a predictable generator are as bad as fixed bases against an adversary. Use the system cryptographic generator.

Using the test in real systems

Random primes for keys. By the prime number theorem, the density of primes near 2¹⁰²⁴ is about 1 / ln(2¹⁰²⁴), roughly 1 in 710, or 1 in 355 among odd numbers. Generators therefore test hundreds of candidates per prime. Almost all are rejected cheaply: trial division by a table of small primes discards most candidates, and only survivors pay for Miller-Rabin. One base-2 round rejects nearly every remaining composite, so the extra rounds run mostly on true primes.

Hash table capacities. Tables that use prime bucket counts, discussed in the article on hash tables, need the next prime above a size. With 64-bit sizes and the deterministic base set, a loop of next-odd-then-test is both exact and fast.

Factoring after compositeness. Miller-Rabin tells you a number is composite but rarely gives a factor; the 561 gcd trick is a lucky special case. When you need factors, the next step is Pollard rho, whose cycle detection is the tortoise-and-hare idea from Floyd cycle detection.

For untrusted input of any size, randomize the bases per call:

import secrets

def is_probable_prime_random(n: int, rounds: int = 40) -> bool:
    # For inputs an adversary may have chosen: fresh random bases each call.
    if n < 5:
        return n in (2, 3)
    bases = [2 + secrets.randbelow(n - 3) for _ in range(rounds)]   # a in [2, n-2]
    return is_probable_prime(n, bases)

What to do next

1. Reproduce the 97, 561 and 221 traces with pow(a, d, n) so the squaring chain is concrete.
2. Implement the test with the seven 64-bit bases, including the reduction and zero skip, and compare it with a sieve for every n below 10⁶.
3. Add regression cases for 2047, 1,373,653, 3,215,031,751, the inputs 0 to 40 and a large known prime such as 2⁶¹ - 1.
4. If you write it in a fixed-width language, add a 128-bit or Montgomery multiply and test n just below 2⁶⁴.
5. Audit every place your code accepts a number from outside: switch those paths to random bases or Baillie-PSW.
6. For key generation, use your cryptographic library rather than your own code, and look up the round count in the standard you follow.