[libmath-prime-util-perl] 18/50: Major rewrite of primes.pl, and add more filters

Partha P. Mukherjee ppm-guest at moszumanska.debian.org
Thu May 21 18:45:34 UTC 2015


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 <dana at 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



More information about the Pkg-perl-cvs-commits mailing list