diff --git a/include/primesieve/pmath.hpp b/include/primesieve/pmath.hpp index 95ee3c0fb..16924fce6 100644 --- a/include/primesieve/pmath.hpp +++ b/include/primesieve/pmath.hpp @@ -72,12 +72,6 @@ inline T floorPow2(T x) #endif } -template -inline int ilog(T x) -{ - return (int) std::log((double) x); -} - template inline T ilog2(T x) { diff --git a/src/nthPrime.cpp b/src/nthPrime.cpp index cea58e493..f43fff09b 100644 --- a/src/nthPrime.cpp +++ b/src/nthPrime.cpp @@ -26,6 +26,20 @@ namespace { /// PrimePi(2^64) const uint64_t max_n = 425656284035217743ull; +/// Average prime gap near n +inline uint64_t avgPrimeGap(uint64_t n) +{ + double x = (double) n; + x = std::max(8.0, x); + double logx = std::log(x); + // When we buffer primes using primesieve::iterator, we + // want to make sure we buffer primes up to the nth + // prime. Therefore we use +2 here, better to buffer + // slightly too many primes than not enough primes. + double primeGap = logx + 2; + return (uint64_t) primeGap; +} + } // namespace namespace primesieve { @@ -46,63 +60,49 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start) setStart(start); auto t1 = std::chrono::system_clock::now(); - uint64_t primeApprox = start; - uint64_t avgPrimeGap = 0; + uint64_t nApprox = checkedAdd(primesApprox(start), n); + nApprox = std::min(nApprox, max_n); + uint64_t primeApprox = nthPrimeApprox(nApprox); + primeApprox = std::max(primeApprox, start); int64_t countApprox = 0; uint64_t prime = 0; - if (start == 0) - primeApprox = nthPrimeApprox(n); - else - { - ASSERT(start > 0); - uint64_t nApprox = checkedAdd(primesApprox(start), n); - nApprox = std::min(nApprox, max_n); - primeApprox = nthPrimeApprox(nApprox); - } - - primeApprox = std::max(primeApprox, start); - if (primeApprox > 0) - avgPrimeGap = ilog(primeApprox) + 2; - // Only use multi-threading if the sieving distance is sufficiently // large. For small n this if statement also avoids calling // countPrimes() and hence the initialization overhead of // O(x^0.5 log log x^0.5) occurs only once (instead of twice) when // using primesieve::iterator further down. - if (primeApprox - start < isqrt(primeApprox) / 10) - primeApprox = start; - else + if (primeApprox - start > isqrt(primeApprox) / 10) { // Count primes > start start = checkedAdd(start, 1); primeApprox = std::max(start, primeApprox); countApprox = countPrimes(start, primeApprox); + start = primeApprox; } // Here we are very close to the nth prime < sqrt(nth_prime), // we simply iterate over the primes until we find it. if (countApprox < n) { - uint64_t start = primeApprox + 1; - uint64_t stop = checkedAdd(start, (n - countApprox) * avgPrimeGap); + start = checkedAdd(start, 1); + uint64_t dist = (n - countApprox) * avgPrimeGap(start); + uint64_t stop = checkedAdd(start, dist); primesieve::iterator iter(start, stop); - int64_t i = countApprox; - while (i < n) - { + for (int64_t i = countApprox; i < n; i++) prime = iter.next_prime(); - i += 1; - if_unlikely(i < n && prime == 18446744073709551557ull) - throw primesieve_error("nth_prime(n) > 2^64 is not supported!"); - } } else // if (countApprox >= n) { - uint64_t start = primeApprox; - uint64_t stop = checkedSub(start, (countApprox - n) * avgPrimeGap); + uint64_t dist = (countApprox - n) * avgPrimeGap(start); + uint64_t stop = checkedSub(start, dist); primesieve::iterator iter(start, stop); for (int64_t i = countApprox; i >= n; i--) + { prime = iter.prev_prime(); + if_unlikely(prime == 0) + throw primesieve_error("nth_prime(n): n is too small, nth_prime(n) < 2!"); + } } auto t2 = std::chrono::system_clock::now(); @@ -125,52 +125,49 @@ uint64_t PrimeSieve::negativeNthPrime(int64_t n, uint64_t start) setStart(start); auto t1 = std::chrono::system_clock::now(); uint64_t nApprox = checkedSub(primesApprox(start), n); + nApprox = std::min(nApprox, max_n); uint64_t primeApprox = nthPrimeApprox(nApprox); - uint64_t avgPrimeGap = 0; + primeApprox = std::min(primeApprox, start); uint64_t prime = 0; int64_t countApprox = 0; - primeApprox = std::min(primeApprox, start); - if (primeApprox > 0) - avgPrimeGap = ilog(primeApprox) + 2; - // Only use multi-threading if the sieving distance is sufficiently // large. For small n this if statement also avoids calling // countPrimes() and hence the initialization overhead of // O(x^0.5 log log x^0.5) occurs only once (instead of twice) when // using primesieve::iterator further down. - if (start - primeApprox < isqrt(start) / 10) - primeApprox = start; - else + if (start - primeApprox > isqrt(start) / 10) { // Count primes < start start = checkedSub(start, 1); primeApprox = std::min(primeApprox, start); countApprox = countPrimes(primeApprox, start); + start = primeApprox; } // Here we are very close to the nth prime < sqrt(nth_prime), // we simply iterate over the primes until we find it. - if (countApprox >= n) - { - uint64_t start = primeApprox; - uint64_t stop = checkedAdd(start, (n - countApprox) * avgPrimeGap); - primesieve::iterator iter(start, stop); - for (int64_t i = countApprox; i >= n; i--) - prime = iter.next_prime(); - } - else // if (countApprox < n) + if (countApprox < n) { - uint64_t start = checkedSub(primeApprox, 1); - uint64_t stop = checkedSub(start, (countApprox - n) * avgPrimeGap); + start = checkedSub(start, 1); + uint64_t dist = (countApprox - n) * avgPrimeGap(start); + uint64_t stop = checkedSub(start, dist); primesieve::iterator iter(start, stop); - for (int64_t i = countApprox; i < n; i += 1) + for (int64_t i = countApprox; i < n; i++) { prime = iter.prev_prime(); if_unlikely(prime == 0) throw primesieve_error("nth_prime(n): n is too small, nth_prime(n) < 2!"); } } + else // if (countApprox >= n) + { + uint64_t dist = (n - countApprox) * avgPrimeGap(start); + uint64_t stop = checkedAdd(start, dist); + primesieve::iterator iter(start, stop); + for (int64_t i = countApprox; i >= n; i--) + prime = iter.next_prime(); + } auto t2 = std::chrono::system_clock::now(); std::chrono::duration seconds = t2 - t1;