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 6c850f8a310be40a48bf437d821b9904fc1c41d1
Author: Dana Jacobsen <d...@acm.org>
Date:   Thu Oct 18 22:56:42 2012 -0600

    Next round of primes.pl mods
---
 bin/primes.pl                     |  50 +++++++++++--
 examples/make-script-test-data.pl | 153 ++++++++++++++++++++++++++++++++++++++
 examples/test-primes-script.pl    | 133 ++++++++++++++++++++++++---------
 lib/Math/Prime/Util.pm            |   2 +-
 4 files changed, 293 insertions(+), 45 deletions(-)

diff --git a/bin/primes.pl b/bin/primes.pl
index f44420d..80bfb7d 100755
--- a/bin/primes.pl
+++ b/bin/primes.pl
@@ -11,6 +11,33 @@ $| = 1;
 #   http://en.wikipedia.org/wiki/List_of_prime_numbers
 #   http://mathworld.wolfram.com/IntegerSequencePrimes.html
 
+# This program shouldn't contain any special knowledge about the series
+# members other than perhaps the start.  It can know patterns, but don't
+# include a static list of the members, for instance.  It should actually
+# compute the entries in a range (though go ahead and be clever about it).
+# Example:
+#   DO     use knowledge that F_k is prime only if k <= 4 or k is prime.
+#   DO     use knowledge that safe primes are <= 7 or congruent to 11 mod 12.
+#   DO NOT use knowledge that fibprime(14) = 19134702400093278081449423917
+
+# The various primorial primes are confusing.  Some things to consider:
+#   1) there are two definitions of primorial: p# and p_n#
+#   2) three sequences:
+#         p where 1+p# is prime
+#         n where 1+p_n# is prime
+#         p_n#+1 where 1+p_n# is prime
+#   3) intersections of sequences (e.g. p_n#+1 and p_n#-1)
+#   4) other sequences like A057705: p where p+1 is an A002110 primorial
+#      plus all the crazy primorial sequences (unlikely to be confused)
+#
+# A005234  p      where p#+1   prime
+# A136351  p#     where p#+1   prime     2,6,30,210,2310,200560490130
+# A014545  n      where p_n#+1 prime     1,2,3,4,5,11,75,171,172
+# A018239  p_n#+1 where p_n#+1 prime
+#
+# A006794  p      where p#-1   prime     3,5,11,13,41,89,317,337
+# A057704  n      where p_n#-1 prime     2,3,5,6,13,24,66,68,167
+
 my $segment_size = 30 * 128_000;   # 128kB
 my %opts;
 GetOptions(\%opts,
@@ -32,9 +59,11 @@ GetOptions(\%opts,
            'pnp1|A005234',
            'pnm1|A006794',
            'euclid|A018239',
+           'nompugmp',   # turn off MPU::GMP for debugging
            'help',
           ) || die_usage();
 die_usage() if exists $opts{'help'};
+Math::Prime::Util::prime_set_config(gmp=>0) if exists $opts{'nompugmp'};
 
 
 # Get the start and end values.  Verify they're positive integers.
@@ -91,17 +120,17 @@ sub lucas_primes {
 sub fibonacci_primes {
   my ($start, $end) = @_;
   my @fprimes;
-  my $prime = 2;
+  my $Fk = 2;
   my $k = 3;
-  while ($prime < $start) {
+  while ($Fk < $start) {
     $k++;
-    $prime = fib($k);
+    $Fk = fib($k);
   }
-  while ($prime <= $end) {
-    push @fprimes, $prime if is_prime($prime);
+  while ($Fk <= $end) {
+    push @fprimes, $Fk if is_prime($Fk);
     # For all but k=4, F_k is prime only when k is prime.
     $k = ($k <= 4)  ?  $k+1  :  next_prime($k);
-    $prime = fib($k);
+    $Fk = fib($k);
   }
   @fprimes;
 }
@@ -163,6 +192,13 @@ sub is_pillai {
   #   foreach my $n (grep { ($p % $_) != 1 } (2 .. $p-1)) {
   #     return 1 if Math::BigInt->new($n)->bfac()->bmod($p) == ($p-1);
   #   }
+  # About 20% faster but sucks memory:
+  #   foreach my $n (grep { ($p % $_) != 1 } (2 .. $p-1)) {
+  #     $fac_c[$n] = Math::BigInt->new($n)->bfac() if !defined $fac_c[$n];
+  #     return 1 if ($fac_c[$n] % $p) == ($p-1);
+  #   }
+  # Once p gets large (say 20,000+), then calculating and storing n! is
+  # unreasonable, and the code below will be much faster.
   my $n_factorial_mod_p = Math::BigInt->bone();
   for my $n (2 .. $p-1) {
     $n_factorial_mod_p->bmul($n)->bmod($p);
@@ -172,7 +208,7 @@ sub is_pillai {
   0;
 }
 
-# Also pretty slow.
+# Not nearly as slow as Pillai, but not fast.
 sub is_good_prime {
   my $p = shift;
   return 0 if $p <= 2; # 2 isn't a good prime
diff --git a/examples/make-script-test-data.pl 
b/examples/make-script-test-data.pl
new file mode 100755
index 0000000..bfee664
--- /dev/null
+++ b/examples/make-script-test-data.pl
@@ -0,0 +1,153 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+use File::Spec::Functions;
+use FindBin;
+use bigint;
+use Data::BitStream::XS;
+use Math::Prime::Util qw/is_prime/;
+$|++; #flush the output buffer after every write() or print() function
+
+
+my @test_data = (
+  [  668, "Mersenne",    "--mersenne",   10**100],
+  [ 7529, "Triplet",     "--triplet",    0],
+  [ 7530, "Quadruplet",  "--quadruplet", 0],
+  [23200, "Cousin",      "--cousin",     0],
+  [23201, "Sexy",        "--sexy",       0],
+  [ 1359, "Twin",        "--twin",       0],
+  [ 5385, "Safe",        "--safe",       0],
+  [ 5384, "SG",          "--sophie",     0],
+  [ 2385, "Palindromic", "--palin",      32_965_656_923],
+  [ 5479, "Lucas",       "--lucas",      0],
+  [ 5478, "Fibonacci",   "--fibonacci",  0],
+  [63980, "Pillai",      "--pillai",     2000],
+  [28388, "Good",        "--good",       20000],
+  [ 2407, "Cuban y+1",   "--cuban1",     0],
+  [ 2648, "Cuban y+2",   "--cuban2",     100_000_000],
+  [ 5234, "Primorial+1", "--pnp1",       2500],
+  [ 6794, "Primorial-1", "--pnm1",       2500],
+  [18239, "Euclid",      "--euclid",     0],
+);
+
+foreach my $test (@test_data) {
+  my $oeis_no = $test->[0];
+  my $filename = sprintf("b%06d.txt", $oeis_no);
+  my $link     = sprintf("http://oeis.org/A%06d/b%06d.txt";, $oeis_no, 
$oeis_no);
+  if (!-r $filename) {
+    warn "Getting $filename from $link\n";
+    qx/wget $link/;
+    die "Could not retrieve.  Bailing\n" unless -r $filename;
+  }
+  my $ref_data = read_oeis(@$test);
+  push @$test, $ref_data;
+}
+
+my $stream = Data::BitStream::XS->new( file => 'script-test-data.bs', mode => 
'w' );
+foreach my $test (@test_data) {
+  test_oeis(@$test);
+}
+$stream->write_close();
+
+sub read_oeis {
+  my($oeis_no, $name, $script_arg, $restrict) = @_;
+  die "Restrict isn't defined for $oeis_no : $name" unless defined $restrict;
+
+  my $filename = sprintf("b%06d.txt", $oeis_no);
+  my $link     = sprintf("http://oeis.org/A%06d/b%06d.txt";, $oeis_no, 
$oeis_no);
+
+  my @ref;
+  {
+    open my $fh, '<', $filename
+        or die "Can't read $filename.\nYou should run:\n  wget $link\n";
+    printf "%12s primes: reading %12s...", $name, $filename;
+    my $char = " ";
+    while (<$fh>) {
+      next unless /^(\d+)\s+(\d+)/;
+      my $v = (length($2) < 20)  ?  $2  :  Math::BigInt->new("$2");
+      if ($restrict > 0 && $v > $restrict) {
+        $char = '*';
+        last;
+      }
+      push @ref, $v;
+    }
+    close $fh;
+    print "$char";
+  }
+  printf " %7d.", scalar @ref;
+  print "  Testing..";
+  if ($ref[-1] > 18446744073709551615) {
+    print ",";
+    # Check for monotonic and primeness
+    foreach my $i (0 .. $#ref) {
+      die "non-prime in $oeis_no $name\n" unless is_prime($ref[$i]);
+      if ($i > 0) {
+        die "non-monotonic sequence in $oeis_no $name ($i $ref[$i-1] 
$ref[$i])\n" if $ref[$i] <= $ref[$i-1];
+        die "even number in $oeis_no $name\n" if ($ref[$i] % 2) == 0;
+      }
+    }
+  } else {
+    no bigint;
+    print ".";
+    # Check for monotonic and primeness
+    foreach my $i (0 .. $#ref) {
+      die "non-prime in $oeis_no $name\n" unless is_prime($ref[$i]);
+      if ($i > 0) {
+        die "non-monotonic sequence in $oeis_no $name\n" if $ref[$i] <= 
$ref[$i-1];
+        die "even number in $oeis_no $name\n" if ($ref[$i] % 2) == 0;
+      }
+    }
+  }
+  print "done\n";
+  return \@ref;
+}
+
+sub test_oeis {
+  my($oeis_no, $name, $script_arg, $restrict, $ref_data) = @_;
+
+  my @ref = @$ref_data;
+  printf "%12s primes: stream..", $name;
+
+  if ($ref[-1] > 18446744073709551615) {
+    print ",";
+    # Store the first two values, then a list of deltas
+    $stream->put_gamma($oeis_no, 1, scalar @ref, $ref[0], $ref[1]);
+    print ".";
+    my @deltas = map { ($ref[$_] - $ref[$_-1] - 2)/2 } (2..$#ref);
+    print ".";
+    # Ugly...  Check for anything really big;
+    my @giant;
+    foreach my $d (@deltas) {
+      if ($d >= 18446744073709551614) {
+        push @giant, $d;
+        $d = 18446744073709551614;
+      }
+    }
+    print ".";
+    my $k = 2;
+    $stream->put_arice($k, @deltas);
+    print ".";
+    # Store giant deltas raw
+    foreach my $d (@giant) {
+      if (ref($d) ne 'Math::BigInt') {
+        warn "big delta $d isn't a bigint.\n";
+        $d = Math::BigInt->new(0);
+      }
+      my $binstr = substr($d->as_bin, 2);
+      $stream->put_gamma(length($binstr));
+      $stream->put_string($binstr);
+    }
+  } else {
+    no bigint;
+    print ".";
+    # Store the first two values, then a list of deltas
+    $stream->put_gamma($oeis_no, 0, scalar @ref, $ref[0], $ref[1]);
+    print ".";
+    my @deltas = map { ($ref[$_] - $ref[$_-1] - 2)/2 } (2..$#ref);
+    print ".";
+    my $k = 2;
+    $stream->put_arice($k, @deltas);
+  }
+
+  print "done\n";
+}
diff --git a/examples/test-primes-script.pl b/examples/test-primes-script.pl
index 12d2090..e91ccdb 100755
--- a/examples/test-primes-script.pl
+++ b/examples/test-primes-script.pl
@@ -3,51 +3,110 @@ use strict;
 use warnings;
 use File::Spec::Functions;
 use FindBin;
+use Time::HiRes qw(gettimeofday tv_interval);
+use bigint;
+use Data::BitStream::XS;
 $|++; #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(63980, "Pillai", "--pillai", 2000);
-test_oeis(28388, "Good", "--good", 20000);
+# Should use a "put_test_string" sort of call on the test data, so we
+# wouldn't need this.
+my @test_data = (
+  [  668, "Mersenne",    "--mersenne",   10**100],
+  [ 7529, "Triplet",     "--triplet",    0],
+  [ 7530, "Quadruplet",  "--quadruplet", 0],
+  [23200, "Cousin",      "--cousin",     0],
+  [23201, "Sexy",        "--sexy",       0],
+  [ 1359, "Twin",        "--twin",       0],
+  [ 5385, "Safe",        "--safe",       0],
+  [ 5384, "SG",          "--sophie",     0],
+  [ 2385, "Palindromic", "--palin",      32_965_656_923],
+  [ 5479, "Lucas",       "--lucas",      0],
+  [ 5478, "Fibonacci",   "--fibonacci",  0],
+  [63980, "Pillai",      "--pillai",     2000],
+  [28388, "Good",        "--good",       20000],
+  [ 2407, "Cuban y+1",   "--cuban1",     0],
+  [ 2648, "Cuban y+2",   "--cuban2",     100_000_000],
+  [ 5234, "Primorial+1", "--pnp1",       2500],
+  [ 6794, "Primorial-1", "--pnm1",       2500],
+  [18239, "Euclid",      "--euclid",     0],
+);
+my %oeis_name = map { $_->[0] => $_->[1] } @test_data;
 
-sub test_oeis {
-  my($oeis_no, $name, $script_arg, $restrict) = @_;
+my $test_data_hash = read_script_data('script-test-data.bs');
 
-  my $filename = sprintf("b%06d.txt", $oeis_no);
-  my $link     = sprintf("http://oeis.org/A%06d/b%06d.txt";, $oeis_no, 
$oeis_no);
+foreach my $test (@test_data) {
+  my $oeis_no = $test->[0];
+  my $test_data = $test_data_hash->{$oeis_no};
+  if (!defined $test_data) {
+    die "No test data found for OEIS $oeis_no : $test->[1] primes\n";
+  }
+  test_oeis(@$test, $test_data);
+}
 
-  my @ref;
-  my @scr;
-  {
-    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 defined $restrict && $2 > $restrict;
-      push @ref, "$2";
+
+
+
+sub read_script_data {
+  my ($filename) = @_;
+
+  my $stream = Data::BitStream::XS->new( file => $filename, mode => 'ro' );
+  my %hash;
+
+  while (!$stream->exhausted) {
+    my ($oeis_no, $is_bigint, $num_entries, @ref) = $stream->get_gamma(5);
+    printf "%12s primes (OEIS A%06d): reading %7d entries..", 
$oeis_name{$oeis_no}, $oeis_no, $num_entries;
+    if ($is_bigint) {
+      print ",";
+      my $k = 2;
+      my @deltas = $stream->get_arice($k, $num_entries-2);
+      print ".";
+      # Check to see if we have any giant deltas
+      foreach my $d (@deltas) {
+        if ($d >= 18446744073709551614) {
+          my $len = $stream->get_gamma;
+          my $binstr = $stream->read_string($len);
+          $d = Math::BigInt->new('0b' . $binstr);
+        }
+      }
+      print ".";
+      my $prev = $ref[1];
+      push @ref, map { $prev = $_*2+$prev+2; } @deltas;
+      $hash{$oeis_no} = \@ref;
+      print ".\n";
+    } else {
+      no bigint;
+      print ".";
+      my $k = 2;
+      my @deltas = $stream->get_arice($k, $num_entries-2);
+      print ".";
+      my $prev = $ref[1];
+      push @ref, map { $prev = $_*2+$prev+2; } @deltas;
+      $hash{$oeis_no} = \@ref;
+      print ".\n";
     }
-    close $fh;
   }
-  printf " %7d.", scalar @ref;
+  \%hash;
+}
+
+sub test_oeis {
+  my($oeis_no, $name, $script_arg, $restrict, $ref_data) = @_;
 
-  printf "  Generating...";
-  @scr = split /\s+/, qx+$FindBin::Bin/../bin/primes.pl $script_arg 1 
$ref[-1]+;
-  printf " %7d.\n", scalar @scr;
+  my @ref = @$ref_data;
+  printf "%12s primes (OEIS A%06d): generating..", $name, $oeis_no;
+
+  #print "\n";
+  #print "reference data:\n";
+  #print "  $_\n" for @ref;
+  #print "primes.pl $script_arg 1 $ref[-1]\n";
+  my $start = [gettimeofday];
+  my @scr = split /\s+/, qx+$FindBin::Bin/../bin/primes.pl $script_arg 1 
$ref[-1]+;
+  {
+    no bigint;
+    my $num_generated = scalar @scr;
+    my $seconds = tv_interval($start);
+    my $msperprime = ($seconds * 1000.0) / $num_generated;
+    printf " %7d. %7.2f ms/prime\n", $num_generated, $msperprime;
+  }
 
   die "Not equal numbers:  ", scalar @ref, " - ", scalar @scr, "\n"
     unless @ref == @scr;
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 4b5a207..b858b03 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -663,7 +663,7 @@ sub primorial {
   }
 
   my $pn = 1;
-  $pn = Math::BigInt->new->bone if defined $bigint::VERSION &&
+  $pn = Math::BigInt->new->bone if defined $Math::BigInt::VERSION &&
         $n >= (($_Config{'maxbits'} == 32) ? 29 : 53);
 
   foreach my $p ( @{ primes($n) } ) {

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

Reply via email to