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?
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: