Finding large pseudoprimes

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

b^p-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

2³⁴⁰ = 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 2^n-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.