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
    

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: