Pseudoprimes To Bases 2 And 3
December 10, 2019
Regular readers know of my affinity for the Online Encyclopedia of Integer Sequences. I was browsing the other day and came across A052155, the sequence of pseudoprimes to both base2 and base3. This sequence is important because it forms the basis of a very fast primality test for numbers in the range 232 to 264, which is a very important range for some large factoring algorithms (particularly SIQS and NFS in their 2LP variants): a number n in range is prime if 2n-1 ≡ 1 (mod n), 3n-1 ≡ 1 (mod n), and n is not in a small list of known pseudoprimes to bases 2 and 3. The two modular exponentiations are reasonably quick, the lookup by binary search is even quicker, and the primality test is fully deterministic.
Your task is to calculate the pseudoprimes to bases 2 and 3 less than a million. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
http://www.cecm.sfu.ca/Pseudoprimes/ has a list of base-2 pseudoprimes less than 2^64.
That list has 118968378 entries. Filtering for base-3 as well leaves 18994162.
Here’s a solution in Python that uses trial division for checking primality.
[sourcode lang=”python”]
def is_prime(n):
x = 2
while x * x <= n:
if n % x == 0:
return False
x += 1 + (x & 1)
return n > 1
def expm(base, exponent, modulus):
out = 1
base %= modulus
while exponent > 0:
if exponent & 1:
out = (out * base) % modulus
exponent >>= 1
base = (base * base) % modulus
return out
for x in range(3, 1_000_000):
if 1 == expm(2, x – 1, x) == expm(3, x – 1, x) and not is_prime(x):
print(x)
[/sourcecode]
Output:
Here it is again, hopefully properly formatted this time.
def is_prime(n): x = 2 while x * x <= n: if n % x == 0: return False x += 1 + (x & 1) return n > 1 def expm(base, exponent, modulus): out = 1 base %= modulus while exponent > 0: if exponent & 1: out = (out * base) % modulus exponent >>= 1 base = (base * base) % modulus return out for x in range(3, 1_000_000): if 1 == expm(2, x - 1, x) == expm(3, x - 1, x) and not is_prime(x): print(x)Output: