diff --git a/include/primecount-internal.hpp b/include/primecount-internal.hpp index d21f064f..b9bd8828 100644 --- a/include/primecount-internal.hpp +++ b/include/primecount-internal.hpp @@ -53,6 +53,7 @@ int64_t Li(int64_t); int64_t Li_inverse(int64_t); int64_t Ri(int64_t); int64_t Ri_inverse(int64_t); +int64_t nth_prime_approx(int64_t n); #ifdef HAVE_INT128_T int128_t pi(int128_t x); diff --git a/src/RiemannR.cpp b/src/RiemannR.cpp index 17edbfa9..b91b8a33 100644 --- a/src/RiemannR.cpp +++ b/src/RiemannR.cpp @@ -427,6 +427,21 @@ int64_t Ri_inverse(int64_t x) return (int64_t) res; } +/// nth_prime_approx(n) is a very accurate approximation of the nth +/// prime with |nth prime - nth_prime_approx(n)| < sqrt(nth prime). +/// Please note that nth_prime_approx(n) may be smaller or larger than +/// the actual nth prime. +/// +int64_t nth_prime_approx(int64_t n) +{ + // Li_inverse(n) is faster but less accurate than Ri_inverse(n). + // For small n speed is more important than accuracy. + if (n < 1e8) + return Li_inverse(n); + else + return Ri_inverse(n); +} + #ifdef HAVE_INT128_T int128_t Li(int128_t x) diff --git a/src/nth_prime.cpp b/src/nth_prime.cpp index 213c0332..6ad182c1 100644 --- a/src/nth_prime.cpp +++ b/src/nth_prime.cpp @@ -98,17 +98,10 @@ int64_t nth_prime(int64_t n, int threads) if (n <= PiTable::pi_cache(PiTable::max_cached())) return binary_search_nth_prime(n); - int64_t prime_approx; - - // Li_inverse(x) is faster but less accurate than Ri_inverse(x). - // For small n speed is more important than accuracy. - if (n < 1e8) - prime_approx = Li_inverse(n); - else - prime_approx = Ri_inverse(n); - - // For large n we use the prime counting function - // and the segmented sieve of Eratosthenes. + // Closely approximate the nth prime using the inverse + // Riemann R function and then count the primes up to this + // approximation using the prime counting function. + int64_t prime_approx = nth_prime_approx(n); int64_t count_approx = pi(prime_approx, threads); int64_t avg_prime_gap = ilog(prime_approx) + 2; int64_t prime = -1;