NUM Numbers and algebra Algorithm
The Miller–Rabin test
Also known as: Miller–Rabin primality test · Rabin–Miller test · strong probable prime test · Fermat test · Carmichael numbers · strong pseudoprimes
How to make sure in a tenth of a millisecond that a twenty-digit number is prime without finding a single divisor of it. The test calls witnesses, numbers called bases, and checks their testimony: a prime cannot confuse them, and a composite is caught by at least three quarters of the bases.
Is 18,446,744,073,709,551,557 prime? Trial division up to the square root is two billion divisions, and by our measurement Python would spend about a minute on them. The Miller–Rabin test answers in a tenth of a millisecond, and finds no divisor along the way. It calls witnesses, numbers called bases, and checks their testimony. A prime cannot confuse them, and a composite will almost certainly give itself away. Below the trial runs in your browser. Choose a defendant (among them are record holders that fool thirteen witnesses in a row) or enter your own number.
In each row one base a testifies in two ways. The column “Fermat” is the old test: the remainder of an−1 modulo n. For a prime n it is always 1, as Fermat’s little theorem says. The column “chain” is the same power, but computed with stops: first ad, then its square, the square of the square and so on up to an−1. Look at 561: the Fermat test is silent on most bases, while the chain gives the number away on the very first base. On ψ₁₃ thirteen bases in a row lie, and only 43 tells the truth. And on a Carmichael number of 61 digits the Fermat test will not notice the trick with any base coprime to n, while the chain not only catches it but also hands over a divisor.
The idea
Modulo a prime, one has exactly two square roots: 1 and −1. Write n − 1 = 2s·d with odd d and square ad repeatedly until reaching an−1. If n is prime, the chain ends in 1, and that 1 can only appear at once or right after −1. If a 1 in the chain arose from anything else, or the chain never reached 1 at all, n is composite, and the base a has proved it.
The Fermat test looks only at the last link of the chain. It is fooled by any composite n with an−1 ≡ 1, and such numbers exist for every base. The chain demands more: not just to arrive at one, but to arrive there by the right road. It can be fooled too, but much less easily.
Call a base a a witness if the chain gives the number away: ad ≢ 1 and none of the numbers ad, a2d, …, a2s−1d equals −1 modulo n. Otherwise n is said to be a strong probable prime to base a. A composite number that passed the check is called a strong pseudoprime to base a, and the base itself a strong liar.
The interrogation step by step
Take the prime 97. Then 96 = 2⁵ · 3: s = 5, d = 3. Base 5: 5³ = 125 ≡ 28, then squares: 28² ≡ 8, 8² = 64, 64² ≡ 22, 22² ≡ 96 ≡ −1, and the last square (−1)² = 1. The chain 28, 8, 64, 22, −1, 1 reached one through −1, as a prime should.
Now 2047 = 23 · 89. Here 2046 = 2 · 1023, two links. Base 3: 31023 ≡ 1565, neither 1 nor −1, and no more squares can be made: the next link is already an−1. A witness. Base 2: 21023 ≡ 1, because 2¹¹ = 2048 ≡ 1, and 1023 is divisible by 11. Two lies. The number 2047 is the smallest strong pseudoprime to base 2, and in all it has 242 liars among 2046 bases, under 12%.
Each base costs one fast exponentiation, about log₂ n modular multiplications, and at most s squarings. The bases are tried in turn, and the first witness ends the trial: “composite” is a proof, not a guess. If no base objected, the answer is “probably prime”, and how probably is decided by the next section.
Why a prime never fails and a composite often does
If p is prime and x² ≡ 1 (mod p), then x ≡ 1 or x ≡ −1 (mod p).
x² − 1 = (x − 1)(x + 1) is divisible by p. A prime that divides a product divides one of the factors; this is Euclid’s lemma. So p divides x − 1 or x + 1.
If n is an odd prime, every base a from 1 to n − 1 passes the check.
Denote the links of the chain x₀ = ad, x₁ = x₀², …, xs = an−1. By Fermat’s little theorem xs ≡ 1. If x₀ ≡ 1, the check is passed. Otherwise take the first link xi equal to one; i ≥ 1, and xi−1 ≢ 1. But xi−1² = xi ≡ 1, and by the lemma xi−1 ≡ −1, a link numbered from 0 to s − 1, exactly what the check requires.
The converse is what makes the test useful: a composite has few liars. Let us prove that they are at most half. We will need one observation about the remainders invertible modulo n, those coprime to n; there are φ(n) of them.
Let a set H of invertible remainders modulo n contain, together with any two elements, their product, and together with each element, its inverse. If at least one invertible remainder b is not in H, then H has at most φ(n)/2 elements.
Multiply each element h of H by b. Different h give different bh: an invertible b can be cancelled. And no bh lies in H: otherwise b = (bh)·h⁻¹ would lie in H too. So outside H there are at least as many invertible remainders as in H, and |H| ≤ φ(n)/2.
If n is an odd composite, then among the bases from 1 to n − 1 there are at most φ(n)/2 strong liars, that is, fewer than half.
First: a liar is coprime to n. For it an−1 ≡ 1, that is, a · an−2 ≡ 1, and an inverse of a exists. It remains in each case to find a set H that contains all the liars and fits the half lemma.
Case 1: n is divisible by the square of a prime p. Let H be the remainders with an−1 ≡ 1, the Fermat liars; products and inverses do not leave it, and the strong liars are in it. Write n = pk·m, where k ≥ 2 and m is not divisible by p, and by the Chinese remainder theorem take b ≡ 1 + p (mod pk), b ≡ 1 (mod m). Such a b is coprime to n. By the binomial theorem (1 + p)n−1 ≡ 1 + (n − 1)p (mod p²): all the other terms are divisible by p². But p divides n, so n − 1 ≡ −1 (mod p) and bn−1 ≡ 1 − p (mod p²), not one. So bn−1 ≢ 1 modulo n either, and b is not in H.
Case 2: n is a product of distinct primes. There are at least two of them, so n = q·m, where q is prime, m > 1, and they are coprime. The base n − 1 gives (−1)d = −1, so there are j for which a2jd ≡ −1 for some a; take the largest such j and such an a. Let H be the remainders x with x2jd ≡ ±1. Products and inverses do not leave it. Every liar is in it: if xd ≡ 1, then x2jd ≡ 1 too; if x2id ≡ −1, then i ≤ j by the choice of j, and x2jd is −1 squared j − i times, that is, ±1. Now take b ≡ a (mod q), b ≡ 1 (mod m). Then b2jd ≡ −1 modulo q and ≡ 1 modulo m. Were it +1 modulo n, we would get −1 ≡ 1 (mod q); were it −1, we would get 1 ≡ −1 (mod m). Both moduli are odd and greater than two, so neither is possible, and b is not in H.
In both cases the half lemma gives at most φ(n)/2 liars.
Hence the main property of the test. A random base catches a composite n with probability greater than ½, and k independent bases are all wrong at once with probability less than 2−k. The error is one-sided: “composite” is always right, and doubt remains only in the answer “prime”. This is a Monte Carlo algorithm, but with an error that can be made smaller than the probability of a hardware fault.
In fact the bound is twice as strict. In 1980 Michael Rabin and, independently, Louis Monier proved that an odd composite n > 9 has at most φ(n)/4 liars, and the error of k rounds is less than 4−k. Monier also found an exact formula for their number. Let n have exactly r distinct prime divisors, p − 1 = 2sp·dp with odd dp, and ν the smallest of the sp. Then the number of strong liars is
$$S(n) = \left(1 + \frac{2^{r\nu} - 1}{2^r - 1}\right) \prod_{p \mid n} \gcd(d, d_p),$$and the number of Fermat liars F(n) = ∏ gcd(n − 1, p − 1). For 2047 = 23 · 89 the formula gives (1 + 3/3) · 11 · 11 = 242, exactly as we counted. The derivation is longer than our proof: it rests on the fact that the invertible remainders modulo a power of an odd prime are powers of a single remainder, and it can be found in Monier’s 1980 paper and in the book by Crandall and Pomerance, “Prime Numbers”. We do not repeat it but check it: the listing below compares the formula with a direct count over all bases for all odd composites up to a thousand.
An atlas of liars
A quarter is the worst case. How do real numbers lie? The scene below factors every odd composite number up to a million, 421,502 of them, and by Monier’s formula computes what share of bases lies to it. Each dot is one number, the height is the share of liars, on a logarithmic scale. Tap a dot to learn what the number is, and send it to trial in the first scene.
The lower edge of the cloud is a dense curve 2/(n − 1): numbers to which only 1 and n − 1 lie. These two bases pass the check for any odd n, and taking them is pointless. There are 261,279 such numbers up to a million, 62%. Each next curve higher up is its own small set of liars: 6, 10, 14, 18 bases. The mean share of liars over all 421,502 numbers is 0.0235%, more than a thousand times less than a quarter.
Twenty numbers rise right up to the bound of ¼. Nineteen of them are built the same way, p·(2p − 1), where p and 2p − 1 are prime and p leaves remainder 3 on division by 4: 15 = 3 · 5, 91 = 7 · 13, 703 = 19 · 37 and on up to 954,271 = 691 · 1381. The twentieth is 8911 = 7 · 19 · 67, a Carmichael number of three primes of the form 4k + 3. The record holders of the first scene are from the same family p(2p − 1): ψ₁₂ = 399,165,290,221 · 798,330,580,441 and ψ₁₃ = 1,287,836,182,261 · 2,575,672,364,521. This is no coincidence: the more liars a number has, the better the chance that all bases in a row lie.
Switch the scene to the Fermat test. The cloud rises, and above the line ½ there are only Carmichael numbers. This is the same half lemma: if n has at least one Fermat witness coprime to n, the Fermat liars, a set closed under multiplication, take up at most half. The exception is numbers with no such witnesses at all.
The difference between the worst case and the average explains why standards count rounds differently. For a number chosen by an adversary the honest estimate is 4−k: OpenSSL does at least 64 rounds for numbers up to 2048 bits, which is 2−128, and 128 rounds for longer numbers. For a random candidate during key generation the error is incomparably smaller, and the OpenSSL documentation estimates the chance of getting a composite instead of a 2048-bit prime at 2−648.
Carmichael numbers: why Fermat alone is not enough
A composite n for which an−1 ≡ 1 (mod n) for all a coprime to n is called a Carmichael number. The Fermat test catches it only if a base sharing a divisor happens to be chosen, and for a number made of three twenty-digit primes that is almost impossible. Korselt’s criterion helps to recognise such numbers from the factorisation.
If a composite n is not divisible by the square of a prime and for every prime divisor p the number p − 1 divides n − 1, then n is a Carmichael number.
Let a be coprime to n and p a prime divisor of n. By Fermat’s little theorem ap−1 ≡ 1 (mod p), and n − 1 = (p − 1)t, so an−1 = (ap−1)t ≡ 1 (mod p). The number an−1 − 1 is divisible by every prime divisor of n, they are distinct, so it is divisible by their product too, that is, by n.
The converse is also true (a Carmichael number has no square divisors and p − 1 divides n − 1), but its proof needs the existence of a primitive root, and we skip it. Alwin Korselt stated the criterion in 1899 without knowing a single example. The first seven numbers, 561, 1105, 1729, 2465, 2821, 6601 and 8911, had been found back in 1885 by Václav Šimerka, but his work went unnoticed, and the numbers were named after Robert Carmichael, who described 561 in 1910. The criterion shows at once that a Carmichael number has at least three prime divisors: if n = pq and p < q, then n − 1 = p(q − 1) + (p − 1), and q − 1 would have to divide the smaller number p − 1.
The criterion makes it possible to build Carmichael numbers of any size. Jack Chernick noticed in 1939: if 6k + 1, 12k + 1 and 18k + 1 are prime, their product is a Carmichael number. Indeed, the product equals 1296k³ + 396k² + 36k + 1, and n − 1 = 36k(36k² + 11k + 1) is divisible by 6k, 12k and 18k. With k = 1 it gives 7 · 13 · 19 = 1729. The button “a new Carmichael” in the first scene looks for such a k of twenty digits, usually in a few thousand or tens of thousands of tries, and gets a number of sixty-odd digits against which the Fermat test is powerless.
The table below gives how many numbers of each kind are below 10k. We recounted them ourselves up to a billion, and the count matched the published tables (Pomerance, Selfridge and Wagstaff, 1980; Pinch for Carmichael numbers).
| below | primes | Fermat pseudoprimes to base 2 | strong pseudoprimes to base 2 | Carmichael numbers |
|---|---|---|---|---|
| 10³ | 168 | 3 | 0 | 1 |
| 10⁴ | 1229 | 22 | 5 | 7 |
| 10⁵ | 9592 | 78 | 16 | 16 |
| 10⁶ | 78,498 | 245 | 46 | 43 |
| 10⁷ | 664,579 | 750 | 162 | 105 |
| 10⁸ | 5,761,455 | 2057 | 488 | 255 |
| 10⁹ | 50,847,534 | 5597 | 1282 | 646 |
Base 2 alone, up to a billion, is wrong 1282 times against fifty million primes. Of the 646 Carmichael numbers below a billion, only 43 pass the strong check to base 2. There are infinitely many Carmichael numbers, as Alford, Granville and Pomerance proved in 1994, but for the strong test they are nothing special: they have at most a quarter liars, like everyone else.
How many bases are really needed
For numbers of bounded length randomness can be removed altogether. Denote by ψk the smallest composite that passes the strong test to all of the first k prime bases. If n < ψk and it passed them all, it is prime, without any probabilities. The values of ψk were sought by searches with clever pruning over several decades:
| bases | ψ: the first number that fools them | found by |
|---|---|---|
| 2 | 2047 | — |
| 2, 3 | 1,373,653 | Pomerance, Selfridge and Wagstaff, 1980 |
| 2, 3, 5 | 25,326,001 | |
| 2, 3, 5, 7 | 3,215,031,751 | |
| 2 … 11 | 2,152,302,898,747 | Jaeschke, 1993 |
| 2 … 13 | 3,474,749,660,383 | |
| 2 … 17 and 2 … 19 | 341,550,071,728,321 | |
| 2 … 23, 2 … 29 and 2 … 31 | 3,825,123,056,546,413,051 | Jiang and Deng, 2014 |
| 2 … 37 | 318,665,857,834,031,151,167,461 | Sorenson and Webster, 2015 |
| 2 … 41 | 3,317,044,064,679,887,385,961,981 |
The number ψ₁₂ ≈ 3.2·10²³ is greater than 2⁶⁴ ≈ 1.8·10¹⁹, so the first twelve primes, from 2 to 37, are enough for any 64-bit number. Jim Sinclair’s 2011 set is shorter: the seven bases 2, 325, 9375, 28178, 450775, 9780504, 1795265022 check all n < 2⁶⁴. Such sets are verified against the list of all pseudoprimes to base 2 below 2⁶⁴, computed by Jan Feitsma. Michal Forišek and Jakub Jančina went further in 2015: for 32-bit numbers one round is enough if the base is chosen from a precomputed table by a hash of the number itself, and for 64-bit numbers two or three.
For numbers of any length a deterministic variant exists only under a hypothesis. Gary Miller showed in 1976 that if the extended Riemann hypothesis is true, it is enough to check all bases up to O(ln² n); Eric Bach refined the bound in 1990 to 2 ln² n. The hypothesis is unproven, and in practice random bases are relied on.
Code
The first program is the trial from the scene in code. The function strong builds the whole chain so that it can be shown; in a working test the loop is broken off as soon as −1 is met. The numbers are long, so Go uses math/big, and Python an ordinary int and the built-in pow(a, d, n), that is, fast modular exponentiation.
number digits fermat strong bases 2..41 say with 43 too 2^64 - 59 20 13/13 13/13 probably prime probably prime 561 3 10/13 0/13 composite (witness 2) composite (witness 2) 3215031751 10 13/13 6/13 composite (witness 11) composite (witness 11) psi12 24 12/13 12/13 composite (witness 41) composite (witness 41) psi13 25 13/13 13/13 probably prime composite (witness 43) chernick 61 13/13 2/13 composite (witness 2) composite (witness 2)
The Fermat test is fooled by the 61-digit Carmichael number on all thirteen bases, the strong test only on two. And then the evidence itself: the chain for base 2 is x x 1 1 1, so a one appeared from a number not equal to ±1. That number x is a square root of one, and n divides (x − 1)(x + 1) without dividing either bracket. So gcd(x − 1, n) and gcd(x + 1, n) are non-trivial divisors. The program gets 12k + 1 = 120,000,000,000,000,138,301 and the product of the other two primes. A test that looked for no divisors found a divisor for free.
The second program checks Monier’s formula against a direct count and builds a small atlas, up to 300,000, to stay within a fraction of a second. The share of numbers to which only 1 and n − 1 lie is 63.1% there, and in the scene up to a million 62.0%: it falls, but very slowly. Eric Bach and Andrew Shallue showed in 2013 that the share of such numbers tends to zero, at a rate inverse to the triple logarithm.
The third is a working version for 64-bit numbers with Sinclair’s bases. In Python there is nothing to fear here, but in Go the product of two 64-bit remainders does not fit in a machine word. The function bits.Mul64 returns all 128 bits of the product, and bits.Div64 divides them by n. Before the exponentiation comes trial division by twelve small primes: 85% of all numbers have a divisor no greater than 37, and they are weeded out for a small fraction of the cost of one round. The program checks itself against a sieve up to a hundred thousand, catches ψ₉, before which the bases 2 … 23 give up, finds the eight largest primes that fit in 64 bits, and counts the primes in a window after 10¹⁸: 2398 of them against 2413 expected by the prime number theorem.
In programming contests the test is written exactly like this: trial division by a few small primes, then the bases 2…37 or Sinclair’s set. Two common mistakes: a base that is a multiple of n must be skipped, not counted as a witness (otherwise 2 and 3 will “turn out” composite with base 2 or 3), and in languages without big numbers the product a * b % n overflows for n around 10¹⁸.
What it costs
One round is one exponentiation: ⌊log₂ n⌋ modular squarings and at most as many multiplications, plus at most s squarings of the chain. If numbers of length L = log n are multiplied by the schoolbook method, a multiplication costs O(L²), and a round O(L³). So k rounds are O(k log³ n). With Karatsuba multiplication the exponent drops to 2.585, with the fast Fourier transform almost to two. Trial division up to √n costs √n = 10L/2 divisions: every two new digits make it ten times longer. A measurement: a prime of d digits, trial division against twenty rounds of the test.
Trial division goes straight up: in Python every two digits multiply the time by ten, and a 16-digit prime already takes almost six seconds to check; a 20-digit one would take about ten minutes, a 30-digit one almost two years. Miller–Rabin spends less than a millisecond on a 20-digit number, and on a number of 512 digits half a second in Python and 166 milliseconds in Go. On long numbers its line has a slope between two and three: doubling the length raised the time 5–8 times in Python, and in Go with a large spread, from 3.5 to 11 times. Theory promises a slope of 3, two from schoolbook multiplication and one from the number of multiplications; at such lengths a noticeable share of the time is eaten by overheads that grow more slowly.
Where the test breaks
Fixed bases against someone who chooses the number. A check “to bases 2…41” is reliable for numbers below ψ₁₃ and entirely unreliable for numbers chosen against it: ψ₁₃ is already such a number. Larger ones can be built too. In 1995 François Arnault published a 397-digit Carmichael number that passes the strong test to all prime bases below 307. In 2018 Martin Albrecht, Jake Massimo, Kenneth Paterson and Juraj Somorovsky examined cryptographic libraries and built composite numbers that the libraries took for primes: for OpenSSL, 2048-bit numbers that in the default settings passed the check with probability 1/16 instead of the promised 2−80; for GNU GMP, 1024-bit numbers that always passed with the minimum recommended number of rounds; in Cryptlib, LibTomCrypt, JavaScript Big Number and WolfSSL, numbers that always passed. The danger is not theoretical: an adversary slips a “prime” into Diffie–Hellman parameters, and the discrete logarithm becomes easy. In 2019 OpenSSL added a checking function that cannot be asked for too few rounds: at least 64 for numbers up to 2048 bits and 128 for longer ones; the change cites the research of Massimo and Paterson.
Predictable “random” bases. The bound 4−k holds only if the bases are chosen randomly and independently of n. Bases from a generator with a known seed are the same fixed bases, just harder to see.
The defence. The authors of that work recommend the Baillie–PSW test: one strong round to base 2 plus a strong Lucas test, a check of an entirely different nature. A number that would fool both parts at once has not been found to this day, and its absence has been verified for all numbers below 2⁶⁴. Go’s documentation says plainly that ProbablyPrime is 100% accurate for numbers below 2⁶⁴, but is not suitable for numbers an adversary might have chosen.
Where it works
- Keys. RSA needs two random primes of 1024 bits and more. A program takes a random odd number, sieves it with small primes and checks it with the Miller–Rabin test until it meets a prime; how the sieving works is told in the article on the sieve of Eratosthenes. OpenSSL checks candidates with at least 64 rounds.
- Libraries. In Go
big.Int.ProbablyPrime(n)makes n rounds with pseudo-random bases and the Baillie–PSW test. In GMPmpz_probab_prime_p(n, reps)does trial divisions, Baillie–PSW and another reps − 24 rounds. Python’s standard library has no primality check, but the test is ten lines on the built-inpow(a, d, n). - Factorisation. Together with Pollard’s rho algorithm the test factors numbers: a piece it calls prime is not divided further.
History
Fermat’s little theorem has been known since 1640, and in 1819 Pierre Sarrus showed that it cannot simply be reversed: 2341 ≡ 2 (mod 341), although 341 = 11 · 31. The test with the chain of squares was proposed by Gary Miller in 1976, deterministic but proved only under the extended Riemann hypothesis. In 1980 Michael Rabin turned it into a probabilistic test that needs no hypotheses and proved the bound of a quarter; Louis Monier obtained it at the same time and independently. The same year Pomerance, Selfridge and Wagstaff counted all pseudoprimes up to 25·10⁹ and found the first ψk.
In 2002 Manindra Agrawal, Neeraj Kayal and Nitin Saxena built a deterministic polynomial-time primality test with no hypotheses at all, and the question “can primality be checked quickly” was closed in theory. In practice dice are still thrown: the Miller–Rabin test is incomparably faster, and its error can be made as small as you like.
Problems
In the first problem you have to find all Carmichael numbers in a window far from zero, and figure out how not to factor every number in the window. In the second, count the liars of a number from its factorisation, that is, derive Monier’s formula yourself, at least for one case.
Carmichael numbers in a window How many bases lie