I was reading a book on computer algebra and ran across an interesting theorem about Carmichael numbers in the one of the exercises. I’ll present that theorem below, but first I’ll back up and say what a pseudoprime is and what a Carmichael number is.

## Fermat’s theorem

If *p* is a prime number, then for any integer *b*, Fermat’s little theorem says

*b*^{p−1} = 1 (mod *p*).

The contrapositive of this theorem is useful for testing whether a number is (not) prime: if for some *b*,

*b*^{p−1} ≠ 1 (mod *p*)

then *p* is not a prime. For example, we could test whether 1003 is prime by computing [1]

2^{1002} (mod 1003).

When we do, we find the result is 990, not 1, and so 1003 is not prime.

## Pseudoprimes

Let’s test whether 341 is prime. If we compute

2^{340} (mod 341)

we get 1. But clearly 341 = 11 × 31. If we go back and read the top of the post carefully, we see that Fermat’s theorem only gives a way to prove a number is **not** prime. We say 341 is a pseudoprime base 2 because it is a composite number that slips past the primality test based on Fermat’s little theorem with *b* = 2.

Strictly speaking we should say 341 is a *Fermat* pseudoprime base 2. When someone says “pseudoprime” without further qualification, they usually mean Fermat pseudo prime. But there are other kinds of pseudoprimes. For any theorem that gives a necessary but not sufficient for a number to be prime, the composite numbers that satisfy the condition are called pseudoprimes. See this post for examples, such as Euler pseudoprimes and Wilson pseudoprimes.

Note that 341 is not a Fermat pseudoprime base 3 because

3^{340} = 56 (mod 341)

If a number passes Fermat’s primality test for several bases, it’s probably prime. (Here’s a blog post that carefully unpacks what “probably prime” means.)

## Carmichael numbers

We saw that 341 is a Fermat pseudoprime for base 2 but not base 3; it slipped past our first check on primality but not our second. Some numbers are more slippery. Carmichael numbers are composite numbers that slip past Fermat’s primality test for every base *b*.

Actually, not *every* base. A Carmichael number *n* passes Fermat’s primality test for any base *b* relatively prime to *n*.

The smallest Carmichael number is 561. There are infinitely many Carmichael numbers, but they’re distributed very thinly.

## Searching for Carmichael numbers

Now I’m finally ready to get to the exercise I mentioned at the top of the post. I’m reading *Computer Algebra* by Wolfram Koepf. The exercise asks the reader to prove that if for a given *k* the three numbers 6*k *+ 1, 12*k *+ 1, and 18*k *+ 1 are all prime, then their product is a Carmichael number. This theorem was discovered by J. Chernick in 1939.

We can use this to search for Carmichael numbers with the following Python code. (Note that this approach doesn’t produce all Carmichael numbers—it doesn’t produce 561, for example.) It is unknown whether it produces an infinite sequence of Carmichael numbers.

from sympy import isprime for k in range(2, 100): if isprime(6*k+1) and isprime(12*k+1) and isprime(18*k+1): print(k, (6*k+1)*(12*k+1)*(18*k+1))

Testing values of *k* up to 100 reveals six Carmichael numbers.

6 294409 35 56052361 45 118901521 51 172947529 55 216821881 56 228842209

## Related posts

- Testing for primes less than a quintillion
- Searching for Fermat primes
- How likely is it that a probable prime is prime?

[1] This may not look like a good idea at first. 2^{1002} is an enormous number. Surely it would be easier to test whether 1003 is prime than to compute such a large number. Except you don’t have to compute 2^{1002} per se: you have to compute 2^{1002} (mod 1003). You can reduce your intermediate results (mod 1003) at every step, so you never have to work with large numbers. You can also use fast exponentiation, and so Fermat’s primality test is much more efficient than it may seem at first glance.

Python’s `pow`

function optionally takes a modulus as a third argument. So, for example,

pow(2, 1002)

computes 2^{1002}, and

pow(2, 1002, 1003)

computes 2^{1002} (mod 1003). It produces the same result as

2**1002 % 1003

but is more efficient. Not that it matters in this example, but it matters more when the exponent is much larger.