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