Edit to above:
There was no need to do the Lucas test at all. If it is a 2PSP it is known to be composite. You could check if it is a 2PSP before doing the PRP test. I would probably do this by adding the list of 2PSPs 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. 
This might work ...
If n is a 2psp (pseudoprime to base 2) and prime p divides n then [CODE]n = p, mod (p*ord(p))[/CODE] where ord(p) is the multiplicative order of 2 mod p ([url]https://www.ams.org/journals/mcom/198035151/S00255718198005728727/S00255718198005728727.pdf[/url]) This can be used to sieve for numbers that are possible 2psp. ord(p) is usually about as large as p so the product is ~p[SUP]2[/SUP]. This makes the pspsieve much faster than a regular sieve. The algorithm: 1. Find all primes up to sqrt(N) = 2[SUP]32.5[/SUP] if we want to check up to 2[SUP]65[/SUP]. 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 pspsieve 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[SUP]64[/SUP] to 2[SUP]65[/SUP]. For these we could precompute a list of possible 2psp instead of sieving. 5. Fermat test. I think this will be faster than SPRP + Lucas. Any thoughts? 
[QUOTE=CraigLo;586706]
The algorithm: 2. Precompute ord(p) for all primes. Any thoughts?[/QUOTE] 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. [CODE] #!/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); [/CODE] 
[QUOTE=CraigLo;586706]The algorithm:
1. Find all primes up to sqrt(N) = 2[SUP]32.5[/SUP] if we want to check up to 2[SUP]65[/SUP]. 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 pspsieve 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[SUP]64[/SUP] to 2[SUP]65[/SUP]. For these we could precompute a list of possible 2psp instead of sieving. 5. Fermat test. [/QUOTE] If you sieve to sqrt(N) you can use a normal sieve, you do not need pspsieve 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[SUP]6[/SUP])[SUP]2[/SUP]. 
But normal sieve to sqrt(N) is slow. That's why we do the Fermat test in the first place. The pspsieve 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 2PRPs 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 pspsieve would still be faster than the 2 PRP approach. 
I guess I do not understand why pspsieve is faster. It looks like you want to pspsieve to 2[SUP]32.5[/SUP] (except the first ~10K primes).
Even if ord(p) are precomputed, why is sieving mod (p*ord(p)) faster than normal sieving mod p ? 
It's faster mainly because there are far fewer memory accesses. The normal sieve needs to eliminate one out of every p numbers. The pspsieve 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 pspsieve 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[SUP]32[/SUP] the normal sieve will still have about 9 memory writes. p*ord(p) will typically be close to 2[SUP]64[/SUP]. 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. 
[QUOTE=Bobby Jacobs;583335]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[SUP]64[/SUP].[/QUOTE]
I decided to take log[SUB]2[/SUB](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[SUB]2[/SUB] 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[SUB]2[/SUB] within 0.01 of an integer are 4302407359 and 4302407713, with a log[SUB]2[/SUB] of 32.002. That is especially interesting because 32 and 64 are both powers of 2. It is an amazing coincidence! 
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.

[QUOTE=CraigLo;586706]This might work ...
If n is a 2psp (pseudoprime to base 2) and prime p divides n then [CODE]n = p, mod (p*ord(p))[/CODE] where ord(p) is the multiplicative order of 2 mod p ([url]https://www.ams.org/journals/mcom/198035151/S00255718198005728727/S00255718198005728727.pdf[/url]) This can be used to sieve for numbers that are possible 2psp. ord(p) is usually about as large as p so the product is ~p[SUP]2[/SUP]. This makes the pspsieve much faster than a regular sieve. The algorithm: 1. Find all primes up to sqrt(N) = 2[SUP]32.5[/SUP] if we want to check up to 2[SUP]65[/SUP]. 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 pspsieve 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[SUP]64[/SUP] to 2[SUP]65[/SUP]. For these we could precompute a list of possible 2psp instead of sieving. 5. Fermat test. I think this will be faster than SPRP + Lucas. Any thoughts?[/QUOTE] 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 [URL="https://github.com/kimwalisch/primesieve"]primesieve[/URL] 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: for line in tqdm.tqdm(psp_f.readlines()): psp2 = int(line.split()[1]) assert pow(2, psp21, 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 [/CODE] This isn't a formal mathematicians sign off but I'm an undergrad math major and this is a moderate experimental verification. 
I don't know about the approach in its entirety, but I can confirm the congruence.
If 2[sup]n1[/sup] == 1 (mod n), and p is a prime for which pn, then 2[sup]n1[/sup] == 1 (mod p). Consequently, if m is the least positive integer for which 2[sup]m[/sup] == 1 (mod p) [m is the multiplicative order of 2 (mod p)] then m divides n1. We thus have n == 1 (mod m). Since 2[sup]p1[/sup] == 1 (mod p) m also divides p1, 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 p1, so m and p are relatively prime. By CRT we have n == p (mod mp). 
All times are UTC. The time now is 14:54. 
Powered by vBulletin® Version 3.8.11
Copyright ©2000  2021, Jelsoft Enterprises Ltd.