[libmath-prime-util-perl] 18/25: Change nth_prime to use inverse Li+Corr instead of inverse R
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:50:39 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.35
in repository libmath-prime-util-perl.
commit dd2d4d7652b71a5eba74f5cab21fd5355e0c7664
Author: Dana Jacobsen <dana at acm.org>
Date: Fri Dec 6 12:54:42 2013 -0800
Change nth_prime to use inverse Li+Corr instead of inverse R
---
Changes | 9 ++++--
t/14-nthprime.t | 8 +++++
util.c | 93 +++++++++++++++++++++++++++------------------------------
util.h | 1 +
4 files changed, 60 insertions(+), 51 deletions(-)
diff --git a/Changes b/Changes
index 0a5732e..12dcce7 100644
--- a/Changes
+++ b/Changes
@@ -1,6 +1,6 @@
Revision history for Perl module Math::Prime::Util
-0.35 2013-11
+0.35 2013-12
[API Changes]
@@ -19,7 +19,8 @@ Revision history for Perl module Math::Prime::Util
[FUNCTIONALITY AND PERFORMANCE]
- Switched to extended LMO algorithm for prime_count. Much better memory
- use and much faster for large values. Speeds up nth_prime also.
+ use and much faster for large values. Speeds up nth_prime also. Huge
+ thanks to Christian Bau's excellent paper and examples.
- Some fixes for 32-bit.
@@ -27,6 +28,10 @@ Revision history for Perl module Math::Prime::Util
- Fixed a problem with Lehmer prime_count introduced in 0.34.
+ - nth_prime changed from RiemannR to inverse Li (with partial addition).
+ This makes some of the big nth_prime calculations (e.g. 10^15, 10^16)
+ run quite a bit faster as they sieve less on average.
+
0.34 2013-11-19
diff --git a/t/14-nthprime.t b/t/14-nthprime.t
index 74fa023..11632bb 100644
--- a/t/14-nthprime.t
+++ b/t/14-nthprime.t
@@ -35,6 +35,14 @@ my %nthprimes32 = (
1000000 => 15485863,
10000000 => 179424673,
100000000 => 2038074743,
+ # Some values that estimate right around the value
+ 6305537 => 110040379,
+ 6305538 => 110040383,
+ 6305539 => 110040391,
+ 6305540 => 110040407,
+ 6305541 => 110040467,
+ 6305542 => 110040499,
+ 6305543 => 110040503,
);
my %nthprimes64 = (
1000000000 => 22801763489,
diff --git a/util.c b/util.c
index 7a62714..f4ca5ae 100644
--- a/util.c
+++ b/util.c
@@ -650,12 +650,8 @@ UV _XS_nth_prime(UV n)
/* For relatively small values, generate a sieve and count the results.
*
- * For larger values, compute a lower bound, use Lehmer's algorithm to get
- * a fast prime count, then start segment sieving from there.
- *
- * For very large values, binary search on Riemann's R function to get a
- * good approximation, use Lehmer's algorithm to get the count, then walk
- * backwards or sieve forwards.
+ * For larger values, compute an approximate low estimate, use our fast
+ * prime count, then segment sieve forwards or backwards for the rest.
*/
if (upper_limit <= get_prime_cache(0, 0) || upper_limit <= 32*1024*30) {
/* Generate a sieve and count. */
@@ -665,50 +661,30 @@ UV _XS_nth_prime(UV n)
count += count_segment_maxcount(cache_sieve, segment_size, target, &p);
release_prime_cache(cache_sieve);
} else {
- double fn = n;
- double flogn = log(fn);
- double flog2n = log(flogn); /* Dusart 2010, page 2, n >= 3 */
- UV lower_limit = fn * (flogn + flog2n - 1.0 + ((flog2n-2.10)/flogn));
-#if BITS_PER_WORD == 32
- if (1) {
-#else
- if (n <= UVCONST(20000000000)) {
-#endif
- /* Calculate lower limit, get count, sieve to that */
- segment_size = lower_limit / 30;
- lower_limit = 30 * segment_size - 1;
- count = _XS_LMO_pi(lower_limit) - 3;
- MPUassert(count <= target, "Pi(nth_prime_lower(n))) > n");
- } else {
- /* Compute approximate nth prime via binary search on R(n) */
- UV lo = lower_limit;
- UV hi = upper_limit;
- double lor = _XS_RiemannR(lo);
- double hir = _XS_RiemannR(hi);
- while (lor < hir) {
- UV mid = (UV) ((lo + hi) / 2);
- double midr = _XS_RiemannR(mid);
- if (midr <= n) { lo = mid+1; lor = _XS_RiemannR(lo); }
- else { hi = mid; hir = midr; }
- }
- /* Bias toward lower, because we want to sieve up if possible */
- lower_limit = (UV) (double)(0.9999999*(lo-1));
- segment_size = lower_limit / 30;
- lower_limit = 30 * segment_size - 1;
- count = _XS_LMO_pi(lower_limit);
- /*
- printf("We've estimated %lu too %s.\n", (count>n)?count-n:n-count, (count>n)?"FAR":"little");
- printf("Our limit %lu %s a prime\n", lower_limit, _XS_is_prime(lower_limit) ? "is" : "is not");
- */
-
- if (count > n) { /* Too far. Walk backwards */
- if (_XS_is_prime(lower_limit)) count--;
- for (p = 0; p <= (count-n); p++)
- lower_limit = _XS_prev_prime(lower_limit);
- return lower_limit;
- }
- count -= 3;
+ /* A binary search on RiemannR is nice, but ends up either often being
+ * being higher (requiring going backwards) or biased and then far too
+ * low. Using the inverse Li is easier and more consistent. */
+ UV lower_limit = _XS_Inverse_Li(n);
+ /* For even better performance, add in half the usual correction, which
+ * will get us even closer, so even less sieving required. However, it
+ * is now possible to get a result higher than the value, so we'll need
+ * to handle that case. It still ends up being a better deal than R,
+ * given that we don't have a fast backward sieve. */
+ lower_limit += _XS_Inverse_Li(isqrt(n))/4;
+ segment_size = lower_limit / 30;
+ lower_limit = 30 * segment_size - 1;
+ count = _XS_LMO_pi(lower_limit);
+
+ /* printf("We've estimated %lu too %s.\n", (count>n)?count-n:n-count, (count>n)?"FAR":"little"); */
+ /* printf("Our limit %lu %s a prime\n", lower_limit, _XS_is_prime(lower_limit) ? "is" : "is not"); */
+
+ if (count >= n) { /* Too far. Walk backwards */
+ if (_XS_is_prime(lower_limit)) count--;
+ for (p = 0; p <= (count-n); p++)
+ lower_limit = _XS_prev_prime(lower_limit);
+ return lower_limit;
}
+ count -= 3;
/* Make sure the segment siever won't have to keep resieving. */
prime_precalc(isqrt(upper_limit));
@@ -1110,6 +1086,25 @@ double _XS_LogarithmicIntegral(double x) {
return _XS_ExponentialIntegral(log(x));
}
+/* Thanks to Kim Walisch for this idea */
+UV _XS_Inverse_Li(UV x) {
+ double n = x;
+ double logn = log(n);
+ UV lo = (UV) (n*logn);
+ UV hi = (UV) (n*logn * 2 + 2);
+
+ if (x < 1) return 0;
+ if (hi <= lo) hi = UV_MAX;
+ while (lo < hi) {
+ UV mid = lo + (hi-lo)/2;
+ if (_XS_LogarithmicIntegral(mid) < x) lo = mid+1;
+ else hi = mid;
+ }
+ return lo;
+}
+
+
+
/*
* Storing the first 10-20 Zeta values makes sense. Past that it is purely
* to avoid making the call to generate them ourselves. We could cache the
diff --git a/util.h b/util.h
index 15fd457..22f19b1 100644
--- a/util.h
+++ b/util.h
@@ -25,6 +25,7 @@ extern double _XS_ExponentialIntegral(double x);
extern double _XS_LogarithmicIntegral(double x);
extern long double ld_riemann_zeta(long double x);
extern double _XS_RiemannR(double x);
+extern UV _XS_Inverse_Li(UV x);
/* Above this value, is_prime will do deterministic Miller-Rabin */
/* With 64-bit math, we can do much faster mulmods from 2^16-2^32 */
--
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