This might work ...

If n is a 2-psp (pseudoprime to base 2) and prime p divides n then

Code:

n = p, mod (p*ord(p))

where ord(p) is the multiplicative order of 2 mod p

(

https://www.ams.org/journals/mcom/19...-0572872-7.pdf)

This can be used to sieve for numbers that are possible 2-psp. ord(p) is usually about as large as p so the product is ~p

^{2}. This makes the psp-sieve much faster than a regular sieve.

The algorithm:

1. Find all primes up to sqrt(N) = 2

^{32.5} if we want to check up to 2

^{65}.

2. Precompute ord(p) for all primes.

3. Use regular sieve with small primes. I currently use ~10k primes which is probably a good starting point. This greatly reduces the number of potential primes that need to be checked with a fermat test.

4. Use the psp-sieve with all remaining primes. This guarantees that a number is prime if it passes the fermat test and isn't removed by either sieve. This should be fast because p*ord(p) >> p. In many cases p*ord(p) will be so large that it will only occur a few times in the interval 2

^{64} to 2

^{65}. For these we could precompute a list of possible 2-psp instead of sieving.

5. Fermat test.

I think this will be faster than SPRP + Lucas.

Any thoughts?