View Single Post
Old 2018-10-16, 16:41   #103
robert44444uk
 
robert44444uk's Avatar
 
Jun 2003
Oxford, UK

111100011002 Posts
Default

After a couple of days faffing around with ATH's code, I have built up an impressive list of offsets at 79# covering 2200 and 53# covering 1650. I've just started to run the best 79# offset and already I have one record - a gap of 2486 - after two hours of searching.

I have developed a multithread perl code that takes one crt offset and runs with it. That's the one I am using to run the 79#. So the variable taken by each thread is A in A*79#+offset+1100. In the programme below A is $i

I don't have the ability to develop multithread code that can take, say, 20 to 1000 top offsets and run then all over ever increasing A. Is anyone up for that perl coding challenge? Bascially the thread would need pick up from a queue with $crtoffset and $i, rather than just $i. Either that or allow each thread to take a $crtoffset and run through the $i's

Below is the basic, but working, multithread code for the one CRT offset - that functions currently for 47# and 1450. It is just a hack of danaj's work, so you need blank merits.txt file in the same subdirectory. It also has a lot of red herring code that is not really needed, so it could do with a prune But the code could be used as the base for the threaded multi offset program. It won't deal well at present with $crtoffsets > 2^64, I can't get the print function to print $crtoffset accurately.

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 = 4;
my $startprim = shift || 15;
my $beg = shift || 1;
my $end = shift || 100_000_000_000;

my $divisor   = 1;
my $delta     = 6;
my $crtoffset = 252585620027971755;
my $middle = 725; 
my $mingap = 500;

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 = (0..29);

my @results :shared;

foreach my $pp ($startprim  .. $startprim) {
  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(12 * $snlog);    # Approximately 12 merits

 # $mingap = 2 if $mingap < 3;
  my $prev_thresh = int( $delta * $snlog );

  my $nth = nth_prime($pp);
  print "crtoffset=$crtoffset ${nth}# mingap=$mingap ($beg..$end) \n";
  my @threads;
  push @threads, threads->create('searchi', $_, $sn,$nth,$prev_thresh,$mingap,$crtoffset,$middle) 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,$crtoffset,$middle) = @_;
  $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/30)) + $acc[$index % 30];
    $index += $nthreads;
    last if $i > $end;
    my $n = $sn * $i+$crtoffset+$middle;
    next if is_prime($n);
    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-%d\n", $merit, $gap, length($n), $i, $nth, $crtoffset, $middle, $dprev);
    }
  }
  1;
}

Last fiddled with by robert44444uk on 2018-10-16 at 16:44
robert44444uk is offline   Reply With Quote