From Mathematics to Generic Programming (2011)

---

13.2 Primality Testing

The problem of distinguishing prime numbers from composite ... is known to be one of the most important and useful in arithmetic.

C. F. Gauss, Disquisitiones Arithmeticae

An important problem in modern cryptography is determining whether an integer is prime. Gauss believed that (1) deciding whether a number is prime or composite is a very hard problem, and so is (2) factoring a number. He was wrong about #1, as we shall see. So far, he seems to be right about #2, which is a good thing for us, since modern cryptosystems are based on this assumption.

To find out if a number n is prime, it helps to have a predicate that tells us whether it’s divisible by a given number i:

Click here to view code image

template <Integer I>
bool divides(const I& i, const I& n) {
return n % i == I(0);
}

We can call this repeatedly to find the smallest divisor of a given number n. Just as we did with the Sieve of Eratosthenes in Chapter 3, our loop for testing divisors will start at 3, advance by 2, and stop when the square of the current candidate reaches n:

Click here to view code image

template <Integer I>
I smallest_divisor(I n) {
// precondition: n > 0
if (even(n)) return I(2);
for (I i(3); i * i <= n; i += I(2)) {
if (divides(i, n)) return i;
}
return n;
}

Now we can create a simple function to determine whether n is prime:

Click here to view code image

template <Integer I>
I is_prime(const I& n) {
return n > I(1) && smallest_divisor(n) == n;
}

This is mathematically correct, but it’s not going to be fast enough. Its complexity is . That is, it’s exponential in the number of digits. If we want to test a 200-digit number, we may be waiting for more time than the life of the universe.

To overcome this problem, we’re going to need a different approach, which will rely on the ability to do modular multiplication. We’ll use a function object that provides what we need:

Click here to view code image

template <Integer I>
struct modulo_multiply {
I modulus;
modulo_multiply(const I& i) : modulus(i) {}

I operator() (const I& n, const I& m) const {
return (n * m) % modulus;
}
};

We’ll also need an identity element:

Click here to view code image

template <Integer I>
I identity_element(const modulo_multiply<I>&) {
return I(1);
}

Now we can compute a multiplicative inverse modulo prime p. It uses the result we showed in Chapter 5 that, as a consequence of Fermat’s Little Theorem, the inverse of an integer a, where 0 < a < p, is a^p−2 (see Section 5.4, right after the proof of Fermat’s Little Theorem). It also uses the power function we created in Chapter 7:

Click here to view code image

template <Integer I>
I multiplicative_inverse_fermat(I a, I p) {
// precondition: p is prime & a > 0
return power_monoid(a, p - 2, modulo_multiply<I>(p));
}

With these pieces, we can now use Fermat’s Little Theorem to test if a number n is prime. Recall that Fermat’s Little Theorem says:

If p is prime, then a^p−1 − 1 is divisible by p for any 0 < a < p.

Equivalently:

If p is prime, then a^p–1 = 1 mod p for any 0 < a < p.

We want to know if n is prime. So we take an arbitrary number a smaller than n, raise it to the n − 1 power using modular multiplication (mod n), and check if the result is 1. (We call the number a we’re using a witness.) If the result is not equal to 1, we know definitely by the contrapositive of the theorem that n is not prime. If the result is equal to 1, we know that there’s a good chance that n is prime, and if we do this for lots of random witnesses, there’s a very good chance that it is prime:

Click here to view code image

template <Integer I>
bool fermat_test(I n, I witness) {
// precondition: 0 < witness < n
I remainder(power_semigroup(witness,
n - I(1),
modulo_multiply<I>(n)));
return remainder == I(1);
}

This time we use power_semigroup instead of power_monoid, because we know we’re not going to be raising anything to the power 0. The Fermat test is very fast, because we have a fast way to raise a number to a power—our O(log n) generalized Egyptian multiplication algorithm from Chapter 7.

* * *

While the Fermat test works the vast majority of the time, it turns out that there are some pathological cases of numbers that will fool it for all witnesses coprime to n; they produce remainders of 1 even though they’re composite. These are called Carmichael numbers.

Definition 13.1. A composite number n > 1 is a Carmichael number if and only if

∀b > 1, coprime(b, n) bⁿ⁻¹ = 1 mod n

172081 is an example of a Carmichael number. Its prime factorization is 7 · 13 · 31 · 61.

Exercise 13.1. Implement the function:

bool is_carmichael(n)

Exercise 13.2. Find the first seven Carmichael numbers using your function from Exercise 13.1.

13.3 The Miller-Rabin Test

To avoid worrying about Carmichael numbers, we’re going to use an improved version of our primality tester, called the Miller-Rabin test; it will again rely on the speed of our power algorithm.

We know that n − 1 is even (it would be awfully silly to run a primality test on an even n), so we can represent n − 1 as the product 2^k · q. The Miller-Rabin test uses a sequence of squares w^20q, w^21q, ..., w^2kq, where w is a random number less than the one we are testing. The last exponent in this sequence is n − 1, the same value the Fermat test uses; we’ll see why this is important shortly.

We’re also going to rely on the self-canceling law (Lemma 5.3), except we’ll write it with new variable names and assuming modular multiplication:

For any 0 < x < n ∧ prime(n), x² = 1 mod n x = 1 ∨ x = −1

Remember that in modular arithmetic, −1 mod n is the same as (n − 1) mod n, a fact that we rely on in the following code. If we find some x² = 1 mod n where x is neither 1 nor −1, then n is not prime.

Now we can make two observations: (1) If x² = 1 mod n, then there’s no point in squaring x again, because the result won’t change; if we reach 1, we’re done. (2) If x² = 1 mod n and x is not −1, then we know n is not prime (since we already ruled out x = 1 earlier in the code).

Here’s the code, which returns true if n is probably prime, and false if it definitely is not:

Click here to view code image

template <Integer I>
bool miller_rabin_test(I n, I q, I k, I w) {

// precondition n > 1 ∧ n – 1 = 2^kq ∧ q is odd

modulo_multiply<I> mmult(n);
I x = power_semigroup(w, q, mmult);
if (x == I(1) || x == n - I(1)) return true;
for (I i(1); i < k; ++i) {

// invariant x = w^{2i – 1 q}

x = mmult(x, x);
if (x == n - I(1)) return true;
if (x == I(1)) return false;
}
return false;
}

Note that we pass in q and k as arguments. Since we’re going to call the function many times with different witnesses, we don’t want to refactor n – 1 every time.

Why can we return true at the beginning if the power_semigroup call returns 1 or – 1? Because we know that squaring the result will give 1, and squaring is equivalent to multiplying the exponent in the power calculation by a factor of 2, and doing this k times will make the exponent n – 1, the value we need for Fermat’s Little Theorem to hold. In other words, if w^q mod n = 1 or –1, then w^2k q mod n = w^{n – 1} mod n = 1.

Let’s look at an example. Suppose we want to know if n = 2793 is prime. We choose a random witness w = 150. We factor n – 1 = 2792 into 2² · 349, so q = 349 and k = 2. We compute

x = w^q mod n = 150³⁴⁹ mod 2793 = 2019

Since the result is neither 1 nor – 1, we start squaring x:

i = 1; x² = 150²¹ ·349 mod 2793 = 1374

i = 2; x² = 150²² ·349 mod 2793 = 2601

Since we haven’t reached 1 or – 1 yet, and i = k, we can stop and return false; 2793 is not prime.

Like the Fermat test, the Miller-Rabin test is right most of the time. Unlike the Fermat test, the Miller-Rabin test has a provable guarantee: it is right at least 75% of the time² for a random witness w. (In practice, it’s even more often.) Randomly choosing, say, 100 witnesses makes the probability of error less than 1 in 2²⁰⁰. As Knuth remarked, “It is much more likely that our computer has dropped a bit, due ... to cosmic radiations.”

² In fact, the requirement that q must be odd is needed for this probability guarantee.