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 <d...@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

_______________________________________________
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