Miller-Rabin Primality Test: Probabilistic Prime Detection Explained

Prime numbers are the building blocks of number theory and modern cryptography. From RSA encryption to hash functions, the ability to quickly determine whether a large number is prime is crucial. But for numbers with hundreds of digits, deterministic primality tests like trial division or the AKS algorithm are either too slow or too complex. Enter the Miller-Rabin primality test: a probabilistic algorithm that trades a tiny chance of error for blazing speed.

Developed by Gary L. Miller and Michael O. Rabin, this test is the workhorse of primality testing in practice. It works by checking whether the number satisfies certain properties of prime numbers. If it fails any test, it's definitely composite. If it passes many tests, it's 'probably prime' with an error probability that can be made arbitrarily small.

In this tutorial, you'll learn the mathematics behind Miller-Rabin, how to implement it in Python, and how to use it in production systems. We'll also cover common pitfalls, debugging strategies, and interview questions. By the end, you'll be able to confidently apply this algorithm to real-world problems.

Mathematical Foundation

The Miller-Rabin test is based on two properties of prime numbers. Let n be an odd prime, and write n-1 = 2^s d, where d is odd. Then for any integer a not divisible by n, either a^d ≡ 1 (mod n) or a^(2^r d) ≡ -1 (mod n) for some 0 ≤ r < s. This is a consequence of Fermat's little theorem and the fact that if x^2 ≡ 1 (mod p) then x ≡ ±1 (mod p) for prime p.

If n is composite, there exists at least one a (called a witness) for which neither condition holds. The algorithm picks random a's and checks these conditions. If any a fails, n is composite. If all pass, n is probably prime.

The probability that a composite number passes a single test is at most 1/4. By repeating k times, the error probability drops to 4^(-k). For cryptographic applications, k=20 is common, giving an error probability of 4^(-20) ≈ 10^(-12).

def miller_rabin_test(n, a):
    # Write n-1 as 2^s * d
    s = 0
    d = n - 1
    while d % 2 == 0:
        s += 1
        d //= 2
    x = pow(a, d, n)
    if x == 1 or x == n - 1:
        return True  # probably prime
    for _ in range(s - 1):
        x = pow(x, 2, n)
        if x == n - 1:
            return True
    return False  # composite

Full Algorithm Implementation

The full Miller-Rabin algorithm repeats the test with multiple bases. For random bases, we choose k random integers a in [2, n-2]. For deterministic variants, we use a fixed set of bases that are known to work for numbers up to a certain bound.

Here's a complete implementation with configurable iterations:

import random

def is_prime_miller_rabin(n, k=10):
    if n < 2:
        return False
    if n == 2 or n == 3:
        return True
    if n % 2 == 0:
        return False
    # Write n-1 as 2^s * d
    s = 0
    d = n - 1
    while d % 2 == 0:
        s += 1
        d //= 2
    # Witness loop
    for _ in range(k):
        a = random.randrange(2, n - 1)
        x = pow(a, d, n)
        if x == 1 or x == n - 1:
            continue
        for _ in range(s - 1):
            x = pow(x, 2, n)
            if x == n - 1:
                break
        else:
            return False  # composite
    return True  # probably prime

Optimization Techniques

Performance is critical when testing many large numbers. Key optimizations include: - Precompute s and d once per n. - Use early exit for small n or even numbers. - For repeated tests on the same n, cache results. - Use bit operations for division by 2. - In Python, the built-in pow() is highly optimized in C.

For bulk primality testing (e.g., generating many primes), consider using a sieve to eliminate obvious composites first.

def is_prime_optimized(n, bases):
    if n < 2:
        return False
    if n in (2, 3):
        return True
    if n % 2 == 0:
        return False
    s = (n - 1) & - (n - 1)  # lowest set bit
    d = (n - 1) // s
    s = s.bit_length() - 1
    for a in bases:
        if a % n == 0:
            continue
        x = pow(a, d, n)
        if x == 1 or x == n - 1:
            continue
        for _ in range(s - 1):
            x = pow(x, 2, n)
            if x == n - 1:
                break
        else:
            return False
    return True

Error Probability and Iteration Count

The error probability after k iterations is at most 4^(-k). However, this bound is worst-case; for most composites, the probability is much lower. In practice, k=10 gives an error probability of about 10^(-6), which is acceptable for many non-cryptographic applications. For cryptography, k=20 or more is standard.

Note that the test only errs in one direction: it may call a composite prime, but never calls a prime composite. This makes it suitable for applications where false positives are tolerable but false negatives are not.

# Example: compute error probability for k iterations
def error_prob(k):
    return 4**(-k)

for k in [1, 5, 10, 20, 40]:
    print(f"k={k}: error <= {error_prob(k):.2e}")

Deterministic Variants

For numbers up to certain bounds, we can use a fixed set of bases that make the test deterministic. These sets have been proven to work for all numbers below the bound. For example: - n < 2^32: bases [2, 7, 61] - n < 2^64: bases [2, 325, 9375, 28178, 450775, 9780504, 1795265022] (or the smaller set [2, 3, 5, 7, 11, 13, 17]) - n < 3,317,044,064,679,887,385,961,981: bases [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]

Using these sets eliminates randomness and guarantees correctness.

def is_prime_deterministic(n):
    if n < 2:
        return False
    small_primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]
    for p in small_primes:
        if n % p == 0:
            return n == p
    # Use deterministic bases for n < 2^64
    bases = [2, 325, 9375, 28178, 450775, 9780504, 1795265022]
    return is_prime_optimized(n, bases)

Real-World Applications

In each case, the speed of Miller-Rabin makes it the go-to choice over deterministic tests like AKS.

import random

def generate_large_prime(bits):
    while True:
        n = random.getrandbits(bits)
        # Ensure n is odd and has the high bit set
        n |= (1 << (bits - 1)) | 1
        if is_prime_miller_rabin(n, k=20):
            return n