New Maximal Gaps

[CraigLo]

Edit to above:
There was no need to do the Lucas test at all. If it is a 2-PSP it is known to be composite. You could check if it is a 2-PSP before doing the PRP test. I would probably do this by adding the list of 2-PSPs into the sieve.

Also, it might be inefficient to do a SPRP test and Lucas test in the GPU. The parallel threads in a block need to execute the same instruction so some threads would probably sit idle while others were doing a Lucas test.

[CraigLo]

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
(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². 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⁶⁵.
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⁶⁴ to 2⁶⁵. 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?

[robert44444uk]

Originally Posted by CraigLo

The algorithm:

2. Precompute ord(p) for all primes.


Any thoughts?

In Perl, the function znorder is quite quick. Up to 2^31 the following Perl program was run on my standard desktop in 304.5 seconds, suggesting about 12 minutes for 2^31.5.

#!/usr/bin/env perl
use warnings;
use strict;
use Math::Prime::Util qw/:all/;
use Time::HiRes qw(gettimeofday tv_interval);

my $start = 2;
my $end = 2**31;
my $starttime = [gettimeofday];
my $x;

forprimes {$x = znorder(2, $_);} $start, $end;

my $duration = tv_interval($starttime);
printf("Execution time: %.2fs\n", $duration);

[ATH]

Originally Posted by CraigLo

The algorithm:
1. Find all primes up to sqrt(N) = 2^32.5 if we want to check up to 2⁶⁵.
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⁶⁴ to 2⁶⁵. For these we could precompute a list of possible 2-psp instead of sieving.

5. Fermat test.

If you sieve to sqrt(N) you can use a normal sieve, you do not need psp-sieve or Fermat test? Everything remaining will be prime.

I was thinking of another idea:
When you jump ahead, for example from p to p+1300, and test backwards until you find a PRP, then what about testing backwards until you have found 2 PRPs? (If first PRP found is still larger than p).

Most of the time those 2 PRPs will be found above p+1000 probably, and in case you get back to p in 2 PRPs and find a false gap, it will be double checked afterwards anyway.
At the price of 2 Fermat tests (+ any tests that are negative), we get a very low risk of missing a prime gap. What are the odds of finding 2 real PRPs so close together blocking a record prime gap? Probably near the square of the probability you calculated earlier? 1 in (8*10⁶)².

[CraigLo]

But normal sieve to sqrt(N) is slow. That's why we do the Fermat test in the first place. The psp-sieve should be very fast. Then we only need to normal sieve with a small number of primes and do a single Fermat test to prove primality.

Is my math and logic correct?

Using 2 different 2-PRPs should be much stronger than MR tests to 2 different bases on the same number. How low do we need probability of missing a gap to be before we can say we have proven a maximal gap?

It's possible the psp-sieve would still be faster than the 2 PRP approach.

[ATH]

I guess I do not understand why psp-sieve is faster. It looks like you want to psp-sieve to 2^32.5 (except the first ~10K primes).
Even if ord(p) are precomputed, why is sieving mod (p*ord(p)) faster than normal sieving mod p ?

[CraigLo]

It's faster mainly because there are far fewer memory accesses. The normal sieve needs to eliminate one out of every p numbers. The psp-sieve only needs to eliminate 1 out of every p*ord(p) numbers.

I currently run in blocks of 36E9. If we did normal sieve up to 100k the first prime in the psp-sieve would be p=100003. This would require 36E9/100003 = 360k memory writes. Ord(p) is 100002. 36E9/(p*ord(p)) = 3.6 memory writes per block.

For large numbers near 2³² the normal sieve will still have about 9 memory writes. p*ord(p) will typically be close to 2⁶⁴. This is big enough that there is no need to even sieve. Just pre compute a list of possible psp for these really large numbers.

So instead of sieving with 200M+ numbers with memory writes in the 10s of billions we would sieve with maybe 50M numbers and probably under a million memory writes.

[Bobby Jacobs]

Originally Posted by Bobby Jacobs

Congratulations on finding another gap! It is amazing that there are so many maximal prime gaps so logarithmically close to the binary round number 2⁶⁴.

I decided to take log₂(p) of all of the primes in the maximal prime gaps.

2, 3, 1, 1.0, 1.5849625007211563
3, 5, 2, 1.5849625007211563, 2.321928094887362
7, 11, 4, 2.807354922057604, 3.4594316186372978
23, 29, 6, 4.523561956057013, 4.857980995127573
89, 97, 8, 6.475733430966398, 6.599912842187128
113, 127, 14, 6.820178962415189, 6.988684686772166
523, 541, 18, 9.030667136246942, 9.079484783826816
887, 907, 20, 9.792790294301064, 9.824958740528523
1129, 1151, 22, 10.140829770773001, 10.16867211813223
1327, 1361, 34, 10.373952655370193, 10.410451351503994
9551, 9587, 36, 13.22143607744528, 13.226863716982413
15683, 15727, 44, 13.93691393843731, 13.940955875611515
19609, 19661, 52, 14.259228343849179, 14.263049081612024
31397, 31469, 72, 14.938339094975541, 14.941643713954356
155921, 156007, 86, 17.250455722905738, 17.2512512383879
360653, 360749, 96, 18.46025189899012, 18.46063586999309
370261, 370373, 112, 18.498183071287098, 18.498619405146496
492113, 492227, 114, 18.908630102645986, 18.90896427018037
1349533, 1349651, 118, 20.364028824642652, 20.364154964998352
1357201, 1357333, 132, 20.372202967479005, 20.37234327572016
2010733, 2010881, 148, 20.939290091967138, 20.939396277626003
4652353, 4652507, 154, 22.149529135616888, 22.149576890238947
17051707, 17051887, 180, 24.023412834966237, 24.02342806415923
20831323, 20831533, 210, 24.312251132247198, 24.31226567594332
47326693, 47326913, 220, 25.496150780051064, 25.49615748646031
122164747, 122164969, 222, 26.864252786761337, 26.86425540845062
189695659, 189695893, 234, 27.499111423558197, 27.49911320320056
191912783, 191913031, 248, 27.515875569415407, 27.51587743374218
387096133, 387096383, 250, 28.528116654614454, 28.528117586356245
436273009, 436273291, 282, 28.700655980036846, 28.700656912571894
1294268491, 1294268779, 288, 30.269489783877926, 30.269490104905696
1453168141, 1453168433, 292, 30.436554495809393, 30.4365547857049
2300942549, 2300942869, 320, 31.099577816119524, 31.099578016760077
3842610773, 3842611109, 336, 31.839439703844857, 31.8394398299949
4302407359, 4302407713, 354, 32.002496981951865, 32.002497100656115
10726904659, 10726905041, 382, 33.32051478290403, 33.3205148342804
20678048297, 20678048681, 384, 34.26738097180326, 34.2673809985947
22367084959, 22367085353, 394, 34.38065819502348, 34.38065822043679
25056082087, 25056082543, 456, 34.54444179308049, 34.54444181933634
42652618343, 42652618807, 464, 35.311915256759846, 35.311915272454314
127976334671, 127976335139, 468, 36.897086096100765, 36.897086101376594
182226896239, 182226896713, 474, 37.406944956835446, 37.40694496058811
241160624143, 241160624629, 486, 37.81120341206439, 37.811203414971786
297501075799, 297501076289, 490, 38.11410392914653, 38.11410393152272
303371455241, 303371455741, 500, 38.14229438999324, 38.142294392371014
304599508537, 304599509051, 514, 38.14812265783897, 38.14812266027346
416608695821, 416608696337, 516, 38.599901996646274, 38.599901998433154
461690510011, 461690510543, 532, 38.748135122042186, 38.748135123704586
614487453523, 614487454057, 534, 39.16059259801241, 39.16059259926613
738832927927, 738832928467, 540, 39.42645720880682, 39.42645720986126
1346294310749, 1346294311331, 582, 40.29213096779726, 40.29213096842094
1408695493609, 1408695494197, 588, 40.35749692819665, 40.357496928798845
1968188556461, 1968188557063, 602, 40.84000557906143, 40.8400055795027
2614941710599, 2614941711251, 652, 41.24991592650985, 41.24991592686957
7177162611713, 7177162612387, 674, 42.70655074661664, 42.70655074675213
13829048559701, 13829048560417, 716, 43.65276713583621, 43.6527671359109
19581334192423, 19581334193189, 766, 44.15454430119517, 44.15454430125161
42842283925351, 42842283926129, 778, 45.284100625854215, 45.28410062588041
90874329411493, 90874329412297, 804, 46.368938046534566, 46.36893804654733
171231342420521, 171231342421327, 806, 47.282940127218566, 47.28294012722535
218209405436543, 218209405437449, 906, 47.63270661562181, 47.6327066156278
1189459969825483, 1189459969826399, 916, 50.07922814332398, 50.07922814332509
1686994940955803, 1686994940956727, 924, 50.583377070516875, 50.58337707051766
1693182318746371, 1693182318747503, 1132, 50.588658751647024, 50.58865875164798
43841547845541059, 43841547845542243, 1184, 55.283148252391904, 55.283148252391946
55350776431903243, 55350776431904441, 1198, 55.61945307272224, 55.61945307272227
80873624627234849, 80873624627236069, 1220, 56.16651879039921, 56.16651879039923
203986478517455989, 203986478517457213, 1224, 57.50125113772148, 57.50125113772149
218034721194214273, 218034721194215521, 1248, 57.59733551004153, 57.59733551004153
305405826521087869, 305405826521089141, 1272, 58.08350519916523, 58.08350519916523
352521223451364323, 352521223451365651, 1328, 58.290487730302445, 58.290487730302445
401429925999153707, 401429925999155063, 1356, 58.477925784547146, 58.47792578454716
418032645936712127, 418032645936713497, 1370, 58.53639322594612, 58.53639322594612
804212830686677669, 804212830686679111, 1442, 59.48035496665738, 59.48035496665738
1425172824437699411, 1425172824437700887, 1476, 60.30584258713822, 60.30584258713822
5733241593241196731, 5733241593241198219, 1488, 62.314056782060014, 62.314056782060014
6787988999657777797, 6787988999657779307, 1510, 62.55768993488138, 62.55768993488138
15570628755536096243, 15570628755536097769, 1526, 63.75546100573405, 63.75546100573405
17678654157568189057, 17678654157568190587, 1530, 63.93864225213069, 63.93864225213069
18361375334787046697, 18361375334787048247, 1550, 63.99330792884274, 63.99330792884274
18470057946260698231, 18470057946260699783, 1552, 64.00182219536605, 64.00182219536605
18571673432051830099, 18571673432051831671, 1572, 64.00973762046758, 64.00973762046758

It turns out that there are 3 maximal prime gaps with a log₂ within 0.01 of 64. Other than the 2 in the gap between 2 and 3, where 2 is exactly a power of 2, the only primes with a log₂ within 0.01 of an integer are 4302407359 and 4302407713, with a log₂ of 32.002. That is especially interesting because 32 and 64 are both powers of 2. It is an amazing coincidence!

[CraigLo]

Can any of the mathematicians here confirm that my approach (post #46) is a valid prime test? If it is I will start to write code to see how fast it is.

[SethTro]

Originally Posted by CraigLo

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
(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². 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⁶⁵.
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⁶⁴ to 2⁶⁵. 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.

[Dr Sardonicus]

I don't know about the approach in its entirety, but I can confirm the congruence.

If 2^n-1 == 1 (mod n), and p is a prime for which p|n, then 2^n-1 == 1 (mod p).

Consequently, if m is the least positive integer for which 2^m == 1 (mod p) [m is the multiplicative order of 2 (mod p)] then m divides n-1.

We thus have n == 1 (mod m). Since 2^p-1 == 1 (mod p) m also divides p-1, so p == 1 (mod m), and n == p (== 1) (mod m).

Since by hypothesis p divides n, we have n == p (== 0) (mod p).

Now m divides p-1, so m and p are relatively prime. By CRT we have n == p (mod mp).