Miller-Rabin Primality Test

The intuition: make a prime contradict itself

Pick a 2048-bit number at random and you need to know one thing before you can use it as an RSA key: is it prime? Dividing by every smaller number would take longer than the age of the universe — there are roughly 2²⁰⁴⁸ candidate divisors. Miller-Rabin sidesteps factoring entirely. Instead of finding a divisor, it interrogates the number with random witnesses and asks each one a question that only a prime can answer consistently.

The question rests on two facts that every prime p obeys. First, Fermat's little theorem: for any base a not divisible by p, we have a^(p−1) ≡ 1 (mod p). Second, in arithmetic modulo a prime, 1 has exactly two square roots: +1 and −1. There's no third number whose square is 1. That second fact is the lever Miller-Rabin pulls that the older Fermat test never touched.

So the test starts from a^(n−1), which should be 1 if n is prime, then walks backwards taking square roots. If at any point it sees a square root of 1 that isn't ±1, the number has just contradicted a property every prime must satisfy — and that contradiction is a certificate that n is composite. No factor needed; the structure alone gives it away.

The exact mechanism and the math

Let n > 2 be the odd number we're testing. Factor out all the twos from n − 1:

n − 1 = 2^s · d,   where d is odd, s ≥ 1

Choose a random base a with 2 ≤ a ≤ n − 2 and compute the sequence by repeated squaring:

a^d,  a^(2d),  a^(4d),  …,  a^(2^(s−1)·d)   (all mod n)

The very last square of this sequence is a^(2^s·d) = a^(n−1). If n were prime, Fermat forces that final value to 1. The base a is declared a strong liar for n — i.e. it fails to expose n — if and only if one of these holds:

• a^d ≡ 1 (mod n), or
• a^(2^r·d) ≡ −1 (mod n) for some 0 ≤ r < s.

If neither holds, a is a witness: the chain reached 1 at the top without ever passing through −1, which means somewhere a number other than ±1 squared to 1 — impossible modulo a prime. So n is definitely composite, and we stop.

The guarantee that makes this useful is the Monier-Rabin theorem (1980): for any odd composite n, at most (n − 1)/4 of the bases in [1, n−1] are strong liars. Equivalently, at least three-quarters of bases are witnesses. So a single random base fails to detect a composite with probability at most ¼; k independent bases fail with probability at most 4⁻ᵏ. The test never errs in the other direction — a true prime is never called composite.

Complexity. One round is dominated by the modular exponentiation a^d mod n, which uses O(log n) modular multiplications via square-and-multiply. With schoolbook multiplication each multiply is O(log² n) bit operations, giving O(log³ n) per round and O(k log³ n) overall. Swap in FFT-based multiplication and a round drops to Õ(log² n).

When to use Miller-Rabin

• Generating cryptographic primes. RSA, DSA, and Diffie-Hellman all need large random primes; you sieve out small factors by trial division, then confirm with a handful of Miller-Rabin rounds.
• Any "is this big number prime?" query where a vanishingly small error probability is acceptable — and you can make it exactly zero with fixed bases below known thresholds.
• Speed-critical code. It's orders of magnitude faster than the only assumption-free deterministic test, AKS, and far simpler to implement correctly than elliptic-curve primality proving (ECPP).
• Deterministic small-range tests. For 64-bit integers, testing the bases {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37} is provably exact — a fast, exact primality check for any number that fits in a machine word.

Reach for something else when you need a proof a referee will accept without trusting randomness or the generalized Riemann hypothesis — then use ECPP, which emits a verifiable certificate. And if you just want all primes up to a few million, a sieve of Eratosthenes is faster than testing each number individually.

Miller-Rabin vs other primality tests

	Miller-Rabin	Fermat	Solovay-Strassen	Baillie-PSW	AKS	ECPP
Type	Probabilistic (or det. w/ fixed bases)	Probabilistic	Probabilistic	Probabilistic (no known counterexample)	Deterministic	Deterministic, certificate
Error per round	≤ ¼	up to 1 (Carmichael)	≤ ½	0 below 2⁶⁴; none known above	0	0
Time	O(k log³ n)	O(k log³ n)	O(k log³ n)	O(log³ n) (few tests)	Õ(log⁶ n)	Õ(log⁴ n) heuristic
Fooled by Carmichael numbers	No	Yes (561, 1105, …)	No	No	No	No
Needs randomness	Yes (or GRH for full det.)	Yes	Yes	No	No	Yes (in practice)
Produces a proof	No	No	No	No	Yes (slow)	Yes (checkable)
Used in practice	OpenSSL, GnuPG, Java, Python	Rarely (pre-filter only)	Largely retired	SymPy, GMP, PARI/GP	Theory only	Record-breaking proofs

The headline takeaway: Miller-Rabin strictly dominates both Fermat (immune to Carmichael numbers) and Solovay-Strassen (a tighter ¼ error bound instead of ½, and faster). Baillie-PSW pairs one Miller-Rabin base-2 round with a strong Lucas test and has no known counterexample at all — which is why GMP and SymPy use it. AKS and ECPP are for when you need certainty, not speed.

What the numbers actually say

• 40 rounds gives error ≤ 4⁻⁴⁰ ≈ 8·10⁻²⁵. That's far below the probability of an undetected hardware bit-flip during the computation, so beyond ~40 rounds you are testing your RAM, not the number.
• For random (non-adversarial) candidates the true error is much smaller than 4⁻ᵏ. Damgård, Landrock and Pomerance showed that for a random k-round test on a random 1024-bit candidate, the probability of accepting a composite is below 2⁻⁸⁰ with just 3 rounds — which is why FIPS 186-5 allows as few as 4–5 rounds.
• Deterministic thresholds are tabulated. Testing bases {2,3} settles everything below 1,373,653; {2,3,5,7,11,13,17,19,23,29,31,37} settles everything below 3.317·10²⁴, covering all 64-bit and most 80-bit integers exactly.
• One round on a 2048-bit number is a single modular exponentiation — roughly 2048 modular multiplications of 2048-bit numbers, a few hundred microseconds on a modern CPU. Generating a full 2048-bit RSA prime, including trial division and ~5 confirming rounds, typically takes tens of milliseconds.

JavaScript implementation

Native BigInt makes this clean. The only subtlety is modular exponentiation: never compute a ** d first — it would be astronomically large. Reduce mod n at every squaring.

// a^e mod m by square-and-multiply — O(log e) multiplications.
function powMod(a, e, m) {
  a %= m;
  let result = 1n;
  while (e > 0n) {
    if (e & 1n) result = (result * a) % m;
    a = (a * a) % m;
    e >>= 1n;
  }
  return result;
}

// One round: returns true if `a` is a witness that n is composite.
function isWitness(a, n, d, s) {
  let x = powMod(a, d, n);
  if (x === 1n || x === n - 1n) return false;   // a is a strong liar this round
  for (let r = 1n; r < s; r++) {
    x = (x * x) % n;                            // square our way up to a^(n-1)
    if (x === n - 1n) return false;             // hit −1 — still could be prime
  }
  return true;                                  // never saw −1 → composite, proven
}

function isProbablePrime(n, k = 40) {
  if (n < 2n) return false;
  for (const p of [2n, 3n, 5n, 7n, 11n, 13n, 17n, 19n, 23n, 29n, 31n, 37n]) {
    if (n === p) return true;
    if (n % p === 0n) return false;             // trial division pre-filter
  }
  // Write n − 1 = 2^s · d with d odd.
  let d = n - 1n, s = 0n;
  while ((d & 1n) === 0n) { d >>= 1n; s++; }
  for (let i = 0; i < k; i++) {
    // Random base in [2, n − 2].
    const a = 2n + randomBigInt(n - 3n);
    if (isWitness(a, n, d, s)) return false;
  }
  return true;                                  // probably prime: P(error) ≤ 4^-k
}

function randomBigInt(range) {                  // uniform in [0, range)
  const bytes = (range.toString(2).length + 7) >> 3;
  let r;
  do {
    const buf = crypto.getRandomValues(new Uint8Array(bytes));
    r = buf.reduce((acc, b) => (acc << 8n) | BigInt(b), 0n);
  } while (r >= range);
  return r;
}

Two details worth flagging. The inner loop returns false (not a witness) the moment it sees n − 1 ≡ −1, because reaching −1 on the way up is exactly the "legal square root of 1" a prime is allowed. And the early return false after the loop is the heart of the test: if the squaring chain climbed all the way to a^(n−1) without ever passing through −1, some non-trivial square root of 1 must exist mod n, which no prime allows.

Python implementation

Python's arbitrary-precision int and built-in three-argument pow(a, d, n) (which does fast modular exponentiation in C) make this both short and fast. Below is the probabilistic version plus the deterministic small-base variant.

import random

def _is_witness(a, n, d, s):
    x = pow(a, d, n)                 # fast modular exponentiation in C
    if x == 1 or x == n - 1:
        return False                 # strong liar this round
    for _ in range(s - 1):
        x = x * x % n
        if x == n - 1:
            return False             # hit −1 — consistent with prime
    return True                      # never saw −1 → composite proven

def is_probable_prime(n, k=40):
    if n < 2:
        return False
    small = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37)
    for p in small:
        if n == p:
            return True
        if n % p == 0:
            return False
    # n − 1 = 2^s · d, d odd
    d, s = n - 1, 0
    while d % 2 == 0:
        d //= 2
        s += 1
    for _ in range(k):
        a = random.randrange(2, n - 1)
        if _is_witness(a, n, d, s):
            return False
    return True                      # P(error) ≤ 4^-k

def is_prime_deterministic(n):
    """Exact for all n < 3.317 * 10**24 — covers every 64-bit integer."""
    if n < 2:
        return False
    bases = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37)
    for p in bases:
        if n == p:
            return True
        if n % p == 0:
            return False
    d, s = n - 1, 0
    while d % 2 == 0:
        d //= 2
        s += 1
    return not any(_is_witness(a, n, d, s) for a in bases)

The deterministic version isn't a different algorithm — it's the same test with the random bases replaced by a fixed witness set that has been verified (by exhaustive search) to expose every composite below the stated threshold. Above the threshold the same base set is no longer guaranteed, so fall back to random bases or Baillie-PSW.

Variants worth knowing

Deterministic Miller-Rabin under GRH. Miller's original 1976 algorithm was deterministic: assuming the generalized Riemann hypothesis, testing every base a from 2 up to 2(ln n)² is provably correct, giving a polynomial-time test with no randomness. Rabin turned it into the practical randomized version in 1980 to drop the unproven GRH dependence.

Baillie-PSW. Run one Miller-Rabin round with base 2, then a strong Lucas probable-prime test. No composite is known to pass both, despite decades of searching, and it's deterministic for all 64-bit inputs. GMP, SymPy, and PARI/GP use it as their default.

Solovay-Strassen. The historical predecessor (1977) using the Jacobi symbol and Euler's criterion. Correct but strictly weaker — error ≤ ½ per round versus Miller-Rabin's ¼ — so it's been retired in favor of Miller-Rabin everywhere.

Frobenius / Quadratic Frobenius test. A stronger probabilistic test with error per round around 1/7710, so fewer rounds buy the same confidence; used where each round is expensive.

AKS and ECPP. When you need a proof rather than overwhelming confidence: AKS (Agrawal-Kayal-Saxena, 2002) is deterministic, polynomial, and assumption-free but slow; ECPP (elliptic-curve primality proving) is faster in practice and emits a short, independently checkable certificate — it's what's used to certify record-size primes.

Common bugs and edge cases

• Computing a ** d before reducing. Even a^d for a 1024-bit d has astronomically many digits. You must reduce mod n at every squaring — use square-and-multiply, never exponentiate then mod.
• Forgetting to handle even n and small n. The reduction n − 1 = 2^s · d assumes n is odd. Special-case 2 and 3, reject all other even numbers, and bail on n < 2 before factoring.
• Picking a base of 0, 1, or n − 1. Bases 1 and n − 1 are always strong liars (they satisfy the conditions trivially), so they waste a round. Draw a uniformly from [2, n − 2].
• Assuming fixed small bases work for large n. The seven- or twelve-base witness sets are exact only below their published thresholds. Use them above the bound and a Carmichael-like number can sail through — for arbitrary-size inputs you need random bases.
• Using a weak RNG for cryptographic primes. If an attacker can predict your bases — or worse, your candidate primes — the 4⁻ᵏ guarantee is meaningless. Use a CSPRNG (crypto.getRandomValues, secrets), not Math.random.
• Confusing "probable prime" with "prime." A pass means probably prime. For protocols that demand certainty against an adversary who chose the number (not just a random candidate), layer in extra rounds, Baillie-PSW, or a true proof.