20201028, 21:41  #100 
"Seth"
Apr 2019
372_{8} Posts 
Do less work, find more records
I've have recently doubled my rate of finding records with two tricks.
I have a new algorithm that speeds up sieving by ~10,000x. This means deeper sieving meaning fewer PRPs. I have two additional tricks that also give a huge speed advantage. For each 'm' I calculate prob(record gap) given the results of deep sieving, this forms a distribution (see attachment bottom row). I only test records from the top ~20% (why waste time on 'm' with low chance of finding a record). Once I have an m with high prob(record gap) I find prev_prime and recalculate prob(record gap  prev prime), 90% of the time prob(record gap  prev prime) < prob(record gap) / 2 so I skip finding next_prime as the chance of record is lower than going to the next m. I get (.1 * prob(record) + .9 * .75 * prob(record)) / (.1 * 2 + .9 * 1) => 70% of the prob(record) in covered in 50% of the time, in practice it's even better because the faster prev_prime are more likely to skip so the denominator is more like (0.1 * (1.5 + 1) + .9 * 0.5). These tricks mean that I'm finding ~54 records / day vs ~25 records / day (testing all m, both sides). 
20201029, 00:00  #101 
"Curtis"
Feb 2005
Riverside, CA
3^{3}·173 Posts 
This search isn't my cup of tea, but those charts and the work they represent are prehhhhhhtty! Nice work, sir.

20201029, 00:15  #102 
"Seth"
Apr 2019
2·5^{3} Posts 

20201029, 08:05  #103  
Jun 2003
Oxford, UK
2^{2}·13·37 Posts 
Quote:
It would be useful, in this and your other post, #98 if you could explain your variables  so what is "m", "d", "p" "X"  generally, for deficient primorial multiples we have tended to use the convention k*p#/d , n*p#/d or A*p#/d but actually it does not matter what you use as long as it is explained. It is also helpful to provide a small worked example where this helps understanding. Which language are you developing the sieve in? On this group we have tended to use perl, because of dana's work providing Perl solutions for sieves, surround_primes and prev_prime and next_prime. or c++ where c++ coded programs have been made executables. But we are all interested in something that speeds up a sieve 10000 times :) Your results pickup is impressive. 

20201029, 09:46  #104  
"Seth"
Apr 2019
FA_{16} Posts 
Quote:
I wrote the sieve and related tools in c++, my first plan was to try and include it in dana's library but it's a little more special purpose (hopefully I can share later this week) and requires many nongeneric functions (e.g. file writing stuff, database management) that doesn't belong in dana's library IMO. If you're interested in testing it out (or anyone) send me a PM. I use the form N = m * p# / d  a., with m instead of your k.  X refers to the offset from N. similar to a above. sieve length = the range on each side of N that I sieve (for an sieve interval of 2*X+1). Looking at a smaller example Code:
d optimizer for P = 101#  large prime=79  sl=505 (5.0 merit) Optimizing  d = 2#  40 remaining, 322 avg gap  SL insufficient 1.871% of time Optimizing  d = 3#  39 remaining, 192 avg gap  SL insufficient 1.964% of time Optimizing  d = 5#  47 remaining, 146 avg gap  SL insufficient 0.793% of time = 1 / log(101#) / ((11/2) * (11/3) * (11/5) ...) ~= 0.094680 for d = 2# = 2, only coprime modulo is [1] so we test N = 1*101#/2 If we look just at the interval 1*101#/2 + [1, 505] we find 40 coprime X (this includes removing the multiples of 2), the lower side ( [1, 505]) has a similar number. call these X_i. Code:
1 * 101#/2  40  2 4 8 16 32 64 128 206 214 218 226 254 256 262 274 278 298 302 314 326 334 346 358 362 382 386 394 398 412 422 428 436 446 452 454 458 466 478 482 502 Code:
avg gap = sum( (1  prob_prime)^(i1) * prob_prime * X_i SL insufficient = chance that next_prime > SL = (1  prob_prime)^(40) = (1  0.094680)^40 = 0.01871 for d = 3# = 6, there are two coprime residuals [1, 5] These have fewer X (so a larger chance of SL being insufficient) but a smaller average value (so a small expected next_prime). Code:
1 * 101#/6  39  2 6 8 12 18 24 32 36 48 54 72 96 108 128 144 162 192 206 216 218 254 278 288 302 314 324 326 362 384 386 398 422 428 432 446 452 458 482 486 5 * 101#/6  39  4 6 12 16 18 24 36 48 54 64 72 96 108 144 162 192 214 216 226 256 262 274 288 298 324 334 346 358 382 384 394 412 432 436 454 466 478 486 502 They have an average of 47 coprime X_i, so a smaller chance of next_prime exceeding SL. Code:
1 * 101#/30  46  4 10 12 18 24 30 40 48 54 60 64 72 90 100 108 120 144 150 160 162 180 192 214 240 250 262 270 274 288 298 300 324 334 358 360 382 384 394 400 412 432 450 454 478 480 502 7 * 101#/30  48  4 6 10 16 18 24 30 36 40 48 54 60 64 90 96 100 108 120 144 150 160 180 214 216 226 240 250 256 270 274 288 298 300 324 334 346 358 360 384 394 400 436 450 454 466 478 480 486 11 * 101#/30  48  2 8 12 18 20 24 30 32 48 50 54 60 72 80 90 108 120 128 144 150 162 180 192 200 218 240 254 270 278 288 300 302 314 320 324 360 362 384 398 422 428 432 450 452 458 480 482 500 13 * 101#/30  50  4 6 10 12 16 24 30 36 40 54 60 64 72 90 96 100 120 144 150 160 162 180 192 214 216 226 240 250 256 262 270 274 300 324 334 346 360 382 384 394 400 412 432 436 450 454 466 480 486 502 17 * 101#/30  45  6 8 18 20 24 30 36 48 50 54 60 80 90 96 108 120 128 144 150 180 200 206 216 218 240 254 270 278 288 300 314 320 324 326 360 384 386 398 428 446 450 458 480 486 500 19 * 101#/30  45  6 10 12 16 18 30 36 40 48 60 72 90 96 100 108 120 150 160 162 180 192 216 226 240 250 256 262 270 288 298 300 346 358 360 382 400 412 432 436 450 466 478 480 486 502 23 * 101#/30  46  2 6 12 20 24 30 32 36 50 54 60 72 80 90 96 120 144 150 162 180 192 200 206 216 240 254 270 300 302 314 320 324 326 360 362 384 386 422 432 446 450 452 480 482 486 500 29 * 101#/30  50  2 6 8 12 18 20 30 32 36 48 50 60 72 80 90 96 108 120 128 150 162 180 192 200 206 216 218 240 270 278 288 300 302 320 326 360 362 386 398 422 428 432 446 450 452 458 480 482 486 500 The highest chance of having a next_prime gap > 9070 is with d = 5#. Code:
d optimizer for P = 907#  SL=9070 (10.0 merit) Optimizing  d = 1 * 1#  973 remaining, 2809 avg gap  SL insufficient 0.000% of time Optimizing  d = 1 * 2#  688 remaining, 4856 avg gap  SL insufficient 0.006% of time Optimizing  d = 1 * 3#  473 remaining, 4535 avg gap  SL insufficient 0.123% of time Optimizing  d = 1 * 5#  423 remaining, 3541 avg gap  SL insufficient 0.245% of time Optimizing  d = 1 * 7#  432 remaining, 2770 avg gap  SL insufficient 0.214% of time Optimizing  d = 1 * 11#  458 remaining, 2328 avg gap  SL insufficient 0.144% of time Note that this maximization depends on SL too but that relationship is larger to explain. Code:
d optimizer for P = 2503#  SL=25030 (10.0 merit) Optimizing  d = 1 * 1#  2396 remaining, 7864 avg gap  SL insufficient 0.000% of time Optimizing  d = 1 * 2#  1662 remaining, 14483 avg gap  SL insufficient 0.007% of time Optimizing  d = 1 * 3#  1092 remaining, 15502 avg gap  SL insufficient 0.189% of time Optimizing  d = 1 * 5#  918 remaining, 13562 avg gap  SL insufficient 0.512% of time Optimizing  d = 1 * 7#  891 remaining, 11227 avg gap  SL insufficient 0.595% of time Optimizing  d = 1 * 11#  923 remaining, 9492 avg gap  SL insufficient 0.491% of time 

20201029, 10:03  #105  
"Seth"
Apr 2019
2·5^{3} Posts 
Quote:
The 9 graphics are Top Row: 1. probability of gap = X 2. observed gaps from 15,000 tests and average gap (16207) 3. cdf of gap <= X Middle Row: 1. Before the run I calculate the probability of gap > 12 merit. orange line is the sum of those probabilities if I perform the tests in order, green if I perform the tests best to work, blue is the observed count of gaps with merit > 12 (notice how nicely it matches!) 2. Histogram of count of m with prob(gap > 12 merit) 3. CDF of over those probabilities (e.g. 100%95% = 5% of m have more than a 0.0501 chance of having gap with merit greater than 12, 10% have greater than a 0.0474 chance. Bottom Row (same as middle but with Prob(gap is record) instead of Prob(gap > 12 merit). 1. Note that 5 records were found (you can see them here: https://github.com/primegaplistpro...ommit/88f092c1), this exceeds the expectation of finding ~2.2 records. I assume that's just some extra luck (my current run has found 34 and expects 40 from 65000 intervals). 2. Ditto middle row 3. My point in #100 was why test the 50% of records with probability of record < 0.0000951 when we could be testing the 5% of numbers with probability > 0.000148. 

20201029, 20:58  #106  
Dec 2008
you know...around...
2×313 Posts 
Quote:
Wow! I feel obsolete  last time I searched for some gap records my laptop ran for two days without finding any... well, the laptop is already quite obsolete itself, but still... 

20201029, 21:52  #107  
"Seth"
Apr 2019
2·5^{3} Posts 
Quote:
Recent status report Code:
13922707 447 < unknowns > 498 660 < gap > 0 tests 77454 (1.97/sec, 7.1 secs/test) 39334, 550284 secs unknowns 73644513 (avg: 950.82), 99.05% composite 49.99% < % > 50.01% prp tests 7510657 (avg: 86.04/side, 96.97/m) (190.946 tests/sec) side skips 67620 (87.3%) (avg prob_record 0.00054/side, 0.000337/m/2, +59.9%) sum(prob_minmerit): 2.5e+03, 5.55e+03/day found: 143 sum(prob_record): 47, 103/day found: 53 fallback prev_gap 310 (0.4%), next_gap 54 (0.1%) merit 21.251 (at m=13880539) 

20201030, 10:53  #108 
Just call me Henry
"David"
Sep 2007
Cambridge (GMT/BST)
2·2,909 Posts 
Very interesting work on optimizing the choice of d#. Can't wait to try this new code.
This has made me realize that I don't know the recommended scripts for gap searching currently. https://mersenneforum.org/showthread.php?t=20872 could do with updating for the post Thomas Nicely era and adding some links to recommended software/scripts/guides would be useful for any new participants. Is https://mersenneforum.org/showpost.p...&postcount=103 still the latest script? Is ATH's program from that thread still recommended? I have improvements to both that I can share if they aren't obsolete. I realize that some of this will probably be replaced by Seth's work. 
20201030, 15:39  #109  
Jun 2003
Oxford, UK
2^{2}·13·37 Posts 
Quote:
Code:
#!/usr/bin/env perl use warnings; use strict; use threads; use threads::shared; #use Timer::Runtime; use Math::GMPz; use Math::BigFloat lib=>"GMP"; use Math::Prime::Util qw/:all/; $=1; my $nthreads = enter the number of threads you machine can cancel here; my $startprim = shift  enter xth prime here; my $endprime = shift  enter yth prime here; my $beg = shift  enter A here; my $end = shift  enter upper end A of search here; my $mingaplen = enter the smallest gap you want to see here as a merit; my $divisor = enter your divisor here; my $delta = we usually put 4 here, but it will not search this further if a prime comes up in a range of merits + or  delta/2 ; my %merit; my %mindig; open(my $merits, '<', 'merits.txt') or die "Cannot open merits file\n"; while (<$merits>) { next unless /^\s*(\d+)\s+(\S+)/; $merit{$1} = $2; $mindig{$1} = int( $1 / ($2 * 2.3025) + 1.01 ); # Slightly large } close($merits); my @acc = (1,5,7,11,13,17,19,23,25,29); #my @acc = (0..29); my @results :shared; foreach my $pp ($startprim .. $endprime) { my $prim = pn_primorial($pp); $prim = Math::BigInt>new("$prim") unless ref($prim) eq 'Math::BigInt'; my ($sn, $rem) = $prim>copy>bdiv($divisor); next unless $rem == 0; next if ($sn * $end) < 4e18; my $snlog = length($sn) * 2.3026; # Approximately one merit my $mingap = int($mingaplen * $snlog); $mingap = 1432 if $mingap < 1432; my $prev_thresh = int( $delta * $snlog ); my $nth = nth_prime($pp); print "$nth $pp $divisor ($beg..$end) \n"; my @threads; push @threads, threads>create('searchi', $_, $sn,$nth,$prev_thresh,$mingap) for 0 .. $nthreads1; $_>join() for (@threads); while (@results) { print shift(@results); } } sub merit { my($n, $gap) = @_; my $fgap = Math::BigFloat>new("$gap"); my $fn = Math::BigFloat>new("$n"); return sprintf("%.6lf", $fgap / log($fn)); } sub searchi { my($tnum,$sn,$nth,$prev_thresh,$mingap) = @_; $sn = Math::GMPz>new("$sn"); my $beg30 = int($beg/30); my $index = $tnum; while (1) { if ($tnum == 0 && @results) { lock(@results); while (@results) { print shift(@results); } } my $i = 30*($beg30 + int($index/10)) + $acc[$index % 10]; # my $i = 30*($beg30 + int($index/30)) + $acc[$index % 30]; $index += $nthreads; last if $i > $end; my $n = $sn * $i; my($dprev,$dnext)=Math::Prime::Util::GMP::surround_primes($n,$prev_thresh); next if $dprev == 0  $dnext == 0; # prime found inside thresh my $gap = $dprev + $dnext; next unless $gap > $mingap; # gap is obviously too small my $gbeg = $n  $dprev; next if defined $mindig{$gap} && length($gbeg) > $mindig{$gap}; my $merit = merit($gbeg, $gap); next if defined $merit{$gap} && $merit{$gap} >= $merit; if (!is_extra_strong_lucas_pseudoprime($gbeg)  !is_extra_strong_lucas_pseudoprime($gbeg+$gap)) { print "# SPSP2 found for $i*$nth#/$divisor\n"; ($dprev,$dnext) = Math::Prime::Util::GMP::surround_primes($n); $gap = $dprev + $dnext; $gbeg = $n  $dprev; $merit = merit($gbeg, $gap); next if defined $merit{$gap} && $merit{$gap} >= $merit; } { lock(@results); push @results, sprintf("%9.6f %8d PRP%d = %d*%d#/%d%d\n", $merit, $gap, length($n), $i, $nth, $divisor, $dprev); } } 1; } 

20201031, 20:24  #110 
Jan 2018
41 Posts 
Impressive results Seth!
I have been looking at different approaches myself and have been testing those at small integer sizes (<10^20) to see how often they deliver maximum prime gaps (a test of the m efficiency if you will). I have to say that my findings are inconclusive. Skipping m values will result in missing max prime gaps. So I am sticking to much slower search routines that include all m (m = all odd multipliers in center = m * p# / d). If you do want to limit the numbers of m to test, my observation is that m = not squarefree gives more maximum prime gaps than searching for m = prime or m = (squarefree and not prime). If you do want to increase merit, I would suggest the following: I do find that adding an extra q to the center value (center = q * m * p# / d) delivers larger merits. Note: m and q are both prime. Just my 2 cents. Kind regards Michiel 
Thread Tools  
Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Basic Number Theory 4: a first look at prime numbers  Nick  Number Theory Discussion Group  6  20161014 19:38 
Before you post your new theory about prime, remember  firejuggler  Math  0  20160711 23:09 
Mersene Prime and Number Theory  Ricie  Miscellaneous Math  24  20090814 15:31 
online tutoring in prime number theory  jasong  Math  3  20050515 04:01 
Prime Theory  clowns789  Miscellaneous Math  5  20040108 17:09 