[libmath-prime-util-perl] 14/15: Fix prime_count_approx being really slow for > 10^36 without Math::MPFR
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 c22623dcb469968d13b221a4d088d44bc0336788
Author: Dana Jacobsen <dana at acm.org>
Date: Thu May 30 23:31:30 2013 -0700
Fix prime_count_approx being really slow for > 10^36 without Math::MPFR
---
Changes | 6 +++++-
TODO | 5 +++--
factor.c | 28 +++++++++++++++-------------
factor.h | 2 +-
lib/Math/Prime/Util.pm | 34 +++++++++++++++++++++++-----------
lib/Math/Prime/Util/PP.pm | 42 +++++++++++++++---------------------------
6 files changed, 62 insertions(+), 55 deletions(-)
diff --git a/Changes b/Changes
index 79942d7..aa6c381 100644
--- a/Changes
+++ b/Changes
@@ -1,6 +1,6 @@
Revision history for Perl extension Math::Prime::Util.
-0.29 xx May 2013
+0.29 30 May 2013
- Fix a signed vs. unsigned char issue in ranged moebius. Thanks to the
Debian testers for finding this.
@@ -13,6 +13,10 @@ Revision history for Perl extension Math::Prime::Util.
- forprimes now uses a segmented sieve. This (1) allows arbitrary 64-bit
ranges with good memory use, and (2) allows nesting on threaded perls.
+ - prime_count_approx for very large values (> 10^36) was very slow without
+ Math::MPFR. Switch to Li+correction for large values if Math::MPFR is
+ not available.
+
- Added:
is_pseudoprime (Fermat probable prime test)
is_lucas_pseudoprime (standard Lucas-Selfridge test)
diff --git a/TODO b/TODO
index 8b364b6..33f6bc5 100644
--- a/TODO
+++ b/TODO
@@ -44,6 +44,7 @@
- Figure out a way to make the internal FOR_EACH_PRIME macros use a segmented
sieve.
-- prime_count_approx on a 100-digit bigint is really slow without MPFR.
-
- Refactor where functions exist in .c. Lots of primality tests in factor.c.
+
+- Write a standalone function that demonstrates the memory leak with MULTICALL,
+ so we can use MULTICALL.
diff --git a/factor.c b/factor.c
index 976d098..1dfcc4d 100644
--- a/factor.c
+++ b/factor.c
@@ -380,7 +380,7 @@ int _SPRP2(UV n)
UV const nm1 = n-1;
UV d = n-1;
UV x;
- int b, r, s = 0;
+ int r, s = 0;
MPUassert(n > 3, "S-PRP-2 called with n <= 3");
if (!(n & 1)) return 0;
@@ -502,7 +502,7 @@ int _XS_is_prob_prime(UV n)
*/
int _XS_is_extra_strong_lucas_pseudoprime(UV n)
{
- UV P, D, Q, U, V, t, t2, d, s, b;
+ UV P, D, Q, U, V, d, s, b;
int const _verbose = _XS_get_verbose();
if (n == 2 || n == 3) return 1;
@@ -537,7 +537,7 @@ int _XS_is_extra_strong_lucas_pseudoprime(UV n)
b--;
if (_verbose>3) printf("U2k=%lu V2k=%lu\n", U, V);
if ( (d >> (b-1)) & UVCONST(1) ) {
- t2 = mulmod(U, D, n);
+ UV t2 = mulmod(U, D, n);
U = muladdmod(U, P, V, n);
if (U & 1) { U = (n>>1) + (U>>1) + 1; } else { U >>= 1; }
V = muladdmod(V, P, t2, n);
@@ -1198,7 +1198,8 @@ int racing_squfof_factor(UV n, UV *factors, UV rounds)
/*
*
* The Frobenius-Underwood test has no known counterexamples below 10^13, but
- * has not been extensively tested above that.
+ * has not been extensively tested above that. This is the Minimal Lambda+2
+ * test from section 9 of "Quadratic Composite Tests" by Paul Underwood.
*
* Given the script:
* time mpu 'forprimes { Math::Prime::Util::_XS_is_frobenius_underwood_pseudoprime($_); Math::Prime::Util::_XS_is_frobenius_underwood_pseudoprime($_+2); } 100_000_000'
@@ -1234,17 +1235,18 @@ int racing_squfof_factor(UV n, UV *factors, UV rounds)
*
* We try to structure the primality test like:
* 1) simple divisibility very fast primes and ~10% of composites
- * 2) M-R with base 2 1 Selfridge most remaining composites gone
+ * 2) M-R with base 2 1 Selfridge primes and .00000000002% comps
* 3) Lucas test 2 Selfridge only primes
*
- * Hence given a composite, it will typically cost 0-1 Selfridges, and for
- * our 64-bit values, the final Lucas test has no false positives. Replacing
- * the Lucas test with the F-U test won't save any time. Replacing the whole
- * thing with the F-U test (assuming it has no false results for all 64-bit
- * values), doesn't help much either -- it's 2/3 the cost for primes, but much
- * more expensive for composites. It seems of interest for > 2^64 as a
- * different test to do in addition to BPSW.
- *
+ * Hence given a random composite, about 90% of the time it costs us almost
+ * nothing. After spending 1 Selfridge on the first MR test, less than 32M
+ * composites remain undecided out of 18 quintillion 64-bit composites. The
+ * final Lucas test has no false positives.
+ * Replacing the Lucas test with the F-U test won't save any time. Replacing
+ * the whole thing with the F-U test (assuming it has no false results for
+ * all 64-bit values, which hasn't been verified), doesn't help either.
+ * It's 2/3 the cost for primes, but much more expensive for composites. It
+ * seems of interest for > 2^64 as a different test to do in addition to BPSW.
*/
diff --git a/factor.h b/factor.h
index 29722bf..d017d49 100644
--- a/factor.h
+++ b/factor.h
@@ -22,7 +22,7 @@ extern int _XS_is_pseudoprime(UV n, UV a);
extern int _XS_miller_rabin(UV n, const UV *bases, int nbases);
extern int _SPRP2(UV n);
extern int _XS_is_prob_prime(UV n);
-extern int _XS_is_strong_lucas_pseudoprime(UV n);
+extern int _XS_is_extra_strong_lucas_pseudoprime(UV n);
extern int _XS_is_frobenius_underwood_pseudoprime(UV n);
extern UV _XS_divisor_sum(UV n);
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index d13e3f9..d3c1c98 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -2087,16 +2087,23 @@ sub prime_count_approx {
# Also consider: http://trac.sagemath.org/sage_trac/ticket/8135
# my $result = int( (prime_count_upper($x) + prime_count_lower($x)) / 2);
-
# my $result = int( LogarithmicIntegral($x) );
-
# my $result = int(LogarithmicIntegral($x) - LogarithmicIntegral(sqrt($x))/2);
+ # my $result = RiemannR($x) + 0.5;
- if (ref($x) eq 'Math::BigFloat') {
- # Make sure we get enough accuracy, and also not too much more than needed
- $x->accuracy(length($x->bfloor->bstr())+2);
+ # Sadly my Perl RiemannR function is really slow for big values. If MPFR
+ # is available, then use it -- it rocks. Otherwise, switch to LiCorr for
+ # very big values. This is hacky and shouldn't be necessary.
+ my $result;
+ if ( $x < 1e36 || Math::Prime::Util::PP::_MPFR_available() ) {
+ if (ref($x) eq 'Math::BigFloat') {
+ # Make sure we get enough accuracy, and also not too much more than needed
+ $x->accuracy(length($x->bfloor->bstr())+2);
+ }
+ $result = RiemannR($x) + 0.5;
+ } else {
+ $result = int(LogarithmicIntegral($x) - LogarithmicIntegral(sqrt($x))/2);
}
- my $result = RiemannR($x) + 0.5;
return Math::BigInt->new($result->bfloor->bstr()) if ref($result) eq 'Math::BigFloat';
return int($result);
@@ -2815,11 +2822,16 @@ the Schoenfeld (1976) bounds are used for large values.
" primes below one quintillion.\n";
Returns an approximation to the C<prime_count> function, without having to
-generate any primes. The current implementation uses the Riemann R function
-which is quite accurate: an error of less than C<0.0005%> is typical for
-input values over C<2^32>. A slightly faster
-(0.1 millisecond versus 1 millisecond) but much less accurate answer
-can be obtained by averaging the upper and lower bounds.
+generate any primes. For values under C<10^36> this uses the Riemann R
+function, which is quite accurate: an error of less than C<0.0005%> is typical
+for input values over C<2^32>, and decreases as the input gets larger. If
+L<Math::MPFR> is installed, the Riemann R function is used for all values, and
+will be very fast. If not, then values of C<10^36> and larger will use the
+approximation C<li(x) - li(sqrt(x))/2>. While not as accurate as the Riemann
+R function, it still should have error less than C<0.00000000000000001%>.
+
+A slightly faster but much less accurate answer can be obtained by averaging
+the upper and lower bounds.
=head2 nth_prime
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 530e211..07ef07e 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -43,7 +43,17 @@ our $_Infinity = 0+'inf';
$_Infinity = 20**20**20 if 65535 > $_Infinity; # E.g. Windows
our $_Neg_Infinity = -$_Infinity;
-my $_have_MPFR = -1;
+{
+ my $_have_MPFR = -1;
+ sub _MPFR_available {
+ if ($_have_MPFR < 0) {
+ $_have_MPFR = 0;
+ $_have_MPFR = 1 if (!defined $ENV{MPU_NO_MPFR} || $ENV{MPU_NO_MPFR} != 1)
+ && eval { require Math::MPFR; $Math::MPFR::VERSION>=2.03; };
+ }
+ return $_have_MPFR;
+ }
+}
my $_precalc_size = 0;
sub prime_precalc {
@@ -2121,14 +2131,8 @@ sub ExponentialIntegral {
return 0 if $x == $_Neg_Infinity;
return $_Infinity if $x == $_Infinity;
- # Use MPFR if possible.
- if ($_have_MPFR < 0) {
- $_have_MPFR = 0;
- $_have_MPFR = 1 if (!defined $ENV{MPU_NO_MPFR} || $ENV{MPU_NO_MPFR} != 1)
- && eval { require Math::MPFR; $Math::MPFR::VERSION>=2.03; };
- }
# Gotcha -- MPFR decided to make negative inputs return NaN. Grrr.
- if ($_have_MPFR && $x > 0) {
+ if ($x > 0 && _MPFR_available()) {
my $wantbf = 0;
my $xdigits = 17;
if (defined $bignum::VERSION || ref($x) =~ /^Math::Big/) {
@@ -2231,14 +2235,8 @@ sub LogarithmicIntegral {
return $_Infinity if $x == $_Infinity;
croak "Invalid input to LogarithmicIntegral: x must be > 0" if $x <= 0;
- # Use MPFR if possible.
- if ($_have_MPFR < 0) {
- $_have_MPFR = 0;
- $_have_MPFR = 1 if (!defined $ENV{MPU_NO_MPFR} || $ENV{MPU_NO_MPFR} != 1)
- && eval { require Math::MPFR; $Math::MPFR::VERSION>=2.03; };
- }
# Remember MPFR eint doesn't handle negative inputs
- if ($_have_MPFR && $x >= 1) {
+ if ($x >= 1 && _MPFR_available()) {
my $wantbf = 0;
my $xdigits = 17;
if (defined $bignum::VERSION || ref($x) =~ /^Math::Big/) {
@@ -2365,12 +2363,7 @@ sub RiemannZeta {
my($x) = @_;
# Use MPFR if possible.
- if ($_have_MPFR < 0) {
- $_have_MPFR = 0;
- $_have_MPFR = 1 if (!defined $ENV{MPU_NO_MPFR} || $ENV{MPU_NO_MPFR} != 1)
- && eval { require Math::MPFR; $Math::MPFR::VERSION>=2.03; };
- }
- if ($_have_MPFR) {
+ if (_MPFR_available()) {
my $wantbf = 0;
my $xdigits = 17;
if (defined $bignum::VERSION || ref($x) =~ /^Math::Big/) {
@@ -2456,12 +2449,7 @@ sub RiemannR {
croak "Invalid input to ReimannR: x must be > 0" if $x <= 0;
# Use MPFR if possible.
- if ($_have_MPFR < 0) {
- $_have_MPFR = 0;
- $_have_MPFR = 1 if (!defined $ENV{MPU_NO_MPFR} || $ENV{MPU_NO_MPFR} != 1)
- && eval { require Math::MPFR; $Math::MPFR::VERSION>=2.03; };
- }
- if ($_have_MPFR) {
+ if (_MPFR_available()) {
my $wantbf = 0;
my $xdigits = 17;
if (defined $bignum::VERSION || ref($x) =~ /^Math::Big/) {
--
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