View Single Post
Old 2020-10-30, 15:39   #109
robert44444uk
 
robert44444uk's Avatar
 
Jun 2003
Oxford, UK

1,933 Posts
Default

Quote:
Originally Posted by henryzz View Post
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.
There are two major approaches at the moment.
  • Deficient primorials of the form A*p#/d form the majority of gap records, given that it is at a midpoint of a relatively smooth range 4 merits wide, and a truly efficient perl code is shown below - thanks to danaj.
  • The other approach is primorial offsets, of the form A*p#+/-offset+/-midpoint where the offset is something that provides for a range if integers that are relatively smooth

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 .. $nthreads-1;
  $_->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;
}
robert44444uk is offline   Reply With Quote