View Single Post
Old 2021-09-22, 19:48   #54
SethTro's Avatar
Apr 2019

1011100012 Posts

Originally Posted by CraigLo View Post
This might work ...
If n is a 2-psp (pseudoprime to base 2) and prime p divides n then
n = p, mod (p*ord(p))
where ord(p) is the multiplicative order of 2 mod p

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 ~p2. This makes the psp-sieve much faster than a regular sieve.

The algorithm:
1. Find all primes up to sqrt(N) = 232.5 if we want to check up to 265.
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 264 to 265. 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?
I read the paper you linked and found the statement that psp2(n=p*r) => n === p mod (p * ord2(p)).

I created a simple python program using the wonderful primesieve and sympy to verify this on all psp2 <= 10^12

import primesieve
import sympy
import tqdm

LIMIT = 10 ** 12

def ord2_fast(p):
    return sympy.ntheory.residue_ntheory.n_order(2, p)

primes = primesieve.primes(3, int(sympy.sqrt(LIMIT) + 1))
prime_index = {p: i for i, p in enumerate(primes)}
s, e = ", ".join(map(str, primes[:3])), ", ".join(map(str, primes[-3:]))
print(f"primes(sqrt({LIMIT})) = {s} ... {e} ({len(primes)})")
order = list(map(ord2_fast, primes))
print("Done calculating ord2(primes)")

with open("b001567.txt") as psp_f:
    for line in tqdm.tqdm(psp_f.readlines()):
        psp2 = int(line.split()[1])
        assert pow(2, psp2-1, psp2) == 1, psp2
        for prime in sympy.factorint(psp2):
            if prime <= primes[-1]:
                pi = prime_index[prime]
                if psp2 % prime == 0:
                    assert psp2 % (prime * order[pi]) == prime
This isn't a formal mathematicians sign off but I'm an undergrad math major and this is a moderate experimental verification.
SethTro is offline   Reply With Quote