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

Pages: 1 2

### 4 Responses to “Pseudoprimes To Bases 2 And 3”

1. matthew said

http://www.cecm.sfu.ca/Pseudoprimes/ has a list of base-2 pseudoprimes less than 2^64.

2. matthew said

That list has 118968378 entries. Filtering for base-3 as well leaves 18994162.

3. Daniel said

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:

```1105
1729
2465
2701
2821
...
825265
838201
852841
873181
997633
```
4. Daniel said

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:

```1105
1729
2465
2701
2821
...
825265
838201
852841
873181
997633
```