Miller-Rabin Primality Test, in depth: from Fermat to strong pseudoprimes, deterministic 64-bit bases, C and Python code, and adversarial composites

Fermat's test and why it is not enough

Fermat's little theorem says that if n is prime and a is not a multiple of n, then a to the power n - 1 is congruent to 1 mod n. The contrapositive is a compositeness test: pick a, compute pow(a, n - 1, n), and if the result is not 1, n is certainly composite. If the result is 1, n is a probable prime to base a.

The Fermat test has a fatal flaw. Carmichael numbers are composites for which a to the n - 1 is 1 mod n for every base a coprime to n. The smallest is 561 = 3 times 11 times 17. Unless you happen to pick a base sharing a factor with n, the Fermat test calls 561 prime every time, and there are infinitely many Carmichael numbers. Miller-Rabin repairs this by looking at the path to 1, not just the destination.

The idea: square roots of one

The key fact: modulo a prime p, the equation x squared = 1 has exactly two solutions, 1 and p - 1, because p divides (x - 1)(x + 1) only if it divides one of the factors. Modulo a composite with at least two distinct odd prime factors, 1 has more than two square roots. For 561, 67 squared is 4,489 = 8 times 561 + 1, so 67 is a non-trivial square root of 1.

Write n - 1 = 2 to the s times d with d odd. Then a to the n - 1 is obtained from a to the d by squaring s times. If n is prime, that sequence ends at 1 (Fermat), and the element just before the first 1 must be a square root of 1, so it must be n - 1. Hence for prime n, either a to the d is already 1, or one of a to the d, a to the 2d, up to a to the 2 to the s - 1 times d equals n - 1. If neither holds, a is a witness that n is composite. If both conditions hold for a composite n, a is called a strong liar and n a strong pseudoprime to base a.

Rabin proved that for an odd composite n, at most one quarter of the bases from 1 to n - 1 are strong liars (Monier gave the same bound independently). Each round with an independent uniformly random base therefore exposes a composite with probability at least three quarters, and k rounds miss it with probability at most 4 to the minus k. The quarter is a worst case; for most composites the fraction of liars is far smaller.

A Python implementation

The implementation is short. Trial division by small primes first is a cheap filter that removes most composites and handles tiny n. The exponentiation uses Python's three-argument pow, which performs square-and-multiply with reduction at every step; see modular exponentiation for how that works.

import secrets

SMALL_PRIMES = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37)

def _strong_probable_prime(n: int, a: int, d: int, s: int) -> bool:
    x = pow(a, d, n)
    if x == 1 or x == n - 1:
        return True
    for _ in range(s - 1):
        x = x * x % n
        if x == n - 1:
            return True
    return False                      # a is a witness: n is composite

def is_prime(n: int, rounds: int = 40) -> bool:
    if n < 2:
        return False
    for p in SMALL_PRIMES:
        if n % p == 0:
            return n == p
    d, s = n - 1, 0
    while d % 2 == 0:
        d //= 2
        s += 1
    if n < 3_317_044_064_679_887_385_961_981:
        bases = SMALL_PRIMES            # deterministic in this range
    else:
        bases = [2 + secrets.randbelow(n - 3) for _ in range(rounds)]
    return all(_strong_probable_prime(n, a, d, s) for a in bases)

Two choices in this code are deliberate. Random bases come from secrets, not random, because the probability bound assumes bases the input's author cannot predict; that matters in the adversarial section below. And the loop squares at most s - 1 times: the s-th square is a to the n - 1 itself, and if that equalled n - 1 rather than 1, Fermat's condition would fail and n would be composite anyway, so it can never rescue the round.

Worked examples

Prime first. Take n = 97, so n - 1 = 96 = 2 to the 5 times 3, giving s = 5 and d = 3. With base 5, 5 cubed mod 97 is 28. Squaring repeatedly gives 8, 64, 22 and then 96, which is n - 1, so the round passes, as it must for a prime.

Now the Carmichael number 561. Here 560 = 2 to the 4 times 35, so s = 4 and d = 35. With base 2, 2 to the 35 mod 561 is 263, and squaring gives 166, then 67, then 1. The sequence reaches 1 without passing through 560, so 2 is a witness and 561 is composite, the case the Fermat test misses. Bases 5 and 7 also pass through 67 before reaching 1. As a bonus, a non-trivial square root of 1 factors n: gcd(67 - 1, 561) = 33, and 561 = 33 times 17.

Finally a strong liar. The smallest strong pseudoprime to base 2 is 2047 = 23 times 89. Here 2046 = 2 times 1023, and because 2 to the 11 is 2048, which is 1 mod 2047, 2 to the 1023 is also 1, so base 2 alone wrongly passes 2047. Base 3 gives 3 to the 1023 mod 2047 = 1565, which is neither 1 nor 2046, so a second base exposes it immediately.

Deterministic bases for bounded n

For bounded inputs, carefully chosen fixed bases make the test deterministic. These thresholds come from exhaustive searches for the smallest strong pseudoprime to each base set. Except in the seven-base row, whose bound is the limit of the search, each number in the right column is the first composite that fools the whole set, so the set is correct for all n below it.

Apart from 2 to the 64, each threshold in this table is itself a strong pseudoprime to every base in its row; that is what makes it the boundary. For example, 1,373,653 = 829 times 1657 passes bases 2 and 3 and is caught by base 5. The 2 to the 64 case is the one most libraries use, because a 64-bit test becomes a fixed seven or twelve exponentiations with no randomness at all.

In C, the obstacle is overflow: multiplying two 64-bit residues needs 128 bits. GCC and Clang provide unsigned __int128; for hot loops, Montgomery multiplication avoids the division entirely.

#include <stdint.h>
#include <stdbool.h>

static uint64_t mulmod(uint64_t a, uint64_t b, uint64_t m) {
    return (uint64_t)((unsigned __int128)a * b % m);
}

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;
}

bool is_prime_u64(uint64_t n) {
    static const uint64_t small[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
    if (n < 2) return false;
    for (int i = 0; i < 12; i++)
        if (n % small[i] == 0) return n == small[i];
    uint64_t d = n - 1;
    int s = 0;
    while ((d & 1) == 0) { d >>= 1; s++; }
    for (int i = 0; i < 12; i++) {            /* first 12 primes: exact for 64 bits */
        uint64_t x = powmod(small[i], d, n);
        if (x == 1 || x == n - 1) continue;
        bool witness = true;
        for (int r = 1; r < s; r++) {
            x = mulmod(x, x, n);
            if (x == n - 1) { witness = false; break; }
        }
        if (witness) return false;
    }
    return true;
}

Cost and prime generation

One round costs one modular exponentiation with a log2(n)-bit exponent, about log2(n) squarings and half as many multiplications, plus at most s - 1 extra squarings. With schoolbook multiplication each step is quadratic in the bit length, so a round is cubic in the number of bits and k rounds cost O(k log cubed n). For a 2048-bit candidate one round is a few milliseconds in pure Python and far less in a native bignum library.

In prime generation most candidates are composite, and the first round rejects nearly all of them, so the cost is dominated by the single round on each failed candidate. That is why generators sieve candidates by small primes first, as in prime sieves: cheap division removes most composites before any exponentiation happens. Only the final prime pays for all k rounds.

Adversarial inputs

The 4 to the minus k bound holds when bases are random and independent of n. If an attacker chooses n and knows your bases, the bound means nothing: Arnault showed in 1995 how to build composites that are strong pseudoprimes to any given set of prime bases. The 2018 paper Prime and Prejudice by Albrecht, Massimo, Paterson and Somorovsky applied this to cryptographic libraries and found several whose primality testing, when checking parameters supplied by an untrusted party such as Diffie-Hellman groups, could be fooled by adversarially constructed composites because the bases were fixed or too few.

Practical rules follow. For numbers you generated yourself from a good random source, a modest number of random-base rounds is fine, because random candidates are almost never strong pseudoprimes. For numbers from an untrusted source, use random bases from a cryptographic generator, use enough rounds for an adversarial setting, and prefer the Baillie-PSW test, which combines a strong base-2 test with a strong Lucas test; no composite passing it is known, and it has been verified to have no counterexamples below 2 to the 64. Many libraries run Baillie-PSW plus some random Miller-Rabin rounds. For round counts in key generation, follow FIPS 186-5 Appendix B rather than inventing your own.

Failure modes

Related algorithms

Miller-Rabin is one part of a toolbox. Pollard's rho uses it to decide when a factor is prime and recursion can stop; RSA key generation calls it on every candidate; and the Solovay-Strassen test, a predecessor, uses the Jacobi symbol instead of square roots of 1 and has a weaker one-half bound per round.

What to do next

1. Implement the Python version and check it against a sieve for every n up to 10 to the 6.
2. Confirm that 561, 2047, 1,373,653 and 3,317,044,064,679,887,385,961,981 behave as the table says for each base set.
3. If you need 64-bit speed, port the C version and test it on known primes such as 2 to the 61 minus 1 and 2 to the 64 minus 59.
4. Audit any primality check that accepts numbers from outside your system, and switch it to random CSPRNG bases or Baillie-PSW.
5. In a prime generator, add small-prime sieving before Miller-Rabin and measure the speedup.
6. Read your crypto library's documentation for its round counts and compare them with FIPS 186-5 Appendix B.