[libmath-prime-util-perl] 30/43: Inverse Riemann R for nth_prime_approx

[Partha P. Mukherjee]

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

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

commit 51eb53406bf6c7a9418223b4175b5e64ac729e5c
Author: Dana Jacobsen <dana at acm.org>
Date:   Sun Apr 6 23:20:18 2014 -0700

    Inverse Riemann R for nth_prime_approx
---
 Changes                   |  4 ++++
 lib/Math/Prime/Util.pm    |  8 +++++---
 lib/Math/Prime/Util/PP.pm |  5 ++---
 util.c                    | 45 ++++++++++++++-------------------------------
 4 files changed, 25 insertions(+), 37 deletions(-)

diff --git a/Changes b/Changes
index 713c556..00611f2 100644
--- a/Changes
+++ b/Changes
@@ -16,6 +16,10 @@ Revision history for Perl module Math::Prime::Util
 
     - Speed up exp_mangoldt.
 
+    - nth_prime_approx uses inverse RiemannR in XS code for better accuracy.
+      Cippola 1902 is still used for PP and large values, with a slightly
+      more accurate third order correction.
+
 
 0.39  2014-03-01
 
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index c083865..d7a01af 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -1333,7 +1333,7 @@ While this method is thousands of times faster than generating primes, and
 doesn't involve big tables of precomputed values, it still can take a fair
 amount of time for large inputs.  Calculating the C<10^12th> prime takes
 about 1 second, the C<10^13th> prime takes under 10 seconds, and the
-C<10^14th> prime (3475385758524527) takes under one minute.  Think about
+C<10^14th> prime (3475385758524527) takes under 30 seconds.  Think about
 whether a bound or approximation would be acceptable, as they can be
 computed analytically.
 
@@ -1365,8 +1365,10 @@ for small C<n>.
   say "The one trillionth prime is ~ ", nth_prime_approx(10**12);
 
 Returns an approximation to the C<nth_prime> function, without having to
-generate any primes.  Uses the Cipolla 1902 approximation with two
-polynomials, plus a correction for small values to reduce the error.
+generate any primes.  For values where the nth prime is smaller than
+C<2^64>, an inverse Riemann R function is used.  For larger values, uses the
+Cipolla 1902 approximation with up to 2nd order terms, plus a third order
+correction.
 
 
 =head2 nth_twin_prime
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index f25fa77..bc3f0bb 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -1159,9 +1159,8 @@ sub nth_prime_approx {
   elsif ($n <      12000) { $approx +=  3.0 * $order; }
   elsif ($n <     150000) { $approx +=  2.1 * $order; }
   elsif ($n <  200000000) { $approx +=  0.0 * $order; }
-  else                    { $approx += -0.010 * $order; }
-  # $approx = -0.025 is better for the last, but it gives problems with some
-  # other code that always wants the asymptotic approximation to be >= actual.
+  else                    { $approx += -0.025 * $order; }
+  # If we want the asymptotic approximation to be >= actual, use -0.010.
 
   return int($approx + 0.5);
 }
diff --git a/util.c b/util.c
index 6a46b47..88a5efa 100644
--- a/util.c
+++ b/util.c
@@ -736,41 +736,24 @@ UV nth_prime_lower(UV n)
 
 UV nth_prime_approx(UV n)
 {
-  long double fn, flogn, flog2n, approx, order;
+  long double fn, flogn;
+  UV lo, hi;
 
   if (n < NPRIMES_SMALL)
     return primes_small[n];
 
-  fn     = (long double) n;
-  flogn  = logl(n);
-  flog2n = logl(flogn);    /* Note distinction between log_2(n) and log^2(n) */
-
-  /* Cipolla 1902:
-   *    m=0   fn * ( flogn + flog2n - 1 );
-   *    m=1   + ((flog2n - 2)/flogn) );
-   *    m=2   - (((flog2n*flog2n) - 6*flog2n + 11) / (2*flogn*flogn))
-   *    + O((flog2n/flogn)^3)
-   */
-
-  approx = fn * (  flogn + flog2n - 1.0
-                 + ((flog2n - 2.0) / flogn)
-                 - (((flog2n*flog2n) - 6.0*flog2n + 11.0) / (2*flogn*flogn))
-                );
-
-  /* Apply a correction */
-  order = flog2n / flogn;
-  order = order * order * order * fn;
-  if      (n <      259) { approx += 10.4  * order; }
-  else if (n <      775) { approx +=  7.52 * order; }
-  else if (n <     1271) { approx +=  5.6  * order; }
-  else if (n <     2000) { approx +=  5.2  * order; }
-  else if (n <     4000) { approx +=  4.3  * order; }
-  else if (n <    12000) { approx +=  3.0  * order; }
-  else if (n <   150000) { approx +=  2.1  * order; }
-  else if (n <200000000) {                          }
-  else                   { approx += -0.01 * order; } /* -0.25 is closer */
-
-  return (UV) floorl(approx + 0.5);
+  /* Binary search for inverse Riemann R */
+  fn    = (long double) n;
+  flogn = logl(n);
+  lo    = (UV) (fn * flogn);
+  hi    = (UV) (fn * flogn * 2 + 2);
+  if (hi <= lo) hi = UV_MAX;
+  while (lo < hi) {
+    UV mid = lo + (hi-lo)/2;
+    if (_XS_RiemannR(mid) < fn) lo = mid+1;
+    else                        hi = mid;
+  }
+  return lo;
 }
 
 

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git