This is an automated email from the git hooks/post-receive script.

ppm-guest pushed a commit to annotated tag v0.29
in repository libmath-prime-util-perl.

commit 83472c086bb7645d0c8de4f9cabce2185a463559
Author: Dana Jacobsen <d...@acm.org>
Date:   Thu May 30 13:00:13 2013 -0700

    Update some benchmarks and examples
---
 MANIFEST                                | 12 ++---
 examples/bench-primearray.pl            |  9 ++--
 examples/bench-primecount.pl            | 40 +++++++++++------
 examples/sophie_germain.pl              | 51 ++++++++++++++++++----
 examples/twin_primes.pl                 | 38 +++++++++++-----
 lib/Math/Prime/Util.pm                  | 15 +++++++
 xt/small-is-next-prev.pl                | 77 +++++++++++++++++++++++++++------
 {examples => xt}/test-bpsw.pl           |  0
 {examples => xt}/test-factor-mpxs.pl    |  8 ++--
 {examples => xt}/test-nthapprox.pl      |  0
 {examples => xt}/test-pcapprox.pl       |  0
 {examples => xt}/test-primes-script.pl  |  0
 {examples => xt}/test-primes-script2.pl |  0
 13 files changed, 187 insertions(+), 63 deletions(-)

diff --git a/MANIFEST b/MANIFEST
index d1c309b..b86ab69 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -47,19 +47,13 @@ examples/bench-mp-nextprime.pl
 examples/bench-mp-psrp.pl
 examples/bench-mp-prime_count.pl
 examples/test-factor-yafu.pl
-examples/test-factor-mpxs.pl
 examples/test-nextprime-yafu.pl
 examples/test-primes-yafu.pl
-examples/test-nthapprox.pl
-examples/test-pcapprox.pl
 examples/sophie_germain.pl
 examples/twin_primes.pl
 examples/find_mr_bases.pl
 examples/parallel_fibprime.pl
-examples/test-bpsw.pl
 examples/test-factor-gnufactor.pl
-examples/test-primes-script.pl
-examples/test-primes-script2.pl
 examples/verify-gmp-ecpp-cert.pl
 examples/verify-sage-ecpp-cert.pl
 bin/primes.pl
@@ -104,4 +98,10 @@ xt/make-script-test-data.pl
 xt/pari-totient-moebius.pl
 xt/rwh_primecount.py
 xt/rwh_primecount_numpy.py
+xt/test-bpsw.pl
+xt/test-factor-mpxs.pl
+xt/test-nthapprox.pl
+xt/test-pcapprox.pl
+xt/test-primes-script.pl
+xt/test-primes-script2.pl
 .travis.yml
diff --git a/examples/bench-primearray.pl b/examples/bench-primearray.pl
index bc2d402..028288e 100755
--- a/examples/bench-primearray.pl
+++ b/examples/bench-primearray.pl
@@ -2,7 +2,7 @@
 use strict;
 use warnings;
 use Math::Prime::Util qw/:all/;
-use Math::Prime::Util::PrimeArray qw/make_prime_iterator/;
+use Math::Prime::Util::PrimeArray;
 use Math::NumSeq::Primes;
 use Math::Prime::TiedArray;
 use Benchmark qw/:all/;
@@ -23,9 +23,6 @@ cmpthese($count,{
   'iterator'  => sub { $s=0; my $it = prime_iterator();
                        $s += $it->() for 0..$ilimit;
                        die unless $s == $expect; },
-  'pa iter'   => sub { $s=0; my $it = make_prime_iterator();
-                       $s += $it->() for 0..$ilimit;
-                       die unless $s == $expect; },
   'pa index'  => sub { $s=0; tie my @primes, "Math::Prime::Util::PrimeArray";
                        $s += $primes[$_] for 0..$ilimit;
                        die unless $s == $expect; },
@@ -82,7 +79,7 @@ cmpthese($count,{
 });
 }
 
-if (0) {
+if (1) {
 print "Walk primes backwards from 1M\n";
 $nlimit = 1_000_000;
 $ilimit = prime_count($nlimit)-1;
@@ -100,7 +97,7 @@ cmpthese($count,{
 });
 }
 
-if (0) {
+if (1) {
 print "Random walk in 1M\n";
 srand(29);
 my @rindex;
diff --git a/examples/bench-primecount.pl b/examples/bench-primecount.pl
index 32b2f0b..3ec4dcd 100755
--- a/examples/bench-primecount.pl
+++ b/examples/bench-primecount.pl
@@ -11,20 +11,34 @@ my $count = shift || -5;
 srand(29);
 my @darray;
 push @darray, [gendigits($_)]  for (2 .. 10);
+my $sum;
 
-  my $sum;
-  cmpthese($count,{
-    ' 2' => sub { $sum += prime_count($_) for @{$darray[2-2]} },
-    ' 3' => sub { $sum += prime_count($_) for @{$darray[3-2]} },
-    ' 4' => sub { $sum += prime_count($_) for @{$darray[4-2]} },
-    ' 5' => sub { $sum += prime_count($_) for @{$darray[5-2]} },
-    ' 6' => sub { $sum += prime_count($_) for @{$darray[6-2]} },
-    ' 7' => sub { $sum += prime_count($_) for @{$darray[7-2]} },
-    ' 8' => sub { $sum += prime_count($_) for @{$darray[8-2]} },
-    #' 9' => sub { $sum += prime_count($_) for @{$darray[9-2]} },
-    #'10' => sub { $sum += prime_count($_) for @{$darray[10-2]} },
-  });
-  print "\n";
+print "Direct sieving:\n";
+cmpthese($count,{
+    ' 2' => sub { $sum += Math::Prime::Util::_XS_prime_count($_) for 
@{$darray[2-2]} },
+    ' 3' => sub { $sum += Math::Prime::Util::_XS_prime_count($_) for 
@{$darray[3-2]} },
+    ' 4' => sub { $sum += Math::Prime::Util::_XS_prime_count($_) for 
@{$darray[4-2]} },
+    ' 5' => sub { $sum += Math::Prime::Util::_XS_prime_count($_) for 
@{$darray[5-2]} },
+    ' 6' => sub { $sum += Math::Prime::Util::_XS_prime_count($_) for 
@{$darray[6-2]} },
+    ' 7' => sub { $sum += Math::Prime::Util::_XS_prime_count($_) for 
@{$darray[7-2]} },
+    ' 8' => sub { $sum += Math::Prime::Util::_XS_prime_count($_) for 
@{$darray[8-2]} },
+    #' 9' => sub { $sum += Math::Prime::Util::_XS_prime_count($_) for 
@{$darray[9-2]} },
+    #'10' => sub { $sum += Math::Prime::Util::_XS_prime_count($_) for 
@{$darray[10-2]} },
+});
+print "\n";
+print "Direct Lehmer:\n";
+cmpthese($count,{
+    ' 2' => sub { $sum += Math::Prime::Util::_XS_lehmer_pi($_) for 
@{$darray[2-2]} },
+    ' 3' => sub { $sum += Math::Prime::Util::_XS_lehmer_pi($_) for 
@{$darray[3-2]} },
+    ' 4' => sub { $sum += Math::Prime::Util::_XS_lehmer_pi($_) for 
@{$darray[4-2]} },
+    ' 5' => sub { $sum += Math::Prime::Util::_XS_lehmer_pi($_) for 
@{$darray[5-2]} },
+    ' 6' => sub { $sum += Math::Prime::Util::_XS_lehmer_pi($_) for 
@{$darray[6-2]} },
+    ' 7' => sub { $sum += Math::Prime::Util::_XS_lehmer_pi($_) for 
@{$darray[7-2]} },
+    ' 8' => sub { $sum += Math::Prime::Util::_XS_lehmer_pi($_) for 
@{$darray[8-2]} },
+    ' 9' => sub { $sum += Math::Prime::Util::_XS_lehmer_pi($_) for 
@{$darray[9-2]} },
+    '10' => sub { $sum += Math::Prime::Util::_XS_lehmer_pi($_) for 
@{$darray[10-2]} },
+});
+print "\n";
 
 sub gendigits {
   my $digits = shift;
diff --git a/examples/sophie_germain.pl b/examples/sophie_germain.pl
index 38dd3bd..4588c7f 100755
--- a/examples/sophie_germain.pl
+++ b/examples/sophie_germain.pl
@@ -2,17 +2,50 @@
 use strict;
 use warnings;
 
-use Math::Prime::Util qw/-nobigint next_prime is_prime/;
+use Math::Prime::Util qw/-nobigint next_prime is_prime prime_iterator/;
 
 my $count = shift || 20;
 
-# Simple little loop looking for Sophie Germain primes (numbers where
-# p and 2p+1 are both prime).  Calculating the first 100k runs 17x faster than
-# Math::NumSeq::SophieGermainPrimes (about 3x faster if the latter's algorithm
-# is changed).
-my $prime = 2;
+# Find Sophie Germain primes (numbers where p and 2p+1 are both prime).
+
+# In this method, we add a filter in front of our iterator, to create a
+# Sophie-Germain-prime iterator.  This isn't the fastest way, but it's still
+# 20x faster than Math::NumSeq::SophieGermainPrimes at 300k.  If we add the
+# two-line precalc shown below, we can get another 4x or so more.
+
+sub get_sophie_germain_iterator {
+  my $p = shift || 2;
+  my $it = prime_iterator($p);
+  return sub {
+    do { $p = $it->() } while !is_prime(2*$p+1);
+    $p;
+  };
+}
+my $sgit = get_sophie_germain_iterator();
 for (1..$count) {
-  $prime = next_prime($prime) while (!is_prime(2*$prime+1));
-  print "$prime\n";
-  $prime = next_prime($prime);
+  print $sgit->(), "\n";
 }
+
+
+# Method 2.  At 300k this is 70x faster than Math::NumSeq::SophieGermainPrimes.
+#
+# my $estimate = 100 + int( nth_prime_upper($count) * 1.6 * log($count) );
+# prime_precalc(2 * $estimate);
+# 
+# my $prime = 2;
+# for (1..$count) {
+#   $prime = next_prime($prime) while (!is_prime(2*$prime+1));
+#   print "$prime\n";
+#   $prime = next_prime($prime);
+# }
+
+# Alternate method, 10-20% faster, would benefit from a tighter estimate.
+#
+# my $numfound = 0;
+# forprimes {
+#   if ($numfound < $count && is_prime(2*$_+1)) {
+#     print "$_\n";
+#     $numfound++;
+#   }
+# } $estimate;
+# die "Estimate too low" unless $numfound >= $count;
diff --git a/examples/twin_primes.pl b/examples/twin_primes.pl
index aa1ec5a..91e5372 100755
--- a/examples/twin_primes.pl
+++ b/examples/twin_primes.pl
@@ -2,20 +2,34 @@
 use strict;
 use warnings;
 
-use Math::Prime::Util qw/-nobigint next_prime is_prime/;
+use Math::Prime::Util qw/prime_iterator nth_prime_upper prime_precalc/;
 
 my $count = shift || 20;
 
-# Simple little loop looking for twin primes (numbers where p and p+2 are
-# both prime).  Print them both.  About 3x faster than 
Math::NumSeq::TwinPrimes.
-my $prime = 2;
-my $next;
+# Find twin primes (numbers where p and p+2 are prime)
+
+# Time for the first 300k:
+#   3m28s  Math::NumSeq::TwinPrimes
+#   2.5s   this iterator
+#   9.1s   this iterator without the precalc
+
+# This speeds things up, but isn't necessary.
+my $estimate = 5000 + int( nth_prime_upper($count) * 1.4 * log($count) );
+prime_precalc($estimate);
+
+sub get_twin_prime_iterator {
+  my $p = shift || 2;
+  my $it = prime_iterator($p);
+  my $prev = $it->();    # prev = 2
+  $p = $it->();          # p = 3
+  return sub {
+    do {
+      ($prev, $p) = ($p, $it->())
+    } while ($p-$prev) != 2;
+    $prev;
+  };
+}
+my $twinit = get_twin_prime_iterator();
 for (1..$count) {
-  while (1) {
-    $next = next_prime($prime);
-    last if ($next-2) == $prime;
-    $prime = $next;
-  }
-  print $prime, ",", $prime+2, "\n";
-  $prime = $next;
+  print $twinit->(), "\n";
 }
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index d9b5d93..d13e3f9 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -3839,6 +3839,21 @@ Examining the η3(x) function of Planat and Solé (2011):
   say prime_count_approx(1000000);
   say nu3(1000000);
 
+Construct and use a Sophie-Germain prime iterator:
+
+  sub make_sophie_germain_iterator {
+    my $p = shift || 2;
+    my $it = prime_iterator($p);
+    return sub {
+      do { $p = $it->() } while !is_prime(2*$p+1);
+      $p;
+    };
+  }
+  my $sgit = make_sophie_germain_iterator();
+  for (1 .. 10000) {
+    print $sgit->(), "\n";
+  }
+
 Project Euler, problem 3 (Largest prime factor):
 
   use Math::Prime::Util qw/factor/;
diff --git a/xt/small-is-next-prev.pl b/xt/small-is-next-prev.pl
index fec2483..5fa87d2 100755
--- a/xt/small-is-next-prev.pl
+++ b/xt/small-is-next-prev.pl
@@ -1,22 +1,73 @@
 #!/usr/bin/env perl
 use strict;
 use warnings;
+use Math::Prime::Util qw/:all/;
+use Time::HiRes qw(gettimeofday tv_interval);
 $| = 1;  # fast pipes
 
-my $limit = shift || 20_000;
+my $mpu_limit = shift || 50_000_000;
+my $mp_limit = shift || 20_000;
 
-use Math::Prime::Util qw/:all/;
-# Use another code base for comparison.
-use Math::Primality;
+# 1. forprimes does a segmented sieve and calls us for each prime.  This is
+#    independent of is_prime and the main sieve.  So for each entry let's
+#    compare next_prime and prev_prime.
+{
+  print "Using MPU forprimes to $mpu_limit\n";
+  my $start_time = [gettimeofday];
+  my $nextprint = 5000000;
+  my $n = 0;
+  forprimes {
+    die "next $n not $_" unless next_prime($n) == $_;
+    die "prev $n" unless prev_prime($_) == $n;
+    $n = $_;
+    if ($n > $nextprint) { print "$n..";  $nextprint += 5000000; }
+  } $mpu_limit;
+  my $seconds = tv_interval($start_time);
+  my $micro_per_call = ($seconds * 1000000) / (2*prime_count($mpu_limit));
+  printf "\nSuccess using forprimes to $mpu_limit.  %6.2f uSec/call\n", 
$micro_per_call;
+}
+
+print "\n";
+
+# 2. Just like before, but now we'll call prime_precalc first.  This makes the
+#    prev_prime and next_prime functions really fast since they just look in
+#    the cached sieve.
+{
+  print "Using MPU forprimes to $mpu_limit with prime_precalc\n";
+  my $start_time = [gettimeofday];
+  prime_precalc($mpu_limit);
+  my $nextprint = 5000000;
+  my $n = 0;
+  forprimes {
+    die "next $n not $_" unless next_prime($n) == $_;
+    die "prev $n" unless prev_prime($_) == $n;
+    $n = $_;
+    if ($n > $nextprint) { print "$n..";  $nextprint += 5000000; }
+  } $mpu_limit;
+  my $seconds = tv_interval($start_time);
+  my $micro_per_call = ($seconds * 1000000) / (2*prime_count($mpu_limit));
+  printf "\nSuccess using forprimes/precalc to $mpu_limit.  %6.2f 
uSec/call\n", $micro_per_call;
+}
+
+print "\n";
 
-foreach my $n (0 .. $limit) {
-  die "next $n" unless next_prime($n) == Math::Primality::next_prime($n);
-  if ($n <= 2) {
-    die "prev $n" unless prev_prime($n) == 0;
-  } else {
-    die "prev $n" unless prev_prime($n) == Math::Primality::prev_prime($n);
+# 3. Now we'll use Math::Primality to compare next_prime, prev_prime, and
+#    is_prime.
+{
+  require Math::Primality;
+  print "Using Math::Primality to $mp_limit\n";
+  my $start_time = [gettimeofday];
+  foreach my $n (0 .. $mp_limit) {
+    die "next $n" unless next_prime($n) == Math::Primality::next_prime($n);
+    if ($n <= 2) {
+      die "prev $n" unless prev_prime($n) == 0;
+    } else {
+      die "prev $n" unless prev_prime($n) == Math::Primality::prev_prime($n);
+    }
+    die "is $n" unless is_prime($n) == Math::Primality::is_prime($n);
+    print "$n.." unless $n % 10000;
   }
-  die "is $n" unless is_prime($n) == Math::Primality::is_prime($n);
-  print "$n.." unless $n % 10000;
+  my $seconds = tv_interval($start_time);
+  my $micro_per_call = ($seconds * 1000000) / (6*prime_count($mp_limit));
+  printf "\nSuccess using Math::Primality to $mp_limit.  %6.2f uSec/call\n", 
$micro_per_call;
 }
-print "Success to $limit!\n";
diff --git a/examples/test-bpsw.pl b/xt/test-bpsw.pl
similarity index 100%
rename from examples/test-bpsw.pl
rename to xt/test-bpsw.pl
diff --git a/examples/test-factor-mpxs.pl b/xt/test-factor-mpxs.pl
similarity index 88%
rename from examples/test-factor-mpxs.pl
rename to xt/test-factor-mpxs.pl
index ab5f0be..331d9ea 100755
--- a/examples/test-factor-mpxs.pl
+++ b/xt/test-factor-mpxs.pl
@@ -33,8 +33,8 @@ print "OK for first 1";
 my $dig = 1;
 my $i = 9;
 foreach my $n (2 .. $nlinear) {
-  my @mfxs = sort { $a<=>$b } prime_factors($n);
-  my @mpu  = sort { $a<=>$b } factor($n);
+  my @mfxs = prime_factors($n);
+  my @mpu  = factor($n);
   die "failure for $n" unless scalar @mfxs == scalar @mpu;
   for (0 .. $#mfxs) { die "failure for $n" unless $mfxs[$_] == $mpu[$_]; }
   if (--$i == 0) {
@@ -48,8 +48,8 @@ print "Testing random numbers from $nlinear to ", $randmax, 
"\n";
 
 while ($nrandom-- > 0) {
   my $n = $nlinear + 1 + $rgen->($randmax - $nlinear);
-  my @mfxs = sort { $a<=>$b } prime_factors($n);
-  my @mpu  = sort { $a<=>$b } factor($n);
+  my @mfxs = prime_factors($n);
+  my @mpu  = factor($n);
   die "failure for $n" unless scalar @mfxs == scalar @mpu;
   for (0 .. $#mfxs) { die "failure for $n" unless $mfxs[$_] == $mpu[$_]; }
   print "." if ($nrandom % 256) == 0;
diff --git a/examples/test-nthapprox.pl b/xt/test-nthapprox.pl
similarity index 100%
rename from examples/test-nthapprox.pl
rename to xt/test-nthapprox.pl
diff --git a/examples/test-pcapprox.pl b/xt/test-pcapprox.pl
similarity index 100%
rename from examples/test-pcapprox.pl
rename to xt/test-pcapprox.pl
diff --git a/examples/test-primes-script.pl b/xt/test-primes-script.pl
similarity index 100%
rename from examples/test-primes-script.pl
rename to xt/test-primes-script.pl
diff --git a/examples/test-primes-script2.pl b/xt/test-primes-script2.pl
similarity index 100%
rename from examples/test-primes-script2.pl
rename to xt/test-primes-script2.pl

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