Improved nth Prime (#142)

The new algorithm uses the inverse Riemann R function to accurately approximate the nth prime. It then counts the primes up to this approximate nth prime using multi-threading. Afterwards our new code uses a primesieve::iterator (single-threaded) for the remaining distance < sqrt(n).

CMakeLists.txt

@@ -70,6 +70,7 @@ set(LIB_SRC src/api-c.cpp
             src/MemoryPool.cpp
             src/PrimeGenerator.cpp
             src/nthPrime.cpp
+            src/nthPrimeApprox.cpp
             src/ParallelSieve.cpp
             src/popcount.cpp
             src/PreSieve.cpp

ChangeLog

@@ -1,6 +1,8 @@
-Changes in version 11.2, 02/01/2024
+Changes in version 11.2, 08/01/2024
 ===================================
 
+* nthPrime.cpp: Rewritten using faster nth prime approximation.
+* nthPrimeApprox.cpp: Logarithmic integral and Riemann R implementations.
 * cmake/libatomic.cmake: Fix failed to find libatomic #141.
 * .github/workflows/ci.yml: Port AppVeyor CI tests to GitHub Actions.
 * doc/C_API.md: Fix off by 1 error in OpenMP example #137.

include/primesieve/PrimeSieve.hpp

@@ -78,6 +78,7 @@ class PrimeSieve
   // nth prime
   uint64_t nthPrime(uint64_t);
   uint64_t nthPrime(int64_t, uint64_t);
+  uint64_t negativeNthPrime(int64_t, uint64_t);
   // Count
   counts_t& getCounts();
   uint64_t getCount(int) const;

include/primesieve/nthPrimeApprox.hpp

@@ -0,0 +1,23 @@
+///
+/// @file  nthPrimeApprox.hpp
+///
+/// Copyright (C) 2024 Kim Walisch, <[email protected]>
+///
+/// This file is distributed under the BSD License. See the COPYING
+/// file in the top level directory.
+///
+
+#include <primesieve/Vector.hpp>
+#include <stdint.h>
+
+namespace primesieve {
+
+Vector<int32_t> generate_moebius(int64_t max);
+uint64_t Li(uint64_t x);
+uint64_t Li_inverse(uint64_t x);
+uint64_t Ri(uint64_t x);
+uint64_t Ri_inverse(uint64_t x);
+uint64_t primesApprox(uint64_t x);
+uint64_t nthPrimeApprox(uint64_t n);
+
+} // namespace

src/nthPrime.cpp

@@ -1,7 +1,7 @@
 ///
 /// @file  nthPrime.cpp
 ///
-/// Copyright (C) 2022 Kim Walisch, <[email protected]>
+/// Copyright (C) 2024 Kim Walisch, <[email protected]>
 ///
 /// This file is distributed under the BSD License. See the COPYING
 /// file in the top level directory.
@@ -13,84 +13,31 @@
 #include <primesieve/macros.hpp>
 #include <primesieve/pmath.hpp>
 #include <primesieve/primesieve_error.hpp>
+#include <primesieve/nthPrimeApprox.hpp>
 
 #include <stdint.h>
 #include <algorithm>
 #include <chrono>
 #include <cmath>
-
-using namespace primesieve;
+#include <string>
 
 namespace {
 
-void checkLimit(uint64_t start)
-{
-  if (start >= get_max_stop())
-    throw primesieve_error("nth prime > 2^64");
-}
-
-void checkLowerLimit(uint64_t stop)
-{
-  if (stop == 0)
-    throw primesieve_error("nth prime < 2 is impossible");
-}
-
-bool sieveBackwards(int64_t n, int64_t count, uint64_t stop)
-{
-  return (count >= n) &&
-        !(count == n && stop < 2);
-}
+/// PrimePi(2^64)
+const uint64_t max_n = 425656284035217743ull;
 
-// Prime count approximation
-int64_t pix(int64_t n)
+/// Average prime gap near n
+inline uint64_t avgPrimeGap(uint64_t n)
 {
   double x = (double) n;
-  x = std::max(4.0, x);
-  double pix = x / std::log(x);
-  return (int64_t) pix;
-}
-
-uint64_t nthPrimeDist(int64_t n, int64_t count, uint64_t start)
-{
-  double x = (double) (n - count);
-
-  x = std::abs(x);
-  x = std::max(x, 4.0);
-
-  // rough pi(x) approximation
+  x = std::max(8.0, x);
   double logx = std::log(x);
-  double loglogx = std::log(logx);
-  double pix = x * (logx + loglogx - 1);
-
-  // correct start if sieving backwards to
-  // get more accurate approximation
-  if (count >= n)
-  {
-    double st = start - pix;
-    st = std::max(0.0, st);
-    start = (uint64_t) st;
-  }
-
-  // approximate the nth prime using:
-  // start + n * log(start + pi(n) / loglog(n))
-  double startPix = start + pix / loglogx;
-  startPix = std::max(4.0, startPix);
-  double logStartPix = std::log(startPix);
-  double dist = std::max(pix, x * logStartPix);
-
-  // ensure start + dist <= nth prime
-  if (count < n)
-    dist -= std::sqrt(dist) * std::log(logStartPix) * 2;
-  // ensure start + dist >= nth prime
-  if (count > n)
-    dist += std::sqrt(dist) * std::log(logStartPix) * 2;
-
-  // if n is very small:
-  // ensure start + dist >= nth prime
-  double primeGap = maxPrimeGap(startPix);
-  dist = std::max(dist, primeGap);
-
-  return (uint64_t) dist;
+  // 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
@@ -104,60 +51,124 @@ uint64_t PrimeSieve::nthPrime(uint64_t n)
 
 uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
 {
+  if (n < 0)
+    return negativeNthPrime(n, start);
+  else if (n == 0)
+    n = 1; // like Mathematica
+  else if ((uint64_t) n > max_n)
+    throw primesieve_error("nth_prime(n): n must be <= " + std::to_string(max_n));
+
   setStart(start);
   auto t1 = std::chrono::system_clock::now();
+  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 (n == 0)
-    n = 1; // like Mathematica
-  else if (n > 0)
+  // 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)
+  {
+    // Count primes > start
     start = checkedAdd(start, 1);
-  else if (n < 0)
-    start = checkedSub(start, 1);
-
-  uint64_t stop = start;
-  uint64_t dist = nthPrimeDist(n, 0, start);
-  uint64_t nthPrimeGuess = checkedAdd(start, dist);
-
-  int64_t count = 0;
-  int64_t tinyN = 100000;
-  tinyN = std::max(tinyN, pix(isqrt(nthPrimeGuess)));
+    primeApprox = std::max(start, primeApprox);
+    countApprox = countPrimes(start, primeApprox);
+    start = primeApprox;
+  }
 
-  while ((n - count) > tinyN ||
-         sieveBackwards(n, count, stop))
+  // 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)
   {
-    if (count < n)
-    {
-      checkLimit(start);
-      dist = nthPrimeDist(n, count, start);
-      stop = checkedAdd(start, dist);
-      count += countPrimes(start, stop);
-      start = checkedAdd(stop, 1);
-    }
-    if (sieveBackwards(n, count, stop))
+    start = checkedAdd(start, 1);
+    uint64_t dist = (n - countApprox) * avgPrimeGap(primeApprox);
+    uint64_t stop = checkedAdd(start, dist);
+    primesieve::iterator iter(start, stop);
+    for (int64_t i = countApprox; i < n; i++)
+      prime = iter.next_prime();
+  }
+  else // if (countApprox >= n)
+  {
+    uint64_t dist = (countApprox - n) * avgPrimeGap(primeApprox);
+    uint64_t stop = checkedSub(start, dist);
+    primesieve::iterator iter(start, stop);
+    for (int64_t i = countApprox; i >= n; i--)
     {
-      checkLowerLimit(stop);
-      dist = nthPrimeDist(n, count, stop);
-      start = checkedSub(start, dist);
-      count -= (int64_t) countPrimes(start, stop);
-      stop = checkedSub(start, 1);
+      prime = iter.prev_prime();
+      if_unlikely(prime == 0)
+        throw primesieve_error("nth_prime(n): invalid n, nth prime < 2 is impossible!");
     }
   }
 
-  if (n < 0)
-    count -= 1;
+  auto t2 = std::chrono::system_clock::now();
+  std::chrono::duration<double> seconds = t2 - t1;
+  seconds_ = seconds.count();
 
-  // here start < nth prime,
-  // hence we can sieve forward the remaining
-  // distance and find the nth prime
-  ASSERT(count < n);
+  return prime;
+}
+
+/// Used for n < 0
+uint64_t PrimeSieve::negativeNthPrime(int64_t n, uint64_t start)
+{
+  ASSERT(n < 0);
+  n = -n;
 
-  checkLimit(start);
-  dist = nthPrimeDist(n, count, start) * 2;
-  stop = checkedAdd(start, dist);
+  if ((uint64_t) n >= start)
+    throw primesieve_error("nth_prime(n): abs(n) must be < start");
+  else if ((uint64_t) n > max_n)
+    throw primesieve_error("nth_prime(n): abs(n) must be <= " + std::to_string(max_n));
+
+  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);
+  primeApprox = std::min(primeApprox, start);
   uint64_t prime = 0;
+  int64_t countApprox = 0;
+
+  // 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)
+  {
+    // Count primes < start
+    start = checkedSub(start, 1);
+    primeApprox = std::min(primeApprox, start);
+    countApprox = countPrimes(primeApprox, start);
+    start = primeApprox;
+  }
 
-  for (primesieve::iterator it(start, stop); count < n; count++)
-    prime = it.next_prime();
+  // 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 dist = (countApprox - n) * 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();
+  }
+  else // if (countApprox < n)
+  {
+    start = checkedSub(start, 1);
+    uint64_t dist = (n - countApprox) * 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): invalid n, nth prime < 2 is impossible!");
+    }
+  }
 
   auto t2 = std::chrono::system_clock::now();
   std::chrono::duration<double> seconds = t2 - t1;