# Finding large pseudoprimes

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

bp−1 = 1 (mod p).

This gives a necessary but not sufficient test for a number to be prime. A number that satisfies the equation above but is not prime is called a pseudoprime [1] for base b. For example, 341 is a pseudoprime base 2 because

2340 = 1 mod 341

and clearly 341 = 11*31.

How might you go about hunting for pseudoprimes? This page gives a way to find pseduoprimes base 2, quoting the following theorem.

If both numbers q and 2q − 1 are primes and n = q(2q − 1) then 2n−1 = 1 (mod n) if and only if q is of the form 12k + 1.

Here’s Python code that will find pseudoprimes of a given size in bits.

```    from secrets import randbits
from sympy import isprime

def find_pseudoprime(numbits):
n = numbits // 2
while True:
q = randbits(n)
if isprime(q) and q%12 == 1 and isprime(2*q-1):
return q*(2*q-1)

print( find_pseudoprime(256) )
```

This returned

47362882341088596725068562696893704769436677460225591859092704246296157080253

in under a second.

Clearly the return value of `find_pseudoprime` is composite since it is computed as the product of two integers. You could confirm it returns a pseudoprime by computing

```    pow(2, p-1, p)
```

on the return value to verify that it equals 1.

***

[1] Sometimes such a number is called a Fermat pseudoprime to distinguish it from other kinds of pseudo primes, i.e. numbers that satisfy other necessary conditions for being prime without actually being prime. For example, Perrin pseudoprimes. See the next post for more examples.