[libmath-prime-util-perl] 12/15: Update some benchmarks and examples
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:48:47 UTC 2015
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 <dana at 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
More information about the Pkg-perl-cvs-commits
mailing list