The first algorithm most programmers meet is the naïve trial-division method. After skipping the even numbers and stopping at √n, the code looks like this:

(define (naive-prime? k)
  (cond [(= k 2) #t]
        [(even? k) #f]
        [else (trial-loop k 3 (ceiling (sqrt k)))]))

(define (trial-loop n d limit)
  (cond [(> d limit) #t]
        [(zero? (modulo n d)) #f]
        [else (trial-loop n (+ d 2) limit)]))

This is already an improvement over testing every integer up to n, but its Θ(√n) complexity still explodes quickly.

Incremental sieve

A slightly bettter idea is to keep a list of primes discovered so far and only divide the next candidate by those primes. The Scheme version below grows the list until it reaches the requested bound:

(define (primes-up-to bound)
  (let next ((c 3) (found '(2)))
    (cond [(> c bound) (reverse found)]
          [(for-all (lambda (p) (not (zero? (modulo c p))))
                    found)
           (next (+ c 2) (cons c found))]
          [else (next (+ c 2) found)])))

Because every new candidate is tested against all previously found primes, the complexity is roughly Θ(n2 / log n), which is still too slow for large inputs.

Fermat’s little theorem

The Fermat test replaces division with modular exponentiation. If p is prime, then for every 1 ≤ a < p we have ap ≡ a (mod p). The converse is not true: some composite (Carmichael numbers) satisfy the congruence for every a coprime to them.

(define (fermat-witness? n a)
  (not (= (expmod a n n) a)))

(define (fermat-probable-prime? n rounds)
  (let test ((k rounds))
    (or (zero? k)
        (and (fermat-witness? n (+ 1 (random (- n 1))))
             (test (- k 1))))))

Unfortunately, Carmichael numbers such as 561, 1105, 1729, … always pass this test, so the Fermat method is unreliable for cryptographic purposes.

Miller–Rabin refinement

The Miller–Rabin test fixes the Carmichael problem by exploiting the fact that any non-trivial square root of 1 modulo n proves n composite. The algorithm writes n − 1 as d·2r and looks for a witness a that violates the required sequence of squares.

(define (miller-rabin-witness? n a)
  ( (let* ((n-1 (- n 1))
          (r (factor-out-twos n-1))
          (d (/ n-1 (expt 2 r)))
          (x (expmod a d n)))
    (or (and (not (= x 1)) (not (= x n-1)))
        (let loop ((i 1) (x x))
          (and (< i r)
               (let ((next (modulo (* x x) n)))
                 (or (and (= next 1) (not (= x 1)) (not (= x n-1)))
                     (loop (+ i 1) next))))))))

(define (factor-out-twos m)
  (if (even? m)
      (+ 1 (factor-out-twos (/ m 2)))
      0))

(define (miller-rabin-prime? n trials)
  (or (= n 2)
      (and (odd? n)
           (let test ((t trials))
             (or (zero? t)
                 (and (not (miller-rabin-witness? n (+ 1 (random (- n 1)))))
                      (test (- t 1))))))))

For k random bases the error probability is at most 4−k. Ten iterations already drive the failure rate below the chance of a cosmic-ray bit-flip.

Quick sanity check

The following snippet verifies that every Carmichael number below 106 is correctly rejected by the Miller–Rabin test with five random bases:

(define carmichaels
  '(561 1105 1729 2465 2821 6601 8911 10585 15841 29341
    41041 46657 52633 62745 63973 75361 101101 115921
    126217 162401 172081 188461 252601 278545 294409))

(for-each
 (lambda (c)
   (when (miller-rabin-prime? c 5)
     (display c) (display " incorrectly accepted\n")))
 carmichaels)

No output means every Carmichael number was properly flagged composite.