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

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

commit 30e22b8d7bc52103f8c994d5c33547c4d7b9b601
Author: Dana Jacobsen <d...@acm.org>
Date:   Fri Mar 8 03:51:50 2013 -0800

    Minor updates for release
---
 Changes                    |  6 ++++--
 MANIFEST                   |  2 ++
 README                     |  2 +-
 bin/factor.pl              | 25 +++++++++++++++++++++++--
 factor.c                   |  4 ++--
 lehmer.c                   | 17 +++++++++--------
 lib/Math/Prime/Util.pm     | 18 +++++++++---------
 lib/Math/Prime/Util/PP.pm  | 15 ++++++---------
 xt/rwh_primecount.py       | 20 ++++++++++++++++++++
 xt/rwh_primecount_numpy.py | 18 ++++++++++++++++++
 10 files changed, 94 insertions(+), 33 deletions(-)

diff --git a/Changes b/Changes
index 3261b2a..df6c7f3 100644
--- a/Changes
+++ b/Changes
@@ -2,15 +2,17 @@ Revision history for Perl extension Math::Prime::Util.
 
 0.24  5 March 2013
 
-    - euler_phi on a range wasn't working right with some ranges.
-
     - Fix compilation with old pre-C99 strict compilers (decl after statement).
 
+    - euler_phi on a range wasn't working right with some ranges.
+
     - More XS prime count improvements to speed and space.
 
     - PP prime count for 10^9 and larger is ~2x faster and uses much less
       memory.  Similar impact for nth_prime 10^8 or larger.
 
+    - Let factor.pl accept expressions just like primes.pl.
+
 0.23  5 March 2013
 
     - Replace XS Zeta for x > 5 with series from Cephes.  It is 1 eps more
diff --git a/MANIFEST b/MANIFEST
index 8df46f1..b68f0ea 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -91,4 +91,6 @@ xt/primality-small.pl
 xt/primality-aks.pl
 xt/factor-holf.pl
 xt/make-script-test-data.pl
+xt/rwh_primecount.py
+xt/rwh_primecount_numpy.py
 .travis.yml
diff --git a/README b/README
index 9c5a969..14302c3 100644
--- a/README
+++ b/README
@@ -1,4 +1,4 @@
-Math::Prime::Util version 0.23
+Math::Prime::Util version 0.24
 
 A set of utilities related to prime numbers.  These include multiple sieving
 methods, is_prime, prime_count, nth_prime, approximations and bounds for
diff --git a/bin/factor.pl b/bin/factor.pl
index 2b585a1..6905038 100755
--- a/bin/factor.pl
+++ b/bin/factor.pl
@@ -3,7 +3,7 @@ use strict;
 use warnings;
 use Getopt::Long;
 use bigint try => 'GMP';
-use Math::Prime::Util qw/factor/;
+use Math::Prime::Util qw/factor nth_prime/;
 $| = 1;
 no bigint;
 
@@ -14,7 +14,7 @@ GetOptions(\%opts,
           ) || die_usage();
 if (exists $opts{'version'}) {
   my $version_str =
-   "factor.pl version 1.0 using Math::Prime::Util $Math::Prime::Util::VERSION";
+   "factor.pl version 1.1 using Math::Prime::Util $Math::Prime::Util::VERSION";
   $version_str .= " and MPU::GMP $Math::Prime::Util::GMP::VERSION"
     if Math::Prime::Util::prime_get_config->{'gmp'};
   $version_str .= "\nWritten by Dana Jacobsen.\n";
@@ -24,16 +24,37 @@ die_usage() if exists $opts{'help'};
 
 if (@ARGV) {
   foreach my $n (@ARGV) {
+    $n = eval_expr($n) unless $n =~ /^\d+$/;
     print "$n: ", join(" ", factor($n)), "\n";
   }
 } else {
   while (<>) {
     chomp;
     foreach my $n (split / /) {
+      $n = eval_expr($n) unless $n =~ /^\d+$/;
       print "$n: ", join(" ", factor($n)), "\n";
     }
   }
 }
+
+# This is rather braindead.  We're going to eval their input so they can give
+# arbitrary expressions.  But we only want to allow math-like strings.
+sub eval_expr {
+  my $expr = shift;
+  die "$expr cannot be evaluated" if $expr =~ /:/;  # Use : for escape
+  $expr =~ s/nth_prime\(/:1(/g;
+  $expr =~ s/log\(/:2(/g;
+  die "$expr cannot be evaluated" if $expr =~ tr|-0123456789+*/() :||c;
+  $expr =~ s/:1/nth_prime/g;
+  $expr =~ s/:2/log/g;
+  $expr =~ s/(\d+)/ Math::BigInt->new($1) /g;
+  my $res = eval $expr; ## no critic
+  die "Cannot eval: $expr\n" if !defined $res;
+  $res = int($res->bstr) if ref($res) eq 'Math::BigInt' && $res <= ~0;
+  $res;
+}
+
+
 sub die_usage {
   die <<EOU;
 Usage: $0 [options] [number] ...
diff --git a/factor.c b/factor.c
index 2bd894d..45d74d4 100644
--- a/factor.c
+++ b/factor.c
@@ -493,7 +493,7 @@ int holf_factor(UV n, UV *factors, UV rounds)
   MPUassert( (n >= 3) && ((n%2) != 0) , "bad n in holf_factor");
 
   for (i = 1; i <= rounds; i++) {
-    s = sqrt( (double)n * (double)i );
+    s = (UV) sqrt( (double)n * (double)i );
     /* Assume s^2 isn't a perfect square.  We're rapidly losing precision
      * so we won't be able to accurately detect it anyway. */
     s++;    /* s = ceil(sqrt(n*i)) */
@@ -1006,7 +1006,7 @@ int racing_squfof_factor(UV n, UV *factors, UV rounds)
     mult_save[i].valid = 1;
 
     mult_save[i].b0 = isqrt(nn64);
-    mult_save[i].imax = sqrt(mult_save[i].b0) / 3;
+    mult_save[i].imax = (UV) (sqrt(mult_save[i].b0) / 3);
     if (mult_save[i].imax < 20)     mult_save[i].imax = 20;
     if (mult_save[i].imax > rounds) mult_save[i].imax = rounds;
 
diff --git a/lehmer.c b/lehmer.c
index bf52305..bf972a6 100644
--- a/lehmer.c
+++ b/lehmer.c
@@ -25,17 +25,18 @@
  *
  * The phi(x,a) calculation is unique, to the best of my knowledge.  It uses
  * two lists of all x values + signed counts for the given 'a' value, and walks
- * 'a' down until it is small enough to calculate directly (either with Mapes
- * or using a calculated table using the primorial/totient method).  This
- * is relatively fast and low memory compared to many other solutions.  As with
- * all Lehmer-Meissel-Legendre algorithms, memory use will be a constraint
- * with large values of x (see the table below).
+ * 'a' down until it is small enough to calculate directly using Mapes' method.
+ * This is relatively fast and low memory compared to many other solutions.
+ * As with all Lehmer-Meissel-Legendre algorithms, memory use will be a
+ * constraint with large values of x (see the table below).
  *
  * If you want something better, I highly recommend the paper "Computing
  * Pi(x): the combinatorial method" (2006) by Tomás Oliveira e Silva.  His
- * implementation is certainly much faster and lower memory than this, but I
- * have not seen any working source code for one of the LMO methods so it is
- * difficult to compare.
+ * implementation is certainly much faster and lower memory than this.  I have
+ * briefly run Christian Bau's LMO implementation which has the big advantage
+ * of using almost no memory.  On the same machine that ran the times below,
+ * I got 10^16 in 36s, 10^17 in 165s, 10^18 in 769s, 10^19 in 4162s.  All
+ * using a single core.  That's 10-50x faster than this code.
  *
  * Using my sieve code with everything running in serial, calculating pi(10^12)
  * is done under 1 second on my computer.  pi(10^14) takes under 30 seconds,
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index c072b9e..09ae7a9 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -6,7 +6,7 @@ use Bytes::Random::Secure;
 
 BEGIN {
   $Math::Prime::Util::AUTHORITY = 'cpan:DANAJ';
-  $Math::Prime::Util::VERSION = '0.23';
+  $Math::Prime::Util::VERSION = '0.24';
 }
 
 # parent is cleaner, and in the Perl 5.10.1 / 5.12.0 core, but not earlier.
@@ -1992,7 +1992,7 @@ Math::Prime::Util - Utilities related to prime numbers, 
including fast sieves an
 
 =head1 VERSION
 
-Version 0.23
+Version 0.24
 
 
 =head1 SYNOPSIS
@@ -3349,9 +3349,9 @@ Perl modules, counting the primes to C<800_000_000> (800 
million):
 
   Time (s)   Module                      Version  Notes
   ---------  --------------------------  -------  -----------
-       0.03  Math::Prime::Util           0.12     using Lehmer's method
-       0.28  Math::Prime::Util           0.17     segmented mod-30 sieve
-       0.47  Math::Prime::Util::PP       0.14     Perl (Lehmer's method)
+       0.007 Math::Prime::Util           0.12     using Lehmer's method
+       0.27  Math::Prime::Util           0.17     segmented mod-30 sieve
+       0.39  Math::Prime::Util::PP       0.24     Perl (Lehmer's method)
        0.9   Math::Prime::Util           0.01     mod-30 sieve
        2.9   Math::Prime::FastSieve      0.12     decent odd-number sieve
       11.7   Math::Prime::XS             0.26     needs some optimization
@@ -3363,10 +3363,10 @@ Perl modules, counting the primes to C<800_000_000> 
(800 million):
   ~11000     Math::Primality             0.04     Perl + Math::GMPz
   >20000     Math::Big                   1.12     Perl, > 26GB RAM used
 
-Python can do this in 2.8s using an RWH numpy function, 14.3s using an RWH
-pure Python function.  However the standard modules are far slower.  mpmath
-v0.17 primepi takes 169.5s and 25+ GB of RAM.  sympi 0.7.1 primepi takes
-292.2s.
+Python's standard modules are very slow: mpmath v0.17 primepi takes 169.5s
+and 25+ GB of RAM.  sympy 0.7.1 primepi takes 292.2s.  However there are very
+fast solutions written by Robert William Hanks (included in the xt/
+directory of this distribution): pure Python in 12.1s and numpy in 2.8s.
 
 
 C<is_prime>: my impressions for various sized inputs:
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 5d68b16..1547766 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -5,7 +5,7 @@ use Carp qw/carp croak confess/;
 
 BEGIN {
   $Math::Prime::Util::PP::AUTHORITY = 'cpan:DANAJ';
-  $Math::Prime::Util::PP::VERSION = '0.20';
+  $Math::Prime::Util::PP::VERSION = '0.24';
 }
 
 # The Pure Perl versions of all the Math::Prime::Util routines.
@@ -197,7 +197,7 @@ sub _sieve_erat_string {
     }
     do { $n += 2 } while substr($sieve, $n>>1, 1);
   }
-  \$sieve;
+  return \$sieve;
 }
 
 # TODO: this should be plugged into precalc, memfree, etc. just like the C code
@@ -215,7 +215,7 @@ sub _sieve_erat_string {
       $primary_sieve_ref = _sieve_erat_string($primary_sieve_size);
     }
     my $sieve = substr($$primary_sieve_ref, 0, ($end+1)>>1);
-    \$sieve;
+    return \$sieve;
   }
 }
 
@@ -444,11 +444,8 @@ sub _legendre_phi {
 sub _sieve_prime_count {
   my $high = shift;
   return (0,0,1,2,2,3,3)[$high] if $high < 7;
-  $high-- if ($high % 2) == 0; # Make high go to odd number.
-
-  my $sieveref = _sieve_erat($high);
-  my $count = 1 + $$sieveref =~ tr/0//;
-  return $count;
+  $high-- unless ($high & 1);
+  return 1 + ${_sieve_erat($high)} =~ tr/0//;
 }
 
 sub _count_with_sieve {
@@ -457,7 +454,7 @@ sub _count_with_sieve {
   my $count = 0;
   if   ($low < 3) { $low = 3; $count++; }
   else            { $low |= 1; }
-  $high-- if ($high % 2) == 0; # Make high go to odd number.
+  $high-- unless ($high & 1);
   return $count if $low > $high;
   my $sbeg = $low >> 1;
   my $send = $high >> 1;
diff --git a/xt/rwh_primecount.py b/xt/rwh_primecount.py
new file mode 100755
index 0000000..9afd6c5
--- /dev/null
+++ b/xt/rwh_primecount.py
@@ -0,0 +1,20 @@
+#!/usr/bin/env python
+from math import sqrt, ceil
+
+def rwh_pc(n):
+    # 
http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
+    """ Input n>=6, Returns a list of primes, 2 <= p < n """
+    correction = (n%6>1)
+    n = {0:n,1:n-1,2:n+4,3:n+3,4:n+2,5:n+1}[n%6]
+    sieve = [True] * (n/3)
+    sieve[0] = False
+    for i in xrange(int(n**0.5)/3+1):
+      if sieve[i]:
+        k=3*i+1|1
+        sieve[      ((k*k)/3)      ::2*k]=[False]*((n/6-(k*k)/6-1)/k+1)
+        
sieve[(k*k+4*k-2*k*(i&1))/3::2*k]=[False]*((n/6-(k*k+4*k-2*k*(i&1))/6-1)/k+1)
+    sieve[n/3-correction] = False
+    return 2 + sum(sieve)
+    #return [2,3] + [3*i+1|1 for i in xrange(1,n/3-correction) if sieve[i]]
+
+print rwh_pc(800000000)
diff --git a/xt/rwh_primecount_numpy.py b/xt/rwh_primecount_numpy.py
new file mode 100755
index 0000000..77e2987
--- /dev/null
+++ b/xt/rwh_primecount_numpy.py
@@ -0,0 +1,18 @@
+#!/usr/bin/env python
+#from math import sqrt, ceil
+import numpy as np
+
+def rwh_pcn(n):
+    # 
http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
+    """ Input n>=6, Returns a list of primes, 2 <= p < n """
+    sieve = np.ones(n/3 + (n%6==2), dtype=np.bool)
+    for i in xrange(1,int(n**0.5)/3+1):
+        if sieve[i]:
+            k=3*i+1|1
+            sieve[       k*k/3     ::2*k] = False
+            sieve[k*(k-2*(i&1)+4)/3::2*k] = False
+    return 1 + np.count_nonzero(sieve)
+    #return np.r_[2,3,((3*np.nonzero(sieve)[0]+1)|1)]
+
+
+print rwh_pcn(800000000)

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