Prime gap programming challenge

[robert44444uk]

Originally Posted by danaj

Which is really just optimizing the next_prime solution. This may be a range where NewPGen or similar could excel.

I don't think NewpGen is a very fast sieve, unless someone has seriously speeded it up.

[robert44444uk]

Originally Posted by danaj

Your original Perl program took 4m30s on my computer. The following took 2m41s, but also changed the parameters (which seems to have made the most difference). ]

Running my revised code has slowed it down. Using 800 and 200 as the jumper and the next prime hurdle respectively now takes 211 seconds on my computer. I am still trying to find the optimal jumper, but 850 and 150 takes only 190 seconds.

Using 1346 and 134 on my code reduces the time on my computer to 133 seconds.

Clearly 1346 was chosen as it is the lowest prime gap which is not classified as first occurrence and 134 is 10% of that. I would suggest the optimal percentage might be a bit higher based on my 850 and 150 experience.

[henryzz]

Originally Posted by danaj

The main issue is that when we're looking for "large" gaps, we really don't care about finding all the primes. Finding the next two primes is one too many. We want to skip as much work as possible. Unlike the large record prime gaps, these really aren't all that big so it's an interesting tradeoff. Re using next_prime+800 that was part of the idea.

My guess is that we should do some mild sieving over the full range first and use that knowledge to selectively apply primality tests (that skip the pretests we're using) in the best way so we do the fewest tests. Which is really just optimizing the next_prime solution. This may be a range where NewPGen or similar could excel.


This is somewhat reminiscent of the large prime gap discussion, which also had the common knowledge that one should do fast sieves, and turned out to be soundly beaten by prev_prime / next_prime. It's worse at larger sizes -- Gapcoin using sieves seems to be about 20x less efficient than prev/next prime (all things are not equal in that comparison).

I agree sieving as a primality test is overkill, however, light sieving plus selective primality tests makes more sense than doing primality tests based upon a mod 30 wheel even if some of the sieving is wasted. Sieving is cheap compared to primality tests.
Sieving large ranges gives the opportunity to cherry pick and attempt to create high merit gaps.

[danaj]

Originally Posted by henryzz

... light sieving plus selective primality tests makes more sense than doing primality tests based upon a mod 30 wheel even if some of the sieving is wasted. Sieving is cheap compared to primality tests.

It's all a tradeoff -- the mod-30 wheel is simple, but the "primality test" does some native mods as well before diving into BPSW. I don't think it's as obviously bad as one might think. The speed of native ops with some asm is quite different than working with large numbers in PFGW gets you used to.

Regardless I think we agree that light sieving is the better answer in the long run for finding maximal gaps (I still think it's the wrong thing for finding record gaps, which is a different problem). Since we know we're doing a big range so we can pre sieve a fairly large swath. Run in segments of size 323323 bytes for instance, and a memcpy at the start has already sieved the whole range up to 19.

For the C code I posted, you could do the same in Perl for a little more overhead. I mostly wanted to get rid of the obvious "OMG, you're dong this in a horrifically slow language so it can never work" responses. The overhead does matter at some point, but for now we have more important (IMO) things like (1) is this correct, (2) is it substantially worse than a demonstrated other solution, (3) is it possible to make it faster and how, etc.

[robert44444uk]

I am useless at perl threads, so I took one of danaj's working thread programs and stripped it back to essentials so that I could at least get a working code, see below.

I have set the bar rather low in terms of reporting results, as I am not confident that I am picking up all gaps over 700. I note that some results will appear more than once because of the fact that the $jumper is less than the target gap so it will jump into the same gap.

I am also not confident at all that the code is efficient. With 7 threads it is not working anywhere near 7x as fast as the stand alone code I posted earlier. Using 800 and 700 as $mingap and $jumper It is working at about 2.12e12 an hour, so for a 50bn range that is maybe 84 seconds, compared to 190 seconds for the non-thread code (albeit that was only looking for gaps of 1000, but the threaded code is certainly nowhere near 7x). Maybe someone could have a look at it and see what I am doing wrong.

#!/usr/bin/env perl
use warnings;
use strict;
use threads;
use threads::shared;

use Math::GMPz;
use Math::BigFloat lib=>"GMP";
use Math::Prime::Util qw/:all/;
$|=1;

my $nthreads = 7;
my $beg = shift || 4_000_000_000_000_000_000;
my $end = shift || 4_001_000_000_000_000_000;
my $jumper = 600;
my $mingap = 700;  
my @results :shared;
  
  my @threads;
  push @threads, threads->create('searchi', $_,$jumper,$mingap) for 0 .. $nthreads-1;
  $_->join() for (@threads);

  while (@results) {
    print shift(@results);
  }

sub searchi {

  my($tnum,$sn,$nth) = @_;
  my $index = $tnum;
  while (1) {
    if ($tnum == 0 && @results) {
      lock(@results);
      while (@results) {
        print shift(@results);
      }
    }
    my $i = $beg+($index*$jumper);
    $index += $nthreads;
    last if $i > $end;
    my $n = $i;
    my $forwardlook = $mingap-$jumper;  
    my $dprev = prev_prime($n);
    my $dnext = next_prime($n);
  next unless $dnext>$n+$forwardlook;
    my $gap = $dnext-$dprev;
  next unless $gap >= $mingap;             
    
   
    {
      lock(@results);
      push @results, sprintf("%d %d\n", $gap, $dprev);
    }
  }
  1;
}

[Antonio]

Originally Posted by robert44444uk

I am useless at perl threads, so I took one of danaj's working thread programs and stripped it back to essentials so that I could at least get a working code, see below.

I have set the bar rather low in terms of reporting results, as I am not confident that I am picking up all gaps over 700. I note that some results will appear more than once because of the fact that the $jumper is less than the target gap so it will jump into the same gap.

I am also not confident at all that the code is efficient. With 7 threads it is not working anywhere near 7x as fast as the stand alone code I posted earlier. Using 800 and 700 as $mingap and $jumper It is working at about 2.12e12 an hour, so for a 50bn range that is maybe 84 seconds, compared to 190 seconds for the non-thread code (albeit that was only looking for gaps of 1000, but the threaded code is certainly nowhere near 7x). Maybe someone could have a look at it and see what I am doing wrong.

Try this which scales correctly to 4 threads on my machine. I added the number of threads as the first item on the command line to speed up testing.

I've used an attachment, the [CODE] wrapper doesn't seem to be working for me at the moment!

[robert44444uk]

Thank you for posting Antonio. I'll check it for speed on my machine as soon as it gets free.

[Antonio]

Here is an alternative version which I estimate would take just over 6 days to cover a 1e15 range looking for gaps >= 1346, This is using 4 cores of my i5-3570k at 4.3GHz.

Each thread tests a portion of the complete range, so results may not be output in order of initiating prime.

#!/usr/bin/env perl
#
use warnings;
use strict;
use threads;
use threads::shared;
use Math::Prime::Util qw/:all/;
$|=1;

my $instances = shift || 4;
$instances = $instances % 100;
my $ps = shift || 4_000_023_100_000_000_000;
my $pe = shift || 4_000_023_150_000_000_000;
$ps = sprintf("%d", $ps);
$pe = sprintf("%d", $pe);
my $target = shift || 1000;

my %pstart;
my %pend;
my @results :shared;

print "\n";
print "Requested Threads.................. = $instances\n";
print "Requested search start............. = $ps\n";
print "Requested search end............... = $pe\n";
print "Requested target gap............... = $target\n\n";

$target = $target - $target % 2; # ensure the gap is even
my $head = int(0.1 * $target);
$head = $head + 1 unless $head % 2; # ensure adder is odd
my $tail = $target - $head;

$ps = prev_prime($ps + 1);
my $step = int(($pe - $ps) / $instances);
$pe = next_prime($pe - 1);
$step = $step + 1 if ($ps + $instances * $step) < $pe;
my $p2e = $ps;

my $instance = 0;
while ($instance < $instances) {
	my $p2 = $ps + $step * ($instance + 1);
	my $p1s = $p2e;
	$p2e = next_prime($p2 - 1);
	$pstart{$instance} = $p1s;
	$pend{$instance} = $p2e;
	$instance++;
}
print "Starting prime to ensure full cover = $pstart{0}\n";
print "Ending prime to ensure full cover.. = $pend{$instances - 1}\n";
print "Searching for gaps................ >= $target\n";
print "Forward skip....................... = $tail\n";
print "Looking for forward gap........... >= $head\n";
print "\n";


my $start = time;
print "started search\n\n";

my @threads;
push @threads, threads->create('searchi', $_) for 0 .. $instances - 1;
$_ -> join() for (@threads);

while (@results) {
	print shift(@results);
}
   
my $duration = time - $start;
print "Execution time: $duration s\n";
my $rate = int(($pe - $ps) / $duration + 0.5);
printf("Rate = %d/sec.\n", $rate);

exit;

sub searchi {
	my $tnum = $_;
	my $n = $pstart{$tnum};
	my $pprime = $n;
	my $nprime = $n;
	my $pe = $pend{$tnum};
	
	if ($tnum){
		while ($n < $pe) {
			$n = $nprime + $tail;
			$nprime = Math::Prime::Util::next_prime( $n );
			next unless ($nprime - $n) >= $head;

			$pprime = Math::Prime::Util::prev_prime( $n );
			next unless ($nprime - $pprime) >= $target;
			
			my $gap = $nprime - $pprime;
			{
				lock(@results);
				push @results, sprintf ("%8d %s\n", $gap, $pprime);
			}
		}
	}else {
		while ($n < $pe) {
			if (@results) {
				lock(@results);
				while (@results) {
					print shift(@results);
				}
			}
			
			$n = $nprime + $tail;
			$nprime = Math::Prime::Util::next_prime( $n );
			next unless ($nprime - $n) >= $head;

			$pprime = Math::Prime::Util::prev_prime( $n );
			next unless ($nprime - $pprime) >= $target;
			
			my $gap = $nprime - $pprime;
			{
				lock(@results);
				push @results, sprintf ("%8d %s\n", $gap, $pprime);
			}
		}
	}
	1;
}

[robert44444uk]

Originally Posted by Antonio

Here is an alternative version which I estimate would take just over 6 days to cover a 1e15 range looking for gaps >= 1346, This is using 4 cores of my i5-3570k at 4.3GHz.

Wow, this took just 37s on my machine using 7 threads on the 4 cores. A massive improvement. I don't understand why the forward skip should be forced to be an odd number. I would have thought it should be an even number. My thinking is that your forward skip could land on a prime and calculate the wrong gap size. Maybe I am missing a trick here. Also the $head percentage of $target should be a variable, as 10% is probably not optimal (I'm thinking it could be nearer 14%)

It would be useful if you could a programme in a soft exit controlled by a key which records in the output file the lowest integer not tested. Alternatively one which outputs the lowest integer not tested every hour say as the gap results are quite rare.

But, thanks Antonio! You should post this code to http://www.mersenneforum.org/showthread.php?t=22187

Why don't you join in the hunt by choosing a range?

[Antonio]

Originally Posted by robert44444uk

My only comment is that the forward skip should be forced to be an even number.

Since the jump forward is being added to the result of next_prime($n) it MUST be odd to ensure that the next value to be tested is a composite, as discussed in http://www.mersenneforum.org/showpos...&postcount=171


As to posting the code, it is still experimental and I am still trying to fine tune and tidy it up, but I will post a version in the next couple of days, and I may join the hunt at the end of the month when I will have some time and resources available.

[robert44444uk]

Originally Posted by Antonio

Since the jump forward is being added to the result of next_prime($n) it MUST be odd to ensure that the next value to be tested is a composite, as discussed in http://www.mersenneforum.org/showpos...&postcount=171


As to posting the code, it is still experimental and I am still trying to fine tune and tidy it up, but I will post a version in the next couple of days, and I may join the hunt at the end of the month when I will have some time and resources available.

Ah so!!

On the code, I think you will need to break up the range that is covered by each thread into smaller chunks given that 1e15 is a large number and there maybe times when program breaks are required. Perhaps each thread should work on say 5 minutes worth of work, given that the call back for new integer ranges does not eat into processing time.