## Pseudoprimes To Bases 2 And 3

### December 10, 2019

Since we already have a primality checker from previous exercises, this is easy:

```(define (psp23? n)
(and (= (expm 2 (- n 1) n) 1)
(= (expm 3 (- n 1) n) 1)
(not (prime? n))))```
```> (filter psp23? (range 2 1000000))
(1105 1729 2465 2701 2821 6601 8911 10585 15841 18721 29341
31621 41041 46657 49141 52633 63973 75361 83333 83665 88561
90751 93961 101101 104653 107185 115921 126217 162401 172081
176149 188461 204001 226801 228241 252601 276013 278545
282133 294409 314821 334153 340561 399001 410041 449065
488881 512461 530881 534061 552721 563473 574561 622909
653333 658801 665281 670033 721801 748657 786961 825265
838201 852841 873181 997633)```

You can run the program at https://ideone.com/chL78D.

I looked for but could not find a list of all `psp23?` less than 264. Does anybody know of one? Or just the count of `psp23?` less than 264?

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
```