Fermat and Solovay-Strassen primality testing in Scheme

Following up on my implementation of RSA in Scheme, I created an implementation of the probabilistic Fermat and Solovay-Strassen primality tests for the purpose of large prime generation.

These tests are simple and elegant, so I've included a short mathematical explanation of how each works.

The code allows both the byte-size of the prime and the confidence that the result is a prime to be chosen.

This is part of my ongoing effort to create benchmarks for Scheme that exercise only a small part of the language, yet do serious computation.

As both a functional language and a language with support for arbitrary precision integers, Scheme is a beautiful language for describing these algorithms.

The Fermat primality algorithm and the Solovay-Strassen primality algorithm are, by themselves, straightforward.

The Fermat test

The Fermat primality test is based on Fermat's little theorem.

Fermat's little theorem implies that for any prime $p$ and any natural $a$ such that $1 \leq a < p$, it must be the case that

$a^{p-1} \mod p = 1\text.$

However, for any given composite number $n$ and some natural $a$ such that $1 \leq a < n$, it is usually not the case that:

$a^{n-1} \mod n = 1\text.$

Given a number $m$, we can repeatedly test whether $a^{m-1} \mod m = 1$ for many different values of $a$.

Each time we get back 1, we become increasingly confident that the natural $m$ is truly prime.

However, if we ever get back any value other than 1, we know immediately that $m$ is composite.

Warning: Charmichael numbers

It must be noted that there are composites, known as Charmichael numbers, that can pass this test for any value of $a$.

In general, Charmichael numbers are sufficiently rare that this test can be used in practice, with caution.

561 is the lowest Charmichael number.

The Solovay-Strassen test

The Solovay-Strassen test uses a slightly more complex test than the Fermat test, but it provides a tight, probabilistic guarantee on the result.

Euler's criterion implies that for any $a$ coprime with some odd prime $p$:

$a^{(p-1)/2} \mod p = \left( \frac{a}{p} \right)$

[The notation $(\frac{a}{m})$ denotes the Jacobi symbol.]

Meanwhile, for a composite $m$, given an arbitrary natural $a$, there is less than a 50% chance that:

$a^{(m-1)/2} \mod m = \left( \frac{a}{m} \right)$

By successively testing the condition with different values of $a$, the probability that $m$ is composite falls by at least half with each successful test.

Thus, it becomes possible to choose a level of confidence that $m$ is prime and perform the necessary number of tests to achieve that.

The complex part of the algorithm is implementing the Jacobi symbol:

(define (jacobi a n)
(cond
((= n 1) 1)
((= a 1) 1)
((not (= (gcd a n) 1)) 0)
((and (= a 2)
(let ((n-mod-8 (modulo n 8)))
(cond
((or (= n-mod-8 1) (= n-mod-8 7)) 1)
((or (= n-mod-8 3) (= n-mod-8 5)) -1)))))
((> a n) (jacobi (modulo a n) n))
((even? a) (* (jacobi (/ a 2) n) (jacobi 2 n)))
((even? n) (* (jacobi a (/ n 2)) (jacobi a 2)))
(else (* (jacobi n a) (if (even? (/ (* (- a 1) (- n 1)) 4)) 1 -1)))))

Code

Warning: If you want to use this in practice, you'll also have to add a weak-key detection pass to make it reasonably secure.