Fermat’s little theorem says that if *p* is a prime number, then for any positive integer *b* < *p* we hve

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

This theorem gives a necessary but not sufficient condition for a number to be prime.

## Fermat’s primality test

The converse of Fermat’s little theorem is not always true, but it’s often true. That is, if there exists some base 1 < *b* < *n* such that

*b*^{n−1} = 1 (mod *n*)

then *n* is likely to be prime. There are examples where the equation above holds for a pair (*b*, *n*) even though *n* is not prime, and in that case *n* is called a **pseudoprime** to the base *b*.

If you’re searching for large primes, say for use in encryption, then you’d begin by applying Fermat’s little theorem with a few small values of *b*. This is because although Fermat’s test can’t prove that a number *is* prime, it can prove that a number *is not* prime.

For a small example, suppose you wanted to test whether 50621 is prime. You could start by applying Fermat’s test with *b* = 2 as in the following Python code.

>>> *n* = 50621

>>> 2**(*n*−1) % *n*

9605

Since the result is not 1, we know 50621 is not prime. This doesn’t tell us what the factors of 50621 are, but we know that it *has* nontrivial factors. We say 2 is a *witness* that the number 50621 is not prime.

Next, let’s see whether 294409 might be prime.

>>> *n* = 294409

>>> 2**(*n – *1) % *n*

1

This tells us 294409 *might* be prime. It has passed a test that filters out a lot of composite numbers. What now? We could try other values of *b*: 3, 5, 7, 11, …. This will not resolve the question of whether 294409 is prime unless we keep going until we try 37. And in fact 37 is the smallest factor of 294409. Our number 294409 is a **Carmichael number**, a composite number *n* that passes Fermat’s primality test for all bases *b* relatively prime to *n*.

Note that it would be more efficient to use pow(*b*, *n *− 1, *n*) rather than 2**(*n *− 1) % *n* because the former takes advantage of the fact that we don’t need to compute 2^{n−1} per se and can reduce all intermediate calculations mod *n*.

## Factoring pseudoprimes

Now suppose we have a number *n* that has passed Fermat’s primality test for some base *b* and we suspect that *n* is a pseudoprime. If we want to (try to) factor *n*, knowing that it is a pseudoprime to the base *b* gives us a **head start**. We can exploit the fact that we know *b* to factor *n* in polynomial time, unless *n* is a **strong pseudoprime**.

Suppose we have a number *n* that we suspect is a pseudoprime to the base *b*, and we’re smart enough to at least check that *n* is an odd number, then we begin by pulling out all the factors of 2 that we can from *n* − 1*:*

*n* − 1 = 2^{e} *f*.

Next consider the set of numbers

*b*^{kf}

for *k* = 1, 2, 4, …, 2^{e}. Let *x* be the smallest of these numbers which is not congruent to 1 mod *n*. The existence of such an *x* is essentially the definition of strong pseudoprime [1].

Then gcd(*x* − 1, *n*) and gcd(*x* + 1, *n*) are factors of *n*. This is theorem 10.4 of [2].

## Python example

Let *n* = 873181. This is a pseudoprime to the base *b* = 3, which we can confirm by seeing that pow(3, *n*−1, *n*) returns 1.

Now 873180 is divisible by 4 but not by 8, so *e* = 2. So the theorem above says we should compute

>>> *b*, e = 3, 2

>>> [pow(*b*, *f**2***k*, *n*) for *k* in range(*e*+1)]

This produces [2643, 1, 1]. So *x* = 2643,

>> *x* = 2643

>>> from sympy import gcd

>>> gcd(*x*−1, *n*)

1321

>>> gcd(*x*+1, *n*)

661

shows that 1321 and 661 are both factors of 873181.

## Related posts

[1] Definition of strong pseudoprime. A strong pseudo prime to base *b* is a composite odd integer *m* such that if *m* − 1 = 2^{e}*f* with *f* odd, then either *b*^{f} = 1 (mod *m*) or *b*^{f2c} ≡ −1 (mod *m*) for some 0 ≤ *c* < *e*.

[2] The Joy of Factoring by Samuel S. Wagstaff, Jr.