This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.23
in repository libmath-prime-util-perl.
commit 0965fab050222388d688b7556c6d9c09c9c37dff
Author: Dana Jacobsen <d...@acm.org>
Date: Thu Feb 28 02:11:04 2013 -0800
Add binary search for nth_prime, for inputs > 2e11
---
Changes | 4 ++++
util.c | 44 ++++++++++++++++++++++++++++++++++++++++----
2 files changed, 44 insertions(+), 4 deletions(-)
diff --git a/Changes b/Changes
index 60d17ca..945d23a 100644
--- a/Changes
+++ b/Changes
@@ -15,6 +15,10 @@ Revision history for Perl extension Math::Prime::Util.
- Add examples for first and second Chebyshev functions.
+ - Implement binary search on RiemannR for XS nth_prime when n > 2e11.
+ Runs ~2x faster for 1e12, 3x faster for 1e13. Thanks to Programming
+ Praxis for the idea and motivation.
+
0.22 26 February 2013
- Move main factor loop out of xs and into factor.c.
diff --git a/util.c b/util.c
index fe2fc19..e6f2bfa 100644
--- a/util.c
+++ b/util.c
@@ -557,10 +557,46 @@ UV _XS_nth_prime(UV 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));
- segment_size = lower_limit / 30;
- lower_limit = 30 * segment_size - 1;
- count = _XS_lehmer_pi(lower_limit) - 3;
- MPUassert(count <= target, "Pi(nth_prime_lower(n))) > n");
+#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_lehmer_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_lehmer_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(sqrt(upper_limit));
--
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
