# [libmath-prime-util-perl] 25/33: Update some examples

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

ppm-guest pushed a commit to annotated tag v0.37
in repository libmath-prime-util-perl.```
```
commit 1f31843759f3d21a0ff9926d2bafcf1eae4dfac0
Author: Dana Jacobsen <d...@acm.org>
Date:   Fri Jan 24 17:04:58 2014 -0800

Update some examples
---
examples/abundant.pl       |   2 +-
examples/sophie_germain.pl | 133 ++++++++++++++++++++++++++++-----------------
examples/twin_primes.pl    |  30 ++++++----
3 files changed, 103 insertions(+), 62 deletions(-)

diff --git a/examples/abundant.pl b/examples/abundant.pl
index 23b3564..8cfe2b4 100755
--- a/examples/abundant.pl
+++ b/examples/abundant.pl
@@ -5,7 +5,6 @@ use warnings;
# Find the first N abundant, deficient, or perfect numbers.

use Math::Prime::Util qw/divisor_sum next_prime is_prime/;
-use Math::BigInt try => "GMP,Pari";

my \$count = shift || 20;
my \$type = lc(shift || 'abundant');
@@ -26,6 +25,7 @@ if (\$type eq 'abundant') {
# We just look for 2^(p-1)*(2^p-1) where 2^p-1 is prime.
# Basically we're just finding Mersenne primes.
# It's possible there are odd perfect numbers larger than 10^1500.
+  do { require Math::BigInt;  Math::BigInt->import(try=>"GMP,Pari"); };
while (\$count-- > 0) {
while (1) {
\$p = next_prime(\$p);
diff --git a/examples/sophie_germain.pl b/examples/sophie_germain.pl
--- a/examples/sophie_germain.pl
+++ b/examples/sophie_germain.pl
@@ -6,59 +6,92 @@ use Math::Prime::Util qw/prime_iterator is_prime
next_prime nth_prime_upper prime_precalc forprimes/;

my \$count = shift || 20;
+my \$method = shift || 'forprimes';
+my \$precalc = 0;   # If set, precalc all the values we'll call is_prime on

# 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 more.
-#
-# Example:
-#  time perl examples/sophie_germain.pl 300000 | md5sum
-#    d380d31256cc9bc54eb5f236b3edc16d  -
-#    9.673s
-#
-#  time perl -MMath::NumSeq::SophieGermainPrimes -E 'my \$seq =
Math::NumSeq::SophieGermainPrimes->new; do { say 0+(\$seq->next)[1] } for
1..300000' | md5sum
-#    d380d31256cc9bc54eb5f236b3edc16d  -
-#    4m11.5s
+# Four methods are shown: forprimes, iter, iter2, and MNS.
+
+# Times for 300k:
#
-# With method 2:
-#  time perl examples/sophie_germain.pl 300000 | md5sum
-#    d380d31256cc9bc54eb5f236b3edc16d  -
-#    1.828s
-
-
-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;
-  };
+#                300k               1M
+# precalc:
+#   forprimes    1.3s   9.0MB       7.1s   21.6MB
+#   iter         2.8s   8.7MB      12.6s   21.4MB
+#   iter2        1.9s   8.7MB       9.4s   21.4MB
+# no precalc:
+#   forprimes    1.5s   4.5MB       5.6s   4.5MB
+#   iter         9.5s   4.3MB      37.5s   4.3MB
+#   iter2        8.5s   4.3MB      33.9s   4.3MB
+#   MNS        254.3s  11.3MB   >1500s   >15 MB
+
+if (\$precalc) {
+  prime_precalc(2 * sg_upper_bound(\$count));
}
-my \$sgit = get_sophie_germain_iterator();
-print \$sgit->(), "\n" for 1 .. \$count;

-# 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;
+
+if (\$method eq 'forprimes') {
+
+  my \$estimate = sg_upper_bound(\$count);
+  my \$numfound = 0;
+  forprimes {
+    if (\$numfound < \$count && is_prime(2*\$_+1)) {
+      print "\$_\n";
+      \$numfound++;
+    }
+  } \$estimate;
+  die "Estimate too low" unless \$numfound >= \$count;
+
+} elsif (\$method eq 'iter') {
+
+  # Wrap the standard iterator
+  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();
+  print \$sgit->(), "\n" for 1 .. \$count;
+
+} elsif (\$method eq 'iter2') {
+
+  # Iterate directly using next_prime
+  my \$prime = 2;
+  for (1 .. \$count) {
+    \$prime = next_prime(\$prime) while !is_prime(2*\$prime+1);
+    print "\$prime\n";
+    \$prime = next_prime(\$prime);
+  }
+
+} elsif (\$method eq 'MNS') {
+
+  # Use Math::NumSeq
+  require Math::NumSeq::SophieGermainPrimes;
+  my \$seq = Math::NumSeq::SophieGermainPrimes->new;
+  for (1 .. \$count) {
+    print 0+(\$seq->next)[1];
+  }
+
+}
+
+# Used for precalc and the forprimes example
+sub sg_upper_bound {
+  my \$count = shift;
+  my \$nth = nth_prime_upper(\$count);
+  # For lack of a better formula, do this step-wise estimate.
+  my \$estimate = (\$count <   5000) ? 150 + int( \$nth * log(\$nth) * 1.2 )
+               : (\$count <  19000) ? int( \$nth * log(\$nth) * 1.135 )
+               : (\$count <  45000) ? int( \$nth * log(\$nth) * 1.10 )
+               : (\$count < 100000) ? int( \$nth * log(\$nth) * 1.08 )
+               : (\$count < 165000) ? int( \$nth * log(\$nth) * 1.06 )
+               : (\$count < 360000) ? int( \$nth * log(\$nth) * 1.05 )
+               : (\$count < 750000) ? int( \$nth * log(\$nth) * 1.04 )
+               : (\$count <1700000) ? int( \$nth * log(\$nth) * 1.03 )
+               :                     int( \$nth * log(\$nth) * 1.02 );
+
+  return \$estimate;
+}
diff --git a/examples/twin_primes.pl b/examples/twin_primes.pl
index 362a737..514cc70 100755
--- a/examples/twin_primes.pl
+++ b/examples/twin_primes.pl
@@ -2,8 +2,7 @@
use strict;
use warnings;

-use Math::Prime::Util qw/-nobigint
-                         prime_iterator  prime_iterator_object
+use Math::Prime::Util qw/prime_iterator  prime_iterator_object
next_prime  is_prime
nth_prime_upper  prime_precalc/;

@@ -12,16 +11,25 @@ my \$count = shift || 20;
# Find twin primes (numbers where p and p+2 are prime)

# Time for the first 300k:
+#
+# Not iterators:
+#   0.6s   forprimes { say \$l if \$l+2==\$_; \$l=\$_; } 64764841
#   1.0s   bin/primes.pl --twin 2 64764839
-#   1.4s   get_twin_prime_iterator2
-#   2.3s   get_twin_prime_iterator1
-#   4.1s   get_twin_prime_iterator3
-#   6.1s   get_twin_prime_iterator3 (object iterator)
-#   7.6s   get_twin_prime_iterator2 without precalc
-#   8.4s   get_twin_prime_iterator1 without precalc
-#  10.9s   get_twin_prime_iterator3 without precalc
-#  13.1s   get_twin_prime_iterator3 without precalc (object iterator)
-# 219.8s   Math::NumSeq::TwinPrimes (Perl 5.19.4 with v66)
+#
+# Iterators with precalc:
+#   1.6s   get_twin_prime_iterator2
+#   2.4s   get_twin_prime_iterator1
+#   4.2s   get_twin_prime_iterator3
+#   4.5s   get_twin_prime_iterator4 (object iterator)
+#
+# Iterators without precalc:
+#   7.7s   get_twin_prime_iterator2
+#   8.5s   get_twin_prime_iterator1
+#  10.8s   get_twin_prime_iterator3
+#  16.7s   get_twin_prime_iterator4 (object iterator)
+#
+# Alternatives:
+# 251.9s   Math::NumSeq::TwinPrimes (Perl 5.19.7, Math::NumSeq 67)

# This speeds things up, but isn't necessary.
my \$estimate = 5000 + int( nth_prime_upper(\$count) * 1.4 * log(\$count) );

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