[libmath-prime-util-perl] 06/13: Minor updates for release
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:47:42 UTC 2015
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 <dana at 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
More information about the Pkg-perl-cvs-commits
mailing list