Introduction to Primality Testing

A prime number is defined as a natural number greater than \(1\) whose only positive divisors are \(1\) and itself. That definition is already an algorithm in disguise: to decide whether \(n\) is prime, search for a divisor \(d\) with \(1<d<n\) using the modulo operator, \(n\bmod d=0\). If such a \(d\) is found, \(n\) is composite; if the search comes up empty, \(n\) is prime. Every optimization below starts from a number-theoretic fact that rules out a whole class of candidates \(d\) as impossible divisors, shrinking the search before a single modulo operation is wasted on it.

Trial Division: Testing Every Candidate

Translating the definition directly gives the most literal possible implementation:

function isPrimeNaive(n) {
    if (n < 2) return false;
    for (let d = 2; d < n; d++) {
        if (n % d === 0) return false;
    }
    return true;
}

This is correct because it checks every integer strictly between \(1\) and \(n\) as a candidate divisor, leaving nothing out. For \(n\) prime, all \(n-2\) candidates are tested and rejected; the running time is \(O(n)\), which is already too slow once \(n\) reaches into the millions.

No Divisor Exceeds \(n/2\)

If \(d\) is a divisor of \(n\) with \(d<n\), then \(n=dq\) for some integer \(q\geq2\) (\(q=1\) would force \(d=n\)). Solving for \(d\),

\[d=\frac nq\leq\frac n2.\]

So no divisor other than \(n\) itself can be larger than \(n/2\), and the loop can stop there:

function isPrimeHalf(n) {
    if (n < 2) return false;
    for (let d = 2; d <= n / 2; d++) {
        if (n % d === 0) return false;
    }
    return true;
}

This roughly halves the work, but the running time is still \(O(n)\) — a constant factor does not change the growth rate. The real improvement comes from looking at how the two factors of a divisor pair relate to each other, not just at how large one of them can individually be.

Bounding the Search by \(\sqrt n\)

If \(n\) is composite, it factors as \(n=ab\) with \(1<a,b<n\). Suppose both \(a>\sqrt n\) and \(b>\sqrt n\) were true at once; then \(ab>\sqrt n\cdot\sqrt n=n\), contradicting \(ab=n\). So at least one of the two factors is at most \(\sqrt n\). Every composite number therefore has a divisor in the range \(2\leq d\leq\sqrt n\), and it suffices to search only that range: if nothing divides \(n\) up to \(\sqrt n\), nothing divides it at all.

function isPrimeSqrt(n) {
    if (n < 2) return false;
    for (let d = 2; d * d <= n; d++) {
        if (n % d === 0) return false;
    }
    return true;
}

The running time drops to \(O(\sqrt n)\): for \(n=10^{12}\), roughly a million candidates instead of half a trillion. In languages with fixed-width integers, d * d can overflow before reaching \(n\); the equivalent, overflow-free condition \(d^2\leq n\iff d\leq n/d\) sidesteps this:

function isPrimeSqrt(n) {
    if (n < 2) return false;
    for (let d = 2; d <= n / d; d++) {
        if (n % d === 0) return false;
    }
    return true;
}

JavaScript's Number represents integers exactly only up to \(2^{53}-1\); everything from here on assumes \(n\) stays within that range unless stated otherwise.

Skipping Even Candidates

Every even number greater than \(2\) is divisible by \(2\) and therefore composite, so \(2\) is the only even prime. Once \(n=2\) has been handled and every other even \(n\) rejected up front, every remaining candidate divisor can be restricted to the odd numbers:

function isPrimeOdd(n) {
    if (n < 2) return false;
    if (n === 2) return true;
    if (n % 2 === 0) return false;
    for (let d = 3; d <= n / d; d += 2) {
        if (n % d === 0) return false;
    }
    return true;
}

This halves the number of candidates again, without changing the \(O(\sqrt n)\) growth rate.

The \(6k\pm1\) Wheel

Every integer leaves one of the remainders \(0,1,2,3,4,5\) modulo \(6\), so every integer has one of the forms \(6k,6k+1,6k+2,6k+3,6k+4,6k+5\). Three of these six classes are eliminated immediately: \(6k\) and \(6k+2=2(3k+1)\) and \(6k+4=2(3k+2)\) are even, and \(6k+3=3(2k+1)\) is a multiple of \(3\). Only \(6k+1\) and \(6k+5=6(k+1)-1\) survive, so

\[p>3\text{ prime}\ \Longrightarrow\ p=6k+1\ \text{or}\ p=6k-1\ \text{for some }k.\]

The converse fails — \(25=6\cdot4+1\) is composite — so this is a necessary, not sufficient, condition; it discards candidates rather than confirming primes. After removing multiples of \(2\) and \(3\) up front, only two candidates per block of six remain to be tested:

function isPrime6k(n) {
    if (n < 2) return false;
    if (n <= 3) return true;
    if (n % 2 === 0 || n % 3 === 0) return false;
    for (let d = 5; d <= n / d; d += 6) {
        if (n % d === 0 || n % (d + 2) === 0) return false;
    }
    return true;
}

The candidate density has dropped to \(2/6=1/3\), down from \(1/2\) for the odd-only version.

Larger Wheels

The \(6k\pm1\) pattern is the special case \(W=2\cdot3=6\) of a general idea: pick a product \(W\) of small primes, and only test candidates \(d\) with \(\gcd(d,W)=1\), since any \(d\) sharing a factor with \(W\) is already divisible by one of those small primes. The count of such candidates in each block of \(W\) consecutive integers is Euler's totient \(\varphi(W)\); for \(W=30=2\cdot3\cdot5\),

\[ \varphi(30)=30\left(1-\frac12\right)\left(1-\frac13\right)\left(1-\frac15\right)=8, \]

and the coprime residues in \([1,30]\) are \(1,7,11,13,17,19,23,29\), giving a candidate density of \(8/30\approx26.7\%\) versus \(1/3\) for \(W=6\). The gains keep shrinking as \(W\) grows, while the bookkeeping for skipping between residues gets more complicated, so wheels beyond \(W=30\) or \(W=210\) are rarely worth the added complexity for a single primality check.

Testing Only Prime Divisors

A wheel discards candidates that are obviously composite by construction; it does not discard candidates that merely happen to be composite. But as shown in the introduction to prime numbers, every divisor \(d\) that is itself composite has a smaller prime divisor \(p\mid d\), and \(p\mid d\mid n\) then gives \(p\mid n\) as well. So a composite candidate never needs to be tested directly — one of its prime factors will already have been tried and found earlier in the search. It is therefore enough to test only prime candidates, given a list of them up to \(\sqrt n\):

function isPrimeWithTable(n, primes) {
    if (n < 2) return false;
    for (const p of primes) {
        if (p > n / p) break;
        if (n % p === 0) return n === p;
    }
    return true;
}

This is only correct if primes contains every prime up to \(\sqrt n\). Since the count of primes below \(x\) grows like \(x/\ln x\), a single check now costs roughly \(\sqrt n/\ln\sqrt n\) divisions instead of \(\sqrt n/3\) — a genuine asymptotic improvement, but only worthwhile once the list of small primes is built once and reused across many checks, which raises the question of how to build that list efficiently.

Generating Many Primes at Once: The Sieve of Eratosthenes

A primality test answers one question, whether one specific \(n\) is prime. The sieve of Eratosthenes answers a different one: which numbers up to a bound \(N\) are prime. Instead of testing each candidate independently, it marks composites by their multiples: starting from the smallest unmarked number, every multiple of it is composite and gets marked, and the next unmarked number is prime.

function sieve(N) {
    const composite = new Uint8Array(N + 1);
    composite[0] = composite[1] = 1;
    for (let p = 2; p <= N / p; p++) {
        if (composite[p]) continue;
        for (let m = p * p; m <= N; m += p) {
            composite[m] = 1;
        }
    }
    return composite;
}

Two design choices here follow directly from facts already established above. The inner loop can start at \(p^2\) rather than \(2p\), because every smaller multiple \(kp\) with \(k<p\) has already been marked when \(k\) or one of its prime factors was processed as the outer loop variable. And the outer loop can stop once \(p>\sqrt N\), by the same argument used for a single \(\sqrt n\) trial-division test: every composite \(m\leq N\) has a prime factor at most \(\sqrt m\leq\sqrt N\), so once every prime up to \(\sqrt N\) has swept through its multiples, every composite up to \(N\) is already marked. The if (composite[p]) continue; guard is not needed for correctness — marking multiples of a composite number only remarks numbers that are already composite — but it is essential for performance, since without it the same numbers get marked repeatedly instead of only once per prime factor.

The running time is \(O(N\log\log N)\), because the inner loop for each prime \(p\) runs \(N/p\) times, and the sum \(\sum_{p\leq N}1/p\) over the primes grows like \(\log\log N\).

Fermat's Primality Test

All the improvements so far still fundamentally search for a divisor. Fermat's little theorem offers a different kind of evidence, one that says nothing about factors: for a prime \(p\) and any integer \(a\) not divisible by \(p\),

\[a^{p-1}\equiv1\pmod p.\]

Contrapositively, if \(a^{n-1}\not\equiv1\pmod n\) for some \(a\), then \(n\) cannot be prime. This gives a test: pick a base \(a\), and compute \(a^{n-1}\bmod n\) using the modpow routine from modular exponentiation, which computes it in \(O(\log n)\) operations without ever forming the full power \(a^{n-1}\):

function passesFermatTest(n, a) {
    if (n < 2) return false;
    return modpow(a, n - 1, n) === 1;
}

The converse of Fermat's little theorem is false: \(a^{n-1}\equiv1\pmod n\) does not prove that \(n\) is prime, only that \(n\) survived this one test with this one base. A composite number that passes is called a Fermat pseudoprime to that base; for example \(341=11\cdot31\) satisfies \(2^{340}\equiv1\pmod{341}\) despite being composite. Worse, a Carmichael number such as \(561=3\cdot11\cdot17\) passes the test for every base coprime to it, so no single base, and no finite hard-coded list of exceptions, can fix the test — infinitely many Carmichael numbers exist. Fermat's test alone can rule primality out with certainty, but it can never rule it in.

The Miller–Rabin Test

Miller–Rabin strengthens Fermat's test by using one more fact about primes: for a prime \(p\), the only square roots of \(1\) modulo \(p\) are \(1\) and \(-1\). Indeed, \(x^2\equiv1\pmod p\) means \(p\mid(x-1)(x+1)\), and since \(p\) is prime, the fact that a prime dividing a product divides one of its factors forces \(p\mid(x-1)\) or \(p\mid(x+1)\), i.e.\ \(x\equiv1\) or \(x\equiv-1\pmod p\). So if repeatedly squaring some value modulo a prime ever lands exactly on \(1\), the step immediately before it must have been \(-1\), never anything else.

Write \(n-1=2^sd\) with \(d\) odd, by pulling out every factor of \(2\). Fermat's theorem gives \(a^{n-1}\equiv1\pmod n\) for a prime \(n\), and this value is reached by repeatedly squaring \(a^d\) exactly \(s\) times:

\[a^d,\ a^{2d},\ a^{4d},\ \ldots,\ a^{2^sd}=a^{n-1}\pmod n.\]

For \(n\) prime, this sequence must end at \(1\); by the square-root argument above, the entry immediately before the first \(1\) can only be \(-1\), unless the sequence starts at \(1\) outright. A base \(a\) is called a witness against \(n\) if neither of these patterns occurs — proof that \(n\) is composite. A single round tests one base:

function millerRabinRound(n, a, d, s) {
    let x = modpow(a, d, n);
    if (x === 1 || x === n - 1) return true;
    for (let r = 1; r < s; r++) {
        x = (x * x) % n;
        if (x === n - 1) return true;
        if (x === 1) return false;
    }
    return false;
}

function isPrimeMillerRabin(n, rounds = 20) {
    if (n < 2) return false;
    if (n % 2 === 0) return n === 2;
    let d = n - 1, s = 0;
    while (d % 2 === 0) { d /= 2; s++; }
    for (let i = 0; i < rounds; i++) {
        const a = 2 + Math.floor(Math.random() * (n - 3));
        if (!millerRabinRound(n, a, d, s)) return false;
    }
    return true;
}

Unlike Fermat pseudoprimes, no number is a Miller–Rabin liar for more than a quarter of the possible bases: for any odd composite \(n\), at most \(1/4\) of the bases in \(1<a<n-1\) fail to expose it. After \(k\) independent random rounds, the probability that a composite number is mistakenly accepted is therefore at most \((1/4)^k\), which is already below \(2^{-32}\) at \(k=16\) and below \(2^{-64}\) at \(k=32\) — far smaller than the odds of a hardware fault during the computation itself.

A Deterministic Test for 64-Bit Integers

Randomized correctness bounds are only necessary because arbitrary bases are being tried against arbitrary \(n\). Restricted to a bounded range of \(n\), a small, fixed set of bases can be verified in advance to expose every composite in that range, turning Miller–Rabin into a deterministic test with no probability of error at all. For every \(n<2^{64}\), testing only the seven bases

\[2,\ 325,\ 9375,\ 28178,\ 450775,\ 9780504,\ 1795265022\]

is known to correctly classify every integer in that range. Since \(2^{64}\) is far beyond the \(2^{53}-1\) limit of an exact Number, this final version switches to BigInt for exact arithmetic:

function modpowBig(base, exp, mod) {
    base %= mod;
    let result = 1n;
    while (exp > 0n) {
        if (exp & 1n) result = (result * base) % mod;
        base = (base * base) % mod;
        exp >>= 1n;
    }
    return result;
}

function millerRabinRoundBig(n, a, d, s) {
    let x = modpowBig(a, d, n);
    if (x === 1n || x === n - 1n) return true;
    for (let r = 1; r < s; r++) {
        x = (x * x) % n;
        if (x === n - 1n) return true;
        if (x === 1n) return false;
    }
    return false;
}

function isPrime64(n) {
    n = BigInt(n);
    if (n < 2n) return false;
    for (const p of [2n, 3n, 5n, 7n, 11n, 13n]) {
        if (n === p) return true;
        if (n % p === 0n) return false;
    }
    let d = n - 1n, s = 0;
    while (d % 2n === 0n) { d /= 2n; s++; }
    const bases = [2n, 325n, 9375n, 28178n, 450775n, 9780504n, 1795265022n];
    for (const a of bases) {
        const base = a % n;
        if (base === 0n) continue;
        if (!millerRabinRoundBig(n, base, d, s)) return false;
    }
    return true;
}

Choosing an Approach

Every step in this chapter eliminated a class of divisors that could not possibly work, whether by an inequality on the size of factors, a residue argument modulo a small number, or the requirement that a divisor's own prime factors be tested first. Miller–Rabin is the point where this search for a witness of compositeness turns from checking specific candidate divisors into checking an algebraic identity that every prime is guaranteed to satisfy — the same shift in strategy, from exhaustive search to algebraic certificate, that also underlies fast GCD computation and modular exponentiation itself.