Thread: New Maximal Gaps View Single Post 2021-09-22, 19:48   #54
SethTro

"Seth"
Apr 2019

1011100012 Posts Quote:
 Originally Posted by CraigLo 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 ~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

Code:
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:
assert psp2 % (prime * order[pi]) == prime  