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  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) )
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.
 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.