1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
|
#!/usr/bin/env perl
# Verify prime gaps, version 1.0
# Dana Jacobsen, 2014
#
# This is an alternative to T.R. Nicely's cglp4 program from:
# http://www.trnicely.net/#Downloads
# This runs 2-4x faster on my machines. If cglp4 can use PFGW, then it will
# cross over speed around 3000 digits, and PFGW is much faster at 10k+.
#
# It will use the extra-strong BPSW test plus a Frobenius-Underwood test
# for the endpoints so is more stringent about endpoint testing (cglp4 uses
# the strong BPSW test).
#
# The gaps are in one of the formats:
# <gapsize> <P1-expression>
# <gapsize> <merit> <P1-expression>
# <merit> <gapsize> PRP#### = <P1-expression>
#
# This program will DIE if an invalid gap is found. I believe this is
# preferable to printing a 0 result in a list which may be thousands of
# lines long, and hence missed. If the gaps have been properly supplied,
# this should never come up.
use warnings;
use strict;
use Math::BigInt lib=>"GMP";
use Math::Prime::Util qw/:all/;
use Math::Prime::Util::GMP; # Ensure we're using this
use Time::HiRes qw(gettimeofday tv_interval);
$|=1;
# TODO: Use a command line argument
my $use_pfgw = 0;
#my $pfgw_exec = "/users/jacobsen/src/pfgw-3.7.10/pfgw64";
my $pfgw_exec = "pfgw64";
my $pfgw_thresh = 2400; # PFGW faster only for this many digits
my $fstart = [gettimeofday];
my $procn = 0;
while (<>) {
chomp;
next if /^#/ || /^\s*$/;
my($mer, $gap, $expr);
if (/^\s*(\d+) (\S+) (\S+)$/) {
($mer, $gap, $expr) = ($2, $1, $3);
} elsif (/^\s*(\S+)\s+(\d+)\s+PRP\d+ = (.*)/) {
($mer, $gap, $expr) = ($1, $2, $3);
} elsif (/^(\d+) (\S+)$/) {
($gap, $expr) = ($1, $2);
} else {
warn "skipping $_\n";
next;
}
$procn++;
my $start = [gettimeofday];
$expr =~ s/^1\*//;
my $orig_expr = $expr;
my $n = numerate($expr);
my $end = $n + $gap;
my $dstr = length($n) . "D";
my $dstr2 = length($end) . "D";
my $log2n = int(length($n) * 3.322); # approx
printf "G=%7d %10.2fs Checking P1 ($dstr)...\r", $gap, tv_interval($start);
die "beg of '$expr' is not prime" unless test($n);
printf "G=%7d %10.2fs Checking P2 ($dstr2)... \r", $gap, tv_interval($start);
die "end of '$expr' is not prime" unless test($end);
my $next;
# To avoid all the overhead of timing and printing, for very small
# gaps we can just call next_prime which will check all the interior
# points. The only downside is that we're losing some manual control.
if (0 && $gap < 15000 && $log2n < 800) {
printf "G=%7d %10.2fs Checking P1 ($dstr) interval... \r", $gap, tv_interval($start);
$next = next_prime($n);
} else {
my $depth = int( 1.2 * $log2n * $log2n * log($log2n) );
printf "G=%7d %10.2fs Sieving to $depth ...%s \r", $gap, tv_interval($start), " " x 30;
my @list = sieve_range($n+1, $gap-1, $depth);
my $gapstart = [gettimeofday];
my $ntests = scalar(@list);
my $i = 0;
my $nexti = 1;
printf "G=%7d %10.2fs Checking P1 ($dstr) + %d... \r", $gap, tv_interval($start), $list[0]-$n;
foreach my $rgap (@list) {
my $pgap = $rgap + 1; # We sieved from $n+1
die "Interior point $expr + $pgap is prime\n" if testint($n+$pgap);
$i++;
if ($i >= $nexti) {
my $startint = tv_interval($start);
my $gaptime = tv_interval($gapstart);
my $est = $startint + ($ntests-$i) * $gaptime/$i;
printf "G=%7d %10.2fs (est %.2fs) Checking P1 ($dstr) + $pgap... \r", $gap, $startint, $est;
my $display_intervals = int(0.4 / ($gaptime/$i));
#$display_intervals = 256 if $display_intervals > 256;
$nexti = $i + $display_intervals;
}
}
$next = $end;
}
if ($next == $end) {
printf "G=%7d P1=%-40sOK BPSW+FU=1 (%.3fs)\n",
$gap, $expr, tv_interval($start);
} else {
die "gap $gap for $expr should be ", $next-$n, "\n";
}
}
printf "\n Errors=0. OK=%d. T=%.3f.\n", $procn, tv_interval($fstart);
sub numerate {
my $expr = shift;
$expr =~ s/\b(\d+)#/primorial($1)/g;
$expr =~ s/\^/**/;
$expr =~ s/(\d+)/ Math::BigInt->new("$1") /g;
my $n = eval $expr;
die "Cannot eval: $expr\n" if !defined $n;
return $n;
}
sub test {
my $n = shift;
return is_bpsw_prime($n) && is_frobenius_underwood_pseudoprime($n);
}
sub testint {
my $n = shift;
if ($use_pfgw && length($n) >= $pfgw_thresh) {
return 0 if system("$pfgw_exec -k -Cquiet -f0 -u0 -q\"$n\" >/dev/null 2>1");
}
return is_bpsw_prime($n) && is_frobenius_underwood_pseudoprime($n);
}
|
