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
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: