This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.13
in repository libmath-prime-util-perl.
commit 1a3ce7e10877bd70c6477558db7ae7b5b961286c
Author: Dana Jacobsen <d...@acm.org>
Date: Thu Oct 18 01:02:14 2012 -0600
Major rewrite of primes.pl, and add more filters
---
bin/primes.pl | 338 ++++++++++++++++++++++++++++++-----------
examples/test-primes-script.pl | 26 +++-
2 files changed, 268 insertions(+), 96 deletions(-)
diff --git a/bin/primes.pl b/bin/primes.pl
index 3cdb6fe..1e08c65 100755
--- a/bin/primes.pl
+++ b/bin/primes.pl
@@ -1,33 +1,41 @@
#!/usr/bin/env perl
use strict;
use warnings;
-use Math::BigInt;
use Getopt::Long;
-use Math::Prime::Util qw/primes next_prime is_prime/;
+use Math::BigInt try => 'GMP';
+use Math::Prime::Util qw/primes is_prime next_prime prev_prime
+ prime_count primorial pn_primorial/;
$| = 1;
# For many more types, see:
# http://en.wikipedia.org/wiki/List_of_prime_numbers
# http://mathworld.wolfram.com/IntegerSequencePrimes.html
-my $show_safe = 0;
-my $show_sophie = 0;
-my $show_twin = 0;
-my $show_lucas = 0;
-my $show_fibonacci = 0;
-my $show_palindromic = 0;
-my $show_usage = 0;
my $segment_size = 30 * 128_000; # 128kB
-
-GetOptions( "safe" => \$show_safe,
- "sophie|sg" => \$show_sophie,
- "twin" => \$show_twin,
- "lucas" => \$show_lucas,
- "fibonacci" => \$show_fibonacci,
- "palindromic|palindrome|palendrome" => \$show_palindromic,
- "help" => \$show_usage,
+my %opts;
+GetOptions(\%opts,
+ 'safe|A005385',
+ 'sophie|sg|A005384',
+ 'twin|A001359',
+ 'lucas|A005479',
+ 'fibonacci|A005478',
+ 'triplet|A007529',
+ 'quadruplet|A007530',
+ 'cousin|A023200',
+ 'sexy|A023201',
+ 'mersenne|A000668',
+ 'palindromic|palindrome|palendrome|A002385',
+ 'pillai|A063980',
+ 'good|A028388',
+ 'cuban1|A002407',
+ 'cuban2|A002648',
+ 'pnp1|A005234',
+ 'pnm1|A006794',
+ 'euclid|A018239',
+ 'help',
) || die_usage();
-die_usage() if $show_usage;
+die_usage() if exists $opts{'help'};
+
# Get the start and end values. Verify they're positive integers.
die_usage() unless @ARGV == 2;
@@ -98,61 +106,137 @@ sub fibonacci_primes {
@fprimes;
}
+sub mersenne_primes {
+ my ($start, $end) = @_;
+ my @mprimes;
+ my $p = 1;
+ while (1) {
+ $p = next_prime($p); # Mp is not prime if p is not prime
+ next if $p > 3 && ($p % 4) == 3 && is_prime(2*$p+1);
+ my $Mp = Math::BigInt->bone->blsft($p)->bsub(1);
+ last if $Mp > $end;
+ # Lucas-Lehmer test would be faster
+ push @mprimes, $Mp if $Mp >= $start && is_prime($Mp);
+ }
+ @mprimes;
+}
-my @p;
-while ($start <= $end) {
+sub euclid_primes {
+ my ($start, $end, $add) = @_;
+ my @eprimes;
+ my $k = 0;
+ while (1) {
+ my $primorial = pn_primorial(Math::BigInt->new($k)) + $add;
+ last if $primorial > $end;
+ push @eprimes, $primorial if $primorial >= $start && is_prime($primorial);
+ $k++;
+ }
+ @eprimes;
+}
- # This will need to be made more generic if we add more options like this.
- # FIXME: It also doesn't work for our twin prime selection.
- if ($show_lucas && $show_fibonacci) {
- my %l;
- undef @l{ lucas_primes($start, $end) };
- @p = grep { exists $l{$_} } fibonacci_primes($start, $end);
- $start = $end+1;
+# This is not a general palindromic digit function!
+sub ndig_palindromes {
+ my $digits = shift;
+ return (2,3,5,7,9) if $digits == 1;
+ return (11) if $digits == 2;
+ return () if ($digits % 2) == 0;
- } elsif ($show_lucas) {
- @p = lucas_primes($start, $end);
- $start = $end+1;
+ my @prefixes = (1,3,7,9);
+ my $inner_digits = ($digits-1) / 2 - 1;
+ foreach my $d (1 .. $inner_digits) {
+ @prefixes = map { ($_.'0', $_.'1', $_.'2', $_.'3', $_.'4',
+ $_.'5', $_.'6', $_.'7', $_.'8', $_.'9',); } @prefixes;
+ }
+ return map { my $r = reverse($_);
+ ($_.'0'.$r, $_.'1'.$r, $_.'2'.$r, $_.'3'.$r,
+ $_.'4'.$r, $_.'5'.$r, $_.'6'.$r, $_.'7'.$r,
+ $_.'8'.$r, $_.'9'.$r,);
+ } @prefixes;
+}
- } elsif ($show_fibonacci) {
- @p = fibonacci_primes($start, $end);
- $start = $end+1;
+# See: http://en.wikipedia.org/wiki/Pillai_prime
+# This is quite slow.
+sub is_pillai {
+ my $p = shift;
+ return 0 if $p < 2;
+ # return 0 unless is_prime($p);
+ my $n_factorial_mod_p = Math::BigInt->bone();
+ for my $n (2 .. $p-1) {
+ $n_factorial_mod_p->bmul($n)->bmod($p);
+ next if $p % $n == 1;
+ return 1 if $n_factorial_mod_p == ($p-1);
+ }
+ 0;
+}
- } else {
- # small segment size if we're doing bigints
- $segment_size = 10000 if $start > ~0;
+# Also slow.
+sub is_good_prime {
+ my $p = shift;
+ return 0 if $p <= 2; # 2 isn't a good prime
+ my $lower = $p;
+ my $upper = $p;
+ while ($lower > 2) {
+ $lower = prev_prime($lower);
+ $upper = next_prime($upper);
+ return 0 if ($p*$p) <= ($upper * $lower);
+ }
+ 1;
+}
- if ( $show_palindromic && $start >= 100 && (length($start) % 2) == 0 ) {
- # anything with an even number of digits is divisible by 11
- $start = 10 ** (length($start)) + 1;
- }
+sub merge_primes {
+ my ($genref, $pref, $name, @primes) = @_;
+ if (!defined $$genref) {
+ @$pref = @primes;
+ $$genref = $name;
+ } else {
+ my %f;
+ undef @f{ @primes };
+ @$pref = grep { exists $f{$_} } @$pref;
+ }
+}
- my $seg_start = $start;
- my $seg_end = $start + $segment_size;
- $seg_end = $end if $end < $seg_end;
- $start = $seg_end+1;
- # Skip ranges where there are no palindromic primes.
- # - anything starting with 24568 won't be a prime when reversed
- if ( $show_palindromic && $seg_start >= 100 &&
- length($seg_start) == length($seg_end)) {
- next if $seg_start =~ /^2/ && $seg_end =~ /^2/;
- next if $seg_start =~ /^8/ && $seg_end =~ /^8/;
- next if $seg_start =~ /^[456]/ && $seg_end =~ /^[456]/;
+# This is used for things that can generate a filtered list faster than
+# searching through all primes in the range.
+sub gen_and_filter {
+ my ($start, $end) = @_;
+ my $gen;
+ my @p;
+
+ if (exists $opts{'lucas'}) {
+ merge_primes(\$gen, \@p, 'lucas', lucas_primes($start, $end));
+ }
+ if (exists $opts{'fibonacci'}) {
+ merge_primes(\$gen, \@p, 'fibonacci', fibonacci_primes($start, $end));
+ }
+ if (exists $opts{'mersenne'}) {
+ merge_primes(\$gen, \@p, 'mersenne', mersenne_primes($start, $end));
+ }
+ if (exists $opts{'euclid'}) {
+ merge_primes(\$gen, \@p, 'euclid', euclid_primes($start, $end, 1));
+ }
+ if (exists $opts{'palindromic'}) {
+ if (!defined $gen) {
+ foreach my $d (length($start) .. length($end)) {
+ push @p, grep { $_ >= $start && $_ <= $end && is_prime($_) }
+ ndig_palindromes($d);
}
- @p = @{primes($seg_start, $seg_end)};
- #warn "-- ", scalar @p, " primes between $start and $seg_end\n";
+ $gen = 'palindromic';
+ } else {
+ @p = grep { $_ eq reverse $_; } @p;
}
+ }
- next unless scalar @p > 0;
-
+ if (!defined $gen) {
+ @p = @{primes($start, $end)};
+ $gen = 'primes';
+ }
- # Restrict to twin primes if requested.
- if ($show_twin) {
- # Add a last element so we can look at it. Simple: push next_prime.
+ if (exists $opts{'twin'}) {
+ if ($gen ne 'primes') {
+ @p = grep { is_prime( $_+2 ); } @p;
+ } elsif (scalar @p > 0) {
+ # All primes in the range are here, so just look in the array.
push @p, is_prime($p[-1]+2) ? $p[-1]+2 : 0;
- # Trivial but slow:
- #@p = grep { is_prime( $_+2 ); } @p;
- # Loop over @p looking for values with prime == next+2
my @twin;
my $prime = shift @p;
foreach my $next (@p) {
@@ -160,34 +244,88 @@ while ($start <= $end) {
$prime = $next;
}
@p = @twin;
- #warn " reduced to ", scalar @p, " twin primes\n";
}
+ }
- # Restrict to safe primes if requested.
- if ($show_safe) {
- @p = grep { is_prime( ($_-1) >> 1 ); }
- grep { ($_ <= 7) || ($_ % 12) == 11; } # Optimization
- @p;
- #warn " reduced to ", scalar @p, " safe primes\n";
- }
+ if (exists $opts{'triplet'}) { # could be optimized like twin
+ @p = grep { is_prime($_+6) && (is_prime($_+2) || is_prime($_+4)); } @p;
+ }
- # Restrict to Sophie Germain primes if requested.
- if ($show_sophie) {
- @p = grep { is_prime( 2*$_+1 ); }
- grep { my $m30 = $_ % 30;
- $_ <= 5 || $m30 == 11 || $m30 == 23 || $m30 == 29 ; }
- @p;
- #warn " reduced to ", scalar @p, " SG primes\n";
- }
+ if (exists $opts{'quadruplet'}) { # could be optimized like twin
+ @p = grep { is_prime($_+2) && is_prime($_+6) && is_prime($_+8); }
+ grep { $_ <= 5 || ($_ % 30) == 11; }
+ @p;
+ }
- # Restrict to Palendromic primes if requested.
- if ($show_palindromic) {
- @p = grep { $_ eq reverse $_; } @p;
- #warn " reduced to ", scalar @p, " Palindromic primes\n";
+ if (exists $opts{'cousin'}) { # could be optimized like twin
+ @p = grep { is_prime($_+4); }
+ grep { ($_ <= 3) || ($_ % 6) == 1; }
+ @p;
+ }
+
+ if (exists $opts{'sexy'}) { # could be optimized like twin
+ @p = grep { is_prime($_+6); } @p;
+ }
+
+ if (exists $opts{'safe'}) {
+ @p = grep { is_prime( ($_-1) >> 1 ); }
+ grep { ($_ <= 7) || ($_ % 12) == 11; }
+ @p;
+ }
+ if (exists $opts{'sophie'}) {
+ @p = grep { is_prime( 2*$_+1 ); }
+ grep { my $m30 = $_ % 30;
+ $_ <= 5 || $m30 == 11 || $m30 == 23 || $m30 == 29 ; }
+ @p;
+ }
+ if (exists $opts{'cuban1'}) {
+ @p = grep { my $n = sqrt((4*$_-1)/3); $n == int($n); } @p;
+ }
+ if (exists $opts{'cuban2'}) {
+ @p = grep { my $n = sqrt(($_-1)/3); $n == int($n); } @p;
+ }
+ if (exists $opts{'pnm1'}) {
+ @p = grep { is_prime( primorial(Math::BigInt->new($_))-1 ) } @p;
+ }
+ if (exists $opts{'pnp1'}) {
+ @p = grep { is_prime( primorial(Math::BigInt->new($_))+1 ) } @p;
+ }
+ if (exists $opts{'pillai'}) {
+ @p = grep { is_pillai($_); } @p;
+ }
+ if (exists $opts{'good'}) {
+ @p = grep { is_good_prime($_); } @p;
+ }
+ @p;
+}
+
+if (exists $opts{'lucas'} ||
+ exists $opts{'fibonacci'} ||
+ exists $opts{'palindromic'} ||
+ exists $opts{'euclid'} ||
+ exists $opts{'mersenne'}) {
+ my @p = gen_and_filter($start, $end);
+ print join("\n", @p), "\n" if scalar @p > 0;
+} else {
+ my @p;
+ while ($start <= $end) {
+
+ # Adjust segment sizes for some cases
+ $segment_size = 10000 if $start > ~0; # small if doing bigints
+ if (exists $opts{'pillai'}) {
+ $segment_size = ($start < 10000) ? 100 : 1000; # very small for Pillai
}
+ my $seg_start = $start;
+ my $seg_end = $start + $segment_size;
+ $seg_end = $end if $end < $seg_end;
+ $start = $seg_end+1;
+
+ @p = gen_and_filter($seg_start, $seg_end);
+
# print this segment
print join("\n", @p), "\n" if scalar @p > 0;
+ }
}
@@ -201,15 +339,35 @@ Displays all primes between the positive integers START
and END, inclusive.
The START and END values must be integers, however the shortcut "x**y" may
be used, which allows very large values (e.g. '10**500' or '2**64')
-Options:
- --help displays this help message
- --twin displays only twin primes, where p+2 is also prime
- --safe displays only safe primes, where (p-1)/2 is also prime
- --sophie displays only Sophie Germain primes, where 2p+1 is also prime
- --lucas displays only Lucas primes, where L_n is prime
- --fibonacci displays only Fibonacci primes, where F_n is prime
- --palindr displays only Palindromic primes, where p eq reverse p
+General options:
+
+ --help displays this help message
+
+Filter options, which will cause the list of primes to be further filtered
+to only those primes additionally meeting these conditions:
+
+ --twin Twin p+2 is prime
+ --triplet Triplet p+6 and (p+2 or p+4) are prime
+ --quadruplet Quadruplet p+2, p+6, and p+8 are prime
+ --cousin Cousin p+4 is prime
+ --sexy Sexy p+6 is prime
+ --safe Safe (p-1)/2 is also prime
+ --sophie Sophie Germain 2p+1 is also prime
+ --lucas Lucas L_p is prime
+ --fibonacci Fibonacci F_p is prime
+ --mersenne Mersenne M_p = 2^p-1 is prime
+ --palindr Palindromic p is equal to p with its base-10 digits
reversed
+ --pillai Pillai n! % p = p-1 and p % n != 1 for some n
+ --good Good p_n^2 > p_{n-i}*p_{n+i} for all i in (1..n-1)
+ --cuban1 Cuban (y+1) p = (x^3 - y^3)/(x-y), x=y+1
+ --cuban2 Cuban (y+2) p = (x^3 - y^3)/(x-y), x=y+2
+ --pnp1 Primorial+1 p#+1 is prime
+ --pnm1 Primorial-1 p#-1 is prime
+ --euclid Euclid pn#+1 is prime
+
+Note that options can be combined, e.g. display only safe twin primes.
+In all cases involving multiples (twin, triplet, etc.), the value returned
+is p -- the least value of the set.
-Note that options can be combined, e.g. display all safe twin primes.
EOU
}
diff --git a/examples/test-primes-script.pl b/examples/test-primes-script.pl
old mode 100644
new mode 100755
index 1947073..12d2090
--- a/examples/test-primes-script.pl
+++ b/examples/test-primes-script.pl
@@ -5,15 +5,27 @@ use File::Spec::Functions;
use FindBin;
$|++; #flush the output buffer after every write() or print() function
+test_oeis(668, "Mersenne", "--mersenne", '1' . '0' x 100);
+#test_oeis(2407, "Cuban y+1", "--cuban1");
+#test_oeis(2648, "Cuban y+2", "--cuban2", 100_000_000);
+test_oeis(5234, "Primorial+1", "--pnp1", 2500);
+test_oeis(6794, "Primorial-1", "--pnm1", 2500);
+test_oeis(18239, "Euclid", "--euclid");
+test_oeis(7529, "Triplet", "--triplet");
+test_oeis(7530, "Quadruplet", "--quadruplet");
+test_oeis(23200, "Cousin", "--cousin");
+test_oeis(23201, "Sexy", "--sexy");
test_oeis(1359, "Twin", "--twin");
test_oeis(5385, "Safe", "--safe");
test_oeis(5384, "SG", "--sophie");
+test_oeis(2385, "Palindromic", "--palin");
test_oeis(5479, "Lucas", "--lucas");
test_oeis(5478, "Fibonacci", "--fibonacci");
-test_oeis(2385, "Palindromic", "--palin");
+test_oeis(63980, "Pillai", "--pillai", 2000);
+test_oeis(28388, "Good", "--good", 20000);
sub test_oeis {
- my($oeis_no, $name, $script_arg) = @_;
+ my($oeis_no, $name, $script_arg, $restrict) = @_;
my $filename = sprintf("b%06d.txt", $oeis_no);
my $link = sprintf("http://oeis.org/A%06d/b%06d.txt";, $oeis_no,
$oeis_no);
@@ -23,17 +35,19 @@ sub test_oeis {
{
open my $fh, '<', $filename
or die "Can't read $filename.\nYou should run:\n wget $link\n";
+ printf "%12s primes: reading %15s...", $name, $filename;
while (<$fh>) {
next unless /^(\d+)\s+(\d+)/;
- last if $2 >= 10_000_000_000 && $oeis_no == 2385;
- push @ref, $2;
+ last if defined $restrict && $2 > $restrict;
+ push @ref, "$2";
}
close $fh;
}
- warn "Read ", scalar @ref, " $name primes from reference file\n";
+ printf " %7d.", scalar @ref;
+ printf " Generating...";
@scr = split /\s+/, qx+$FindBin::Bin/../bin/primes.pl $script_arg 1
$ref[-1]+;
- warn "Generated ", scalar @scr, " $name primes using primes.pl\n";
+ printf " %7d.\n", scalar @scr;
die "Not equal numbers: ", scalar @ref, " - ", scalar @scr, "\n"
unless @ref == @scr;
--
Alioth's /usr/local/bin/git-commit-notice on
/srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git
_______________________________________________
Pkg-perl-cvs-commits mailing list
Pkg-perl-cvs-commits@lists.alioth.debian.org
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits
