From 7b7623bdd931b50f1e31cc31e6434271c1bd8b55 Mon Sep 17 00:00:00 2001
From: Kim Walisch <kim.walisch@gmail.com>
Date: Sun, 7 Jan 2024 10:31:34 +0100
Subject: [PATCH 01/22] Initial commit

---
 CMakeLists.txt                        |   1 +
 include/primesieve/nthPrimeApprox.hpp |  23 ++
 include/primesieve/pmath.hpp          |   6 +
 src/nthPrime.cpp                      | 151 +++----------
 src/nthPrimeApprox.cpp                | 295 ++++++++++++++++++++++++++
 5 files changed, 357 insertions(+), 119 deletions(-)
 create mode 100644 include/primesieve/nthPrimeApprox.hpp
 create mode 100644 src/nthPrimeApprox.cpp

diff --git a/CMakeLists.txt b/CMakeLists.txt
index fe639d821..d92fe472c 100644
--- a/CMakeLists.txt
+++ b/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
diff --git a/include/primesieve/nthPrimeApprox.hpp b/include/primesieve/nthPrimeApprox.hpp
new file mode 100644
index 000000000..58e68891c
--- /dev/null
+++ b/include/primesieve/nthPrimeApprox.hpp
@@ -0,0 +1,23 @@
+///
+/// @file  nthPrimeApprox.hpp
+///
+/// Copyright (C) 2024 Kim Walisch, <kim.walisch@gmail.com>
+///
+/// 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
diff --git a/include/primesieve/pmath.hpp b/include/primesieve/pmath.hpp
index 16924fce6..95ee3c0fb 100644
--- a/include/primesieve/pmath.hpp
+++ b/include/primesieve/pmath.hpp
@@ -72,6 +72,12 @@ inline T floorPow2(T x)
 #endif
 }
 
+template <typename T>
+inline int ilog(T x)
+{
+  return (int) std::log((double) x);
+}
+
 template <typename T>
 inline T ilog2(T x)
 {
diff --git a/src/nthPrime.cpp b/src/nthPrime.cpp
index 631ced9cd..fc04e1e98 100644
--- a/src/nthPrime.cpp
+++ b/src/nthPrime.cpp
@@ -1,7 +1,7 @@
 ///
 /// @file  nthPrime.cpp
 ///
-/// Copyright (C) 2022 Kim Walisch, <kim.walisch@gmail.com>
+/// Copyright (C) 2024 Kim Walisch, <kim.walisch@gmail.com>
 ///
 /// This file is distributed under the BSD License. See the COPYING
 /// file in the top level directory.
@@ -13,88 +13,13 @@
 #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;
-
-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);
-}
-
-// Prime count approximation
-int64_t pix(int64_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
-  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;
-}
-
-} // namespace
-
 namespace primesieve {
 
 uint64_t PrimeSieve::nthPrime(uint64_t n)
@@ -109,55 +34,43 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
 
   if (n == 0)
     n = 1; // like Mathematica
-  else if (n > 0)
+  //if (n > max_n)
+  //  throw primecount_error("nth_prime(n): n must be <= " + std::to_string(max_n));
+
+  uint64_t prime_approx;
+
+  if (start == 0) 
+    prime_approx = nthPrimeApprox(n);
+  else // start > 0
+    prime_approx = nthPrimeApprox(primesApprox(start) + n);
+
+  if (n > 0)
     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)));
+  int64_t count_approx = countPrimes(start, prime_approx);
+  int64_t avg_prime_gap =  ilog(prime_approx) + 2;
+  int64_t prime = -1;
 
-  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 (count_approx < 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))
-    {
-      checkLowerLimit(stop);
-      dist = nthPrimeDist(n, count, stop);
-      start = checkedSub(start, dist);
-      count -= (int64_t) countPrimes(start, stop);
-      stop = checkedSub(start, 1);
-    }
+    uint64_t start = prime_approx + 1;
+    uint64_t stop = start + (n - count_approx) * avg_prime_gap;
+    primesieve::iterator iter(start, stop);
+    for (int64_t i = count_approx; i < n; i++)
+      prime = iter.next_prime();
+  }
+  else // if (count_approx >= n)
+  {
+    uint64_t start = prime_approx;
+    uint64_t stop = start - (count_approx - n) * avg_prime_gap;
+    primesieve::iterator iter(start, stop);
+    for (int64_t i = count_approx; i + 1 > n; i--)
+      prime = iter.prev_prime();
   }
-
-  if (n < 0)
-    count -= 1;
-
-  // here start < nth prime,
-  // hence we can sieve forward the remaining
-  // distance and find the nth prime
-  ASSERT(count < n);
-
-  checkLimit(start);
-  dist = nthPrimeDist(n, count, start) * 2;
-  stop = checkedAdd(start, dist);
-  uint64_t prime = 0;
-
-  for (primesieve::iterator it(start, stop); count < n; count++)
-    prime = it.next_prime();
 
   auto t2 = std::chrono::system_clock::now();
   std::chrono::duration<double> seconds = t2 - t1;
diff --git a/src/nthPrimeApprox.cpp b/src/nthPrimeApprox.cpp
new file mode 100644
index 000000000..38ceb69fd
--- /dev/null
+++ b/src/nthPrimeApprox.cpp
@@ -0,0 +1,295 @@
+///
+/// @file  nthPrimeApprox.cpp
+///
+/// Copyright (C) 2024 Kim Walisch, <kim.walisch@gmail.com>
+///
+/// This file is distributed under the BSD License. See the COPYING
+/// file in the top level directory.
+///
+
+#include <primesieve/pmath.hpp>
+#include <primesieve/Vector.hpp>
+
+#include <stdint.h>
+#include <cmath>
+#include <limits>
+
+namespace primesieve {
+
+/// Generate a vector with Möbius function values.
+/// This implementation is based on code by Rick Sladkey:
+/// https://mathoverflow.net/q/99545
+///
+Vector<int32_t> generate_moebius(int64_t max)
+{
+  int64_t sqrt = isqrt(max);
+  int64_t size = max + 1;
+  Vector<int32_t> mu(size);
+  std::fill(mu.begin(), mu.end(), 1);
+
+  for (int64_t i = 2; i <= sqrt; i++)
+  {
+    if (mu[i] == 1)
+    {
+      for (int64_t j = i; j < size; j += i)
+        mu[j] *= (int32_t) -i;
+      for (int64_t j = i * i; j < size; j += i * i)
+        mu[j] = 0;
+    }
+  }
+
+  for (int64_t i = 2; i < size; i++)
+  {
+    if (mu[i] == i)
+      mu[i] = 1;
+    else if (mu[i] == -i)
+      mu[i] = -1;
+    else if (mu[i] < 0)
+      mu[i] = 1;
+    else if (mu[i] > 0)
+      mu[i] = -1;
+  }
+
+  return mu;
+}
+
+} // namespace
+
+namespace {
+
+/// Calculate the logarithmic integral using
+/// Ramanujan's formula:
+/// https://en.wikipedia.org/wiki/Logarithmic_integral_function#Series_representation
+///
+long double li(long double x)
+{
+  if (x <= 1)
+    return 0;
+
+  long double gamma = 0.577215664901532860606512090082402431L;
+  long double sum = 0;
+  long double inner_sum = 0;
+  long double factorial = 1;
+  long double p = -1;
+  long double q = 0;
+  long double power2 = 1;
+  long double logx = std::log(x);
+  int k = 0;
+
+  for (int n = 1; true; n++)
+  {
+    p *= -logx;
+    factorial *= n;
+    q = factorial * power2;
+    power2 *= 2;
+
+    for (; k <= (n - 1) / 2; k++)
+      inner_sum += 1.0L / (2 * k + 1);
+
+    auto old_sum = sum;
+    sum += (p / q) * inner_sum;
+
+    // Not converging anymore
+    if (std::abs(sum - old_sum) < std::numeric_limits<long double>::epsilon())
+      break;
+  }
+
+  return gamma + std::log(logx) + std::sqrt(x) * sum;
+}
+
+/// Calculate the offset logarithmic integral which is a very
+/// accurate approximation of the number of primes <= x.
+/// Li(x) > pi(x) for 24 <= x <= ~ 10^316
+///
+long double Li(long double x)
+{
+  long double li2 = 1.045163780117492784844588889194613136L;
+
+  if (x <= li2)
+    return 0;
+  else
+    return li(x) - li2;
+}
+
+/// Calculate the inverse offset logarithmic integral which
+/// is a very accurate approximation of the nth prime.
+/// Li^-1(x) < nth_prime(x) for 7 <= x <= 10^316
+///
+/// This implementation computes Li^-1(x) as the zero of the
+/// function f(z) = Li(z) - x using the Newton–Raphson method.
+/// Note that Li'(z) = 1 / log(z).
+/// https://math.stackexchange.com/a/853192
+///
+/// Newton–Raphson method:
+/// zn+1 = zn - (f(zn) / f'(zn)).
+/// zn+1 = zn - (Li(zn) - x) / (1 / log(zn))
+/// zn+1 = zn - (Li(zn) - x) * log(zn)
+///
+long double Li_inverse(long double x)
+{
+  if (x < 2)
+    return 0;
+
+  long double t = x * std::log(x);
+  long double old_term = std::numeric_limits<long double>::infinity();
+
+  while (true)
+  {
+    long double term = (Li(t) - x) * std::log(t);
+
+    // Not converging anymore
+    if (std::abs(term) >= std::abs(old_term))
+      break;
+
+    t -= term;
+    old_term = term;
+  }
+
+  return t;
+}
+
+/// Calculate the Riemann R function which is a very accurate
+/// approximation of the number of primes below x.
+/// RiemannR(x) = \sum_{n=1}^{∞} μ(n)/n * li(x^(1/n))
+/// http://mathworld.wolfram.com/RiemannPrimeCountingFunction.html
+///
+long double Ri(long double x)
+{
+  if (x <= 1)
+    return 0;
+
+  long double sum = 0;
+  long double old_term = std::numeric_limits<long double>::infinity();
+  auto terms = (int) (std::log2(x) * 2 + 10);
+  auto mu = primesieve::generate_moebius(terms);
+
+  for (int n = 1; n < terms; n++)
+  {
+    if (mu[n])
+    {
+      long double root = std::pow(x, 1.0L / n);
+      long double term = (li(root) * mu[n]) / n;
+
+      // Not converging anymore
+      if (std::abs(term) >= std::abs(old_term))
+        break;
+
+      sum += term;
+      old_term = term;
+    }
+  }
+
+  return sum;
+}
+
+/// Calculate the inverse Riemann R function which is a very
+/// accurate approximation of the nth prime.
+/// This implementation computes Ri^-1(x) as the zero of the
+/// function f(z) = Ri(z) - x using the Newton–Raphson method.
+/// Note that Ri'(z) = 1 / log(z).
+/// https://math.stackexchange.com/a/853192
+///
+/// Newton–Raphson method:
+/// zn+1 = zn - (f(zn) / f'(zn)).
+/// zn+1 = zn - (Ri(zn) - x) / (1 / log(zn))
+/// zn+1 = zn - (Ri(zn) - x) * log(zn)
+///
+long double Ri_inverse(long double x)
+{
+  if (x < 2)
+    return 0;
+
+  long double t = Li_inverse(x);
+  long double old_term = std::numeric_limits<long double>::infinity();
+
+  while (true)
+  {
+    long double term = (Ri(t) - x) * std::log(t);
+
+    // Not converging anymore
+    if (std::abs(term) >= std::abs(old_term))
+      break;
+
+    t -= term;
+    old_term = term;
+  }
+
+  return t;
+}
+
+} // namespace
+
+namespace {
+
+using primesieve::Array;
+
+// tinyPrimes[1] = 2, tinyPrimes[2] = 3, ...
+const Array<int16_t, 170> tinyPrimes =
+{
+    0,   2,   3,   5,   7,  11,  13,  17,  19,  23, 
+   29,  31,  37,  41,  43,  47,  53,  59,  61,  67, 
+   71,  73,  79,  83,  89,  97, 101, 103, 107, 109, 
+  113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 
+  173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 
+  229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 
+  281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 
+  349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 
+  409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 
+  463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 
+  541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 
+  601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 
+  659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 
+  733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 
+  809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 
+  863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 
+  941, 947, 953, 967, 971, 977, 983, 991, 997, 1009
+};
+
+} //namespace
+
+namespace primesieve {
+
+uint64_t Li(uint64_t x)
+{
+  return (uint64_t) ::Li((long double) x);
+}
+
+uint64_t Li_inverse(uint64_t x)
+{
+  return (uint64_t) ::Li_inverse((long double) x);
+}
+
+uint64_t Ri(uint64_t x)
+{
+  return (uint64_t) ::Ri((long double) x);
+}
+
+uint64_t Ri_inverse(uint64_t x)
+{
+  return (uint64_t) ::Ri_inverse((long double) x);
+}
+
+uint64_t primesApprox(uint64_t x)
+{
+  // Li(x) is faster but less accurate than Ri(x).
+  // For small n speed is more important than accuracy.
+  if (x < 1e10)
+    return Li(x);
+  else
+    return Ri(x);
+}
+
+uint64_t nthPrimeApprox(uint64_t n)
+{
+  if (n < tinyPrimes.size())
+    return tinyPrimes[n];
+
+  // Li_inverse(x) is faster but less accurate than Ri_inverse(x).
+  // For small n speed is more important than accuracy.
+  if (n < 1e8)
+    return Li_inverse(n);
+  else
+    return Ri_inverse(n);
+}
+
+} // namespace

From 50fb93b6c01478c3aa54a6bfe22ff3b30dea2e91 Mon Sep 17 00:00:00 2001
From: Kim Walisch <kim.walisch@gmail.com>
Date: Sun, 7 Jan 2024 10:46:21 +0100
Subject: [PATCH 02/22] Add max_n check

---
 src/nthPrime.cpp | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/src/nthPrime.cpp b/src/nthPrime.cpp
index fc04e1e98..d77f0f6ee 100644
--- a/src/nthPrime.cpp
+++ b/src/nthPrime.cpp
@@ -20,6 +20,13 @@
 #include <chrono>
 #include <cmath>
 
+namespace {
+
+/// PrimePi(2^64)
+const int64_t max_n = 425656284035217743ull;
+
+} // namespace
+
 namespace primesieve {
 
 uint64_t PrimeSieve::nthPrime(uint64_t n)
@@ -34,8 +41,8 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
 
   if (n == 0)
     n = 1; // like Mathematica
-  //if (n > max_n)
-  //  throw primecount_error("nth_prime(n): n must be <= " + std::to_string(max_n));
+  if (n > max_n)
+    throw primesieve_error("nth_prime(n): n must be <= " + std::to_string(max_n));
 
   uint64_t prime_approx;
 

From 506628d67044b7b71e9f0a098e192335c7e4d148 Mon Sep 17 00:00:00 2001
From: Kim Walisch <kim.walisch@gmail.com>
Date: Sun, 7 Jan 2024 10:56:53 +0100
Subject: [PATCH 03/22] Fix int types

---
 src/nthPrime.cpp | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/nthPrime.cpp b/src/nthPrime.cpp
index d77f0f6ee..81bf5a1f6 100644
--- a/src/nthPrime.cpp
+++ b/src/nthPrime.cpp
@@ -57,15 +57,15 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
     start = checkedSub(start, 1);
 
   int64_t count_approx = countPrimes(start, prime_approx);
-  int64_t avg_prime_gap =  ilog(prime_approx) + 2;
-  int64_t prime = -1;
+  uint64_t avg_prime_gap =  ilog(prime_approx) + 2;
+  uint64_t prime = -1;
 
   // Here we are very close to the nth prime < sqrt(nth_prime),
   // we simply iterate over the primes until we find it.
   if (count_approx < n)
   {
     uint64_t start = prime_approx + 1;
-    uint64_t stop = start + (n - count_approx) * avg_prime_gap;
+    uint64_t stop = checkedAdd(start, (n - count_approx) * avg_prime_gap);
     primesieve::iterator iter(start, stop);
     for (int64_t i = count_approx; i < n; i++)
       prime = iter.next_prime();
@@ -73,7 +73,7 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
   else // if (count_approx >= n)
   {
     uint64_t start = prime_approx;
-    uint64_t stop = start - (count_approx - n) * avg_prime_gap;
+    uint64_t stop = checkedSub(start, (count_approx - n) * avg_prime_gap);
     primesieve::iterator iter(start, stop);
     for (int64_t i = count_approx; i + 1 > n; i--)
       prime = iter.prev_prime();

From 406fdefba91f8dc7e1e3c0e01b712dcc9de88b3f Mon Sep 17 00:00:00 2001
From: Kim Walisch <kim.walisch@gmail.com>
Date: Sun, 7 Jan 2024 11:29:34 +0100
Subject: [PATCH 04/22] Improve error handling

---
 src/nthPrime.cpp | 17 ++++++++++++++---
 1 file changed, 14 insertions(+), 3 deletions(-)

diff --git a/src/nthPrime.cpp b/src/nthPrime.cpp
index 81bf5a1f6..6b7d6f8cf 100644
--- a/src/nthPrime.cpp
+++ b/src/nthPrime.cpp
@@ -48,8 +48,13 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
 
   if (start == 0) 
     prime_approx = nthPrimeApprox(n);
-  else // start > 0
-    prime_approx = nthPrimeApprox(primesApprox(start) + n);
+  else
+  {
+    ASSERT(start > 0);
+    uint64_t new_n = checkedAdd(primesApprox(start), n);
+    new_n = std::min(new_n, max_n);
+    prime_approx = nthPrimeApprox(new_n);
+  }
 
   if (n > 0)
     start = checkedAdd(start, 1);
@@ -67,8 +72,14 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
     uint64_t start = prime_approx + 1;
     uint64_t stop = checkedAdd(start, (n - count_approx) * avg_prime_gap);
     primesieve::iterator iter(start, stop);
-    for (int64_t i = count_approx; i < n; i++)
+    int64_t i = count_approx;
+    while (i < n)
+    {
       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 (count_approx >= n)
   {

From a842d93e6b7d64db89912b9edd5261582fcccfb7 Mon Sep 17 00:00:00 2001
From: Kim Walisch <kim.walisch@gmail.com>
Date: Sun, 7 Jan 2024 11:43:03 +0100
Subject: [PATCH 05/22] Refactor

---
 src/nthPrime.cpp | 19 ++++++++-----------
 1 file changed, 8 insertions(+), 11 deletions(-)

diff --git a/src/nthPrime.cpp b/src/nthPrime.cpp
index 6b7d6f8cf..9a18447d5 100644
--- a/src/nthPrime.cpp
+++ b/src/nthPrime.cpp
@@ -36,14 +36,15 @@ uint64_t PrimeSieve::nthPrime(uint64_t n)
 
 uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
 {
-  setStart(start);
-  auto t1 = std::chrono::system_clock::now();
-
-  if (n == 0)
+  if (n < 0)
+    return negativeNthPrime(n, start);
+  else if (n == 0)
     n = 1; // like Mathematica
-  if (n > max_n)
+  else if (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 prime_approx;
 
   if (start == 0) 
@@ -56,13 +57,9 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
     prime_approx = nthPrimeApprox(new_n);
   }
 
-  if (n > 0)
-    start = checkedAdd(start, 1);
-  else if (n < 0)
-    start = checkedSub(start, 1);
-
+  start = checkedAdd(start, 1);
   int64_t count_approx = countPrimes(start, prime_approx);
-  uint64_t avg_prime_gap =  ilog(prime_approx) + 2;
+  uint64_t avg_prime_gap = ilog(prime_approx) + 2;
   uint64_t prime = -1;
 
   // Here we are very close to the nth prime < sqrt(nth_prime),

From cfd196f29f4e1e45ca55ddadd7f345da69fc2b03 Mon Sep 17 00:00:00 2001
From: Kim Walisch <kim.walisch@gmail.com>
Date: Sun, 7 Jan 2024 11:45:53 +0100
Subject: [PATCH 06/22] Include <string>

---
 src/nthPrime.cpp | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/nthPrime.cpp b/src/nthPrime.cpp
index 9a18447d5..79dc7cce2 100644
--- a/src/nthPrime.cpp
+++ b/src/nthPrime.cpp
@@ -19,6 +19,7 @@
 #include <algorithm>
 #include <chrono>
 #include <cmath>
+#include <string>
 
 namespace {
 

From 6618afcea0953420b192826b562b7ee787d71642 Mon Sep 17 00:00:00 2001
From: Kim Walisch <kim.walisch@gmail.com>
Date: Sun, 7 Jan 2024 17:53:05 +0100
Subject: [PATCH 07/22] More changes

---
 include/primesieve/PrimeSieve.hpp |  1 +
 src/nthPrime.cpp                  | 90 +++++++++++++++++++++++++------
 2 files changed, 74 insertions(+), 17 deletions(-)

diff --git a/include/primesieve/PrimeSieve.hpp b/include/primesieve/PrimeSieve.hpp
index 59b58d911..ba5fb529c 100644
--- a/include/primesieve/PrimeSieve.hpp
+++ b/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;
diff --git a/src/nthPrime.cpp b/src/nthPrime.cpp
index 79dc7cce2..fd3ed6699 100644
--- a/src/nthPrime.cpp
+++ b/src/nthPrime.cpp
@@ -24,7 +24,7 @@
 namespace {
 
 /// PrimePi(2^64)
-const int64_t max_n = 425656284035217743ull;
+const uint64_t max_n = 425656284035217743ull;
 
 } // namespace
 
@@ -41,36 +41,37 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
     return negativeNthPrime(n, start);
   else if (n == 0)
     n = 1; // like Mathematica
-  else if (n > max_n)
+  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 prime_approx;
+  uint64_t primeApprox;
 
   if (start == 0) 
-    prime_approx = nthPrimeApprox(n);
+    primeApprox = nthPrimeApprox(n);
   else
   {
     ASSERT(start > 0);
-    uint64_t new_n = checkedAdd(primesApprox(start), n);
-    new_n = std::min(new_n, max_n);
-    prime_approx = nthPrimeApprox(new_n);
+    uint64_t nApprox = checkedAdd(primesApprox(start), n);
+    nApprox = std::min(nApprox, max_n);
+    primeApprox = nthPrimeApprox(nApprox);
   }
 
+  // Count primes > start
   start = checkedAdd(start, 1);
-  int64_t count_approx = countPrimes(start, prime_approx);
-  uint64_t avg_prime_gap = ilog(prime_approx) + 2;
+  int64_t countApprox = countPrimes(start, primeApprox);
+  uint64_t avgPrimeGap = ilog(primeApprox) + 2;
   uint64_t prime = -1;
 
   // Here we are very close to the nth prime < sqrt(nth_prime),
   // we simply iterate over the primes until we find it.
-  if (count_approx < n)
+  if (countApprox < n)
   {
-    uint64_t start = prime_approx + 1;
-    uint64_t stop = checkedAdd(start, (n - count_approx) * avg_prime_gap);
+    uint64_t start = primeApprox + 1;
+    uint64_t stop = checkedAdd(start, (n - countApprox) * avgPrimeGap);
     primesieve::iterator iter(start, stop);
-    int64_t i = count_approx;
+    int64_t i = countApprox;
     while (i < n)
     {
       prime = iter.next_prime();
@@ -79,12 +80,67 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
         throw primesieve_error("nth_prime(n) > 2^64 is not supported!");
     }
   }
-  else // if (count_approx >= n)
+  else // if (countApprox >= n)
   {
-    uint64_t start = prime_approx;
-    uint64_t stop = checkedSub(start, (count_approx - n) * avg_prime_gap);
+    uint64_t start = primeApprox;
+    uint64_t stop = checkedSub(start, (countApprox - n) * avgPrimeGap);
     primesieve::iterator iter(start, stop);
-    for (int64_t i = count_approx; i + 1 > n; i--)
+    for (int64_t i = countApprox; i + 1 > n; i--)
+      prime = iter.prev_prime();
+  }
+
+  auto t2 = std::chrono::system_clock::now();
+  std::chrono::duration<double> seconds = t2 - t1;
+  seconds_ = seconds.count();
+
+  return prime;
+}
+
+uint64_t PrimeSieve::negativeNthPrime(int64_t n, uint64_t start)
+{
+  ASSERT(n < 0);
+  n = -n;
+
+  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): n must be <= " + std::to_string(max_n));
+
+  setStart(start);
+  auto t1 = std::chrono::system_clock::now();
+  uint64_t nApprox = checkedSub(primesApprox(start), n);
+  uint64_t primeApprox = nthPrimeApprox(nApprox);
+
+  // Count primes < start
+  start = checkedSub(start, 1);
+  primeApprox = std::min(primeApprox, start);
+
+  int64_t countApprox = countPrimes(primeApprox, start);
+  uint64_t avgPrimeGap = ilog(primeApprox) + 2;
+  uint64_t prime = -1;
+
+  // 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);
+    int64_t i = countApprox;
+    while (i < n)
+    {
+      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 = checkedSub(primeApprox, 1);
+    uint64_t stop = checkedSub(start, (countApprox - n) * avgPrimeGap);
+    primesieve::iterator iter(start, stop);
+    for (int64_t i = countApprox; i + 1 > n; i--)
       prime = iter.prev_prime();
   }
 

From e20ff0ac58fd38f9a206866fde8f6bfb698e8720 Mon Sep 17 00:00:00 2001
From: Kim Walisch <kim.walisch@gmail.com>
Date: Sun, 7 Jan 2024 19:03:06 +0100
Subject: [PATCH 08/22] Fix issues

---
 src/nthPrime.cpp       | 42 ++++++++++++++++++++++--------------------
 src/nthPrimeApprox.cpp | 31 -------------------------------
 2 files changed, 22 insertions(+), 51 deletions(-)

diff --git a/src/nthPrime.cpp b/src/nthPrime.cpp
index fd3ed6699..952d86f22 100644
--- a/src/nthPrime.cpp
+++ b/src/nthPrime.cpp
@@ -46,7 +46,9 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
 
   setStart(start);
   auto t1 = std::chrono::system_clock::now();
-  uint64_t primeApprox;
+  uint64_t primeApprox = 0;
+  uint64_t avgPrimeGap = 0;
+  uint64_t prime = 0;
 
   if (start == 0) 
     primeApprox = nthPrimeApprox(n);
@@ -60,9 +62,10 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
 
   // Count primes > start
   start = checkedAdd(start, 1);
+  primeApprox = std::max(start, primeApprox);
   int64_t countApprox = countPrimes(start, primeApprox);
-  uint64_t avgPrimeGap = ilog(primeApprox) + 2;
-  uint64_t prime = -1;
+  if (primeApprox > 0)
+    avgPrimeGap = ilog(primeApprox) + 2;
 
   // Here we are very close to the nth prime < sqrt(nth_prime),
   // we simply iterate over the primes until we find it.
@@ -82,10 +85,10 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
   }
   else // if (countApprox >= n)
   {
-    uint64_t start = primeApprox;
-    uint64_t stop = checkedSub(start, (countApprox - n) * avgPrimeGap);
+    uint64_t stop = primeApprox;
+    uint64_t start = checkedSub(stop, (countApprox - n) * avgPrimeGap);
     primesieve::iterator iter(start, stop);
-    for (int64_t i = countApprox; i + 1 > n; i--)
+    for (int64_t i = countApprox; i >= n; i--)
       prime = iter.prev_prime();
   }
 
@@ -104,44 +107,43 @@ uint64_t PrimeSieve::negativeNthPrime(int64_t n, uint64_t start)
   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): n must be <= " + std::to_string(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);
   uint64_t primeApprox = nthPrimeApprox(nApprox);
+  uint64_t avgPrimeGap = 0;
+  uint64_t prime = 0;
 
   // Count primes < start
   start = checkedSub(start, 1);
   primeApprox = std::min(primeApprox, start);
-
   int64_t countApprox = countPrimes(primeApprox, start);
-  uint64_t avgPrimeGap = ilog(primeApprox) + 2;
-  uint64_t prime = -1;
+  if (primeApprox > 0)
+    avgPrimeGap = ilog(primeApprox) + 2;
 
   // 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 (countApprox >= n)
   {
     uint64_t start = primeApprox;
     uint64_t stop = checkedAdd(start, (n - countApprox) * avgPrimeGap);
     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)
+  else // if (countApprox < n)
   {
     uint64_t start = checkedSub(primeApprox, 1);
     uint64_t stop = checkedSub(start, (countApprox - n) * avgPrimeGap);
     primesieve::iterator iter(start, stop);
-    for (int64_t i = countApprox; i + 1 > n; i--)
+    for (int64_t i = countApprox; i < n; i += 1)
+    {
       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();
diff --git a/src/nthPrimeApprox.cpp b/src/nthPrimeApprox.cpp
index 38ceb69fd..c65f6b816 100644
--- a/src/nthPrimeApprox.cpp
+++ b/src/nthPrimeApprox.cpp
@@ -219,34 +219,6 @@ long double Ri_inverse(long double x)
 
 } // namespace
 
-namespace {
-
-using primesieve::Array;
-
-// tinyPrimes[1] = 2, tinyPrimes[2] = 3, ...
-const Array<int16_t, 170> tinyPrimes =
-{
-    0,   2,   3,   5,   7,  11,  13,  17,  19,  23, 
-   29,  31,  37,  41,  43,  47,  53,  59,  61,  67, 
-   71,  73,  79,  83,  89,  97, 101, 103, 107, 109, 
-  113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 
-  173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 
-  229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 
-  281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 
-  349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 
-  409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 
-  463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 
-  541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 
-  601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 
-  659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 
-  733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 
-  809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 
-  863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 
-  941, 947, 953, 967, 971, 977, 983, 991, 997, 1009
-};
-
-} //namespace
-
 namespace primesieve {
 
 uint64_t Li(uint64_t x)
@@ -281,9 +253,6 @@ uint64_t primesApprox(uint64_t x)
 
 uint64_t nthPrimeApprox(uint64_t n)
 {
-  if (n < tinyPrimes.size())
-    return tinyPrimes[n];
-
   // Li_inverse(x) is faster but less accurate than Ri_inverse(x).
   // For small n speed is more important than accuracy.
   if (n < 1e8)

From d90f8440886588b92889fdd651e24bf827e64cd1 Mon Sep 17 00:00:00 2001
From: Kim Walisch <kim.walisch@gmail.com>
Date: Sun, 7 Jan 2024 19:10:20 +0100
Subject: [PATCH 09/22] Fix issue

---
 src/nthPrime.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/nthPrime.cpp b/src/nthPrime.cpp
index 952d86f22..7868b756d 100644
--- a/src/nthPrime.cpp
+++ b/src/nthPrime.cpp
@@ -85,8 +85,8 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
   }
   else // if (countApprox >= n)
   {
-    uint64_t stop = primeApprox;
-    uint64_t start = checkedSub(stop, (countApprox - n) * avgPrimeGap);
+    uint64_t start = primeApprox;
+    uint64_t stop = checkedSub(start, (countApprox - n) * avgPrimeGap);
     primesieve::iterator iter(start, stop);
     for (int64_t i = countApprox; i >= n; i--)
       prime = iter.prev_prime();

From e4600af5d818f95da035ea7685a38169174b8027 Mon Sep 17 00:00:00 2001
From: Kim Walisch <kim.walisch@gmail.com>
Date: Sun, 7 Jan 2024 19:58:58 +0100
Subject: [PATCH 10/22] Faster

---
 src/nthPrime.cpp | 21 ++++++++++++++++-----
 1 file changed, 16 insertions(+), 5 deletions(-)

diff --git a/src/nthPrime.cpp b/src/nthPrime.cpp
index 7868b756d..d4af64744 100644
--- a/src/nthPrime.cpp
+++ b/src/nthPrime.cpp
@@ -46,11 +46,12 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
 
   setStart(start);
   auto t1 = std::chrono::system_clock::now();
+  int64_t countApprox = 0;
   uint64_t primeApprox = 0;
   uint64_t avgPrimeGap = 0;
   uint64_t prime = 0;
 
-  if (start == 0) 
+  if (start == 0)
     primeApprox = nthPrimeApprox(n);
   else
   {
@@ -60,13 +61,23 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
     primeApprox = nthPrimeApprox(nApprox);
   }
 
-  // Count primes > start
-  start = checkedAdd(start, 1);
-  primeApprox = std::max(start, primeApprox);
-  int64_t countApprox = countPrimes(start, primeApprox);
   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 initailization 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 &&
+      primeApprox - start > isqrt(primeApprox - start) / 10)
+  {
+    // Count primes > start
+    start = checkedAdd(start, 1);
+    primeApprox = std::max(start, primeApprox);
+    countApprox = countPrimes(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)

From 71e2303ecc6d09857768cd7a91456ffb4254bd0e Mon Sep 17 00:00:00 2001
From: Kim Walisch <kim.walisch@gmail.com>
Date: Sun, 7 Jan 2024 20:26:05 +0100
Subject: [PATCH 11/22] Fix threshold

---
 src/nthPrime.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/nthPrime.cpp b/src/nthPrime.cpp
index d4af64744..f6185ba42 100644
--- a/src/nthPrime.cpp
+++ b/src/nthPrime.cpp
@@ -70,7 +70,7 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
   // O(x^0.5 log log x^0.5) occurs only once (instead of twice) when
   // using primesieve::iterator further down.
   if (start < primeApprox &&
-      primeApprox - start > isqrt(primeApprox - start) / 10)
+      primeApprox - start > isqrt(primeApprox) / 10)
   {
     // Count primes > start
     start = checkedAdd(start, 1);

From 5bf8fcee873b3f6d1465541cf7bc48ad2b5b80ba Mon Sep 17 00:00:00 2001
From: Kim Walisch <kim.walisch@gmail.com>
Date: Sun, 7 Jan 2024 21:11:26 +0100
Subject: [PATCH 12/22] Fix minimum primeApprox

---
 src/nthPrime.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/nthPrime.cpp b/src/nthPrime.cpp
index f6185ba42..fae603437 100644
--- a/src/nthPrime.cpp
+++ b/src/nthPrime.cpp
@@ -61,6 +61,7 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
     primeApprox = nthPrimeApprox(nApprox);
   }
 
+  primeApprox = std::max(primeApprox, start);
   if (primeApprox > 0)
     avgPrimeGap = ilog(primeApprox) + 2;
 
@@ -69,8 +70,7 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
   // countPrimes() and hence the initailization 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 &&
-      primeApprox - start > isqrt(primeApprox) / 10)
+  if (primeApprox - start > isqrt(primeApprox) / 10)
   {
     // Count primes > start
     start = checkedAdd(start, 1);

From 1ea8fe94d56c81f16e0e6ad4b44142ca6aab78c7 Mon Sep 17 00:00:00 2001
From: Kim Walisch <kim.walisch@gmail.com>
Date: Sun, 7 Jan 2024 21:28:41 +0100
Subject: [PATCH 13/22] Fix typo

---
 src/nthPrime.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/nthPrime.cpp b/src/nthPrime.cpp
index fae603437..69c42af68 100644
--- a/src/nthPrime.cpp
+++ b/src/nthPrime.cpp
@@ -67,7 +67,7 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
 
   // Only use multi-threading if the sieving distance is sufficiently
   // large. For small n this if statement also avoids calling
-  // countPrimes() and hence the initailization overhead of
+  // 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)

From 535ab5f2285c55599e87f5eb60566c54649565c7 Mon Sep 17 00:00:00 2001
From: Kim Walisch <kim.walisch@gmail.com>
Date: Sun, 7 Jan 2024 21:49:18 +0100
Subject: [PATCH 14/22] Fix issues

---
 src/nthPrime.cpp | 27 +++++++++++++++++++++------
 1 file changed, 21 insertions(+), 6 deletions(-)

diff --git a/src/nthPrime.cpp b/src/nthPrime.cpp
index 69c42af68..cea58e493 100644
--- a/src/nthPrime.cpp
+++ b/src/nthPrime.cpp
@@ -46,9 +46,9 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
 
   setStart(start);
   auto t1 = std::chrono::system_clock::now();
-  int64_t countApprox = 0;
-  uint64_t primeApprox = 0;
+  uint64_t primeApprox = start;
   uint64_t avgPrimeGap = 0;
+  int64_t countApprox = 0;
   uint64_t prime = 0;
 
   if (start == 0)
@@ -70,7 +70,9 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
   // 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)
+  if (primeApprox - start < isqrt(primeApprox) / 10)
+    primeApprox = start;
+  else
   {
     // Count primes > start
     start = checkedAdd(start, 1);
@@ -126,14 +128,27 @@ uint64_t PrimeSieve::negativeNthPrime(int64_t n, uint64_t start)
   uint64_t primeApprox = nthPrimeApprox(nApprox);
   uint64_t avgPrimeGap = 0;
   uint64_t prime = 0;
+  int64_t countApprox = 0;
 
-  // Count primes < start
-  start = checkedSub(start, 1);
   primeApprox = std::min(primeApprox, start);
-  int64_t countApprox = countPrimes(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
+  {
+    // Count primes < start
+    start = checkedSub(start, 1);
+    primeApprox = std::min(primeApprox, start);
+    countApprox = countPrimes(primeApprox, start);
+  }
+
   // 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)

From 5bae46c9e8e7f2e08b08355433524dff0b9bd09c Mon Sep 17 00:00:00 2001
From: Kim Walisch <kim.walisch@gmail.com>
Date: Mon, 8 Jan 2024 08:55:12 +0100
Subject: [PATCH 15/22] Refactor

---
 include/primesieve/pmath.hpp |  6 ---
 src/nthPrime.cpp             | 97 +++++++++++++++++-------------------
 2 files changed, 47 insertions(+), 56 deletions(-)

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 <typename T>
-inline int ilog(T x)
-{
-  return (int) std::log((double) x);
-}
-
 template <typename T>
 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<double> seconds = t2 - t1;

From 9275650ea971b819adf9d1eacce083806f5763fb Mon Sep 17 00:00:00 2001
From: Kim Walisch <kim.walisch@gmail.com>
Date: Mon, 8 Jan 2024 09:05:04 +0100
Subject: [PATCH 16/22] Improve bounds

---
 src/nthPrime.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/nthPrime.cpp b/src/nthPrime.cpp
index f43fff09b..4c1b68a81 100644
--- a/src/nthPrime.cpp
+++ b/src/nthPrime.cpp
@@ -86,7 +86,7 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
   if (countApprox < n)
   {
     start = checkedAdd(start, 1);
-    uint64_t dist = (n - countApprox) * avgPrimeGap(start);
+    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++)
@@ -94,7 +94,7 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
   }
   else // if (countApprox >= n)
   {
-    uint64_t dist = (countApprox - n) * avgPrimeGap(start);
+    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--)

From b777c81ea93103bb606587139b4f863b4bf9b50f Mon Sep 17 00:00:00 2001
From: Kim Walisch <kim.walisch@gmail.com>
Date: Mon, 8 Jan 2024 09:35:55 +0100
Subject: [PATCH 17/22] Add comment

---
 src/nthPrimeApprox.cpp | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/src/nthPrimeApprox.cpp b/src/nthPrimeApprox.cpp
index c65f6b816..b55f214b3 100644
--- a/src/nthPrimeApprox.cpp
+++ b/src/nthPrimeApprox.cpp
@@ -1,5 +1,18 @@
 ///
 /// @file  nthPrimeApprox.cpp
+///        This file contains implementations of the logarithmic
+///        integral and the Riemann R function which are very
+///        accurate approximations of PrimePi(x). Please note that
+///        most of this code has been copied from the primecount
+///        project: https://github.com/kimwalisch/primecount
+///
+///        Note that while the Riemann R function is extremely
+///        accurate it is much slower than other simpler PrimePi(x)
+///        approximations. When speed matters, e.g. for allocating
+///        a vector of primes, we avoid using the functions defined
+///        in this file. Currently, the functions defined in this
+///        file are only used in nthPrime.cpp where accuracy is of
+///        utmost importance.
 ///
 /// Copyright (C) 2024 Kim Walisch, <kim.walisch@gmail.com>
 ///

From b3927f40f84277d85f51fa6a05f7ad2c0bf64b1b Mon Sep 17 00:00:00 2001
From: Kim Walisch <kim.walisch@gmail.com>
Date: Mon, 8 Jan 2024 10:11:06 +0100
Subject: [PATCH 18/22] Add more tests

---
 test/Li.cpp      | 101 +++++++++++++++++++++++++++++++++++++
 test/Ri.cpp      | 127 +++++++++++++++++++++++++++++++++++++++++++++++
 test/moebius.cpp |  99 ++++++++++++++++++++++++++++++++++++
 3 files changed, 327 insertions(+)
 create mode 100644 test/Li.cpp
 create mode 100644 test/Ri.cpp
 create mode 100644 test/moebius.cpp

diff --git a/test/Li.cpp b/test/Li.cpp
new file mode 100644
index 000000000..f1f674a69
--- /dev/null
+++ b/test/Li.cpp
@@ -0,0 +1,101 @@
+///
+/// @file   Li.cpp
+/// @brief  Test the offset logarithmic integral function.
+///         Li(x) = li(x) - li(2)
+///
+/// Copyright (C) 2024 Kim Walisch, <kim.walisch@gmail.com>
+///
+/// This file is distributed under the BSD License. See the COPYING
+/// file in the top level directory.
+///
+
+#include <primesieve/nthPrimeApprox.hpp>
+
+#include <stdint.h>
+#include <iostream>
+#include <cmath>
+#include <cstdlib>
+#include <vector>
+
+using std::max;
+using std::size_t;
+using namespace primesieve;
+
+std::vector<uint64_t> Li_table =
+{
+                     5, // Li(10^1)
+                    29, // Li(10^2)
+                   176, // Li(10^3)
+                  1245, // Li(10^4)
+                  9628, // Li(10^5)
+                 78626, // Li(10^6)
+                664917, // Li(10^7)
+               5762208, // Li(10^8)
+              50849233, // Li(10^9)
+             455055613, // Li(10^10)
+          4118066399ll, // Li(10^11)
+         37607950279ll, // Li(10^12)
+        346065645809ll, // Li(10^13)
+       3204942065690ll, // Li(10^14)
+      29844571475286ll  // Li(10^15)
+};
+
+void check(bool OK)
+{
+  std::cout << "   " << (OK ? "OK" : "ERROR") << "\n";
+  if (!OK)
+    std::exit(1);
+}
+
+int main()
+{
+  uint64_t x = 1;
+  for (size_t i = 0; i < Li_table.size(); i++)
+  {
+    x *= 10;
+    std::cout << "Li(" << x << ") = " << Li(x);
+    check(Li(x) == Li_table[i]);
+  }
+
+  x = 1;
+  for (size_t i = 0; i < Li_table.size(); i++)
+  {
+    x *= 10;
+    std::cout << "Li_inverse(" << Li_table[i] << ") = " << Li_inverse(Li_table[i]);
+    check(Li_inverse(Li_table[i]) <= x &&
+          Li_inverse(Li_table[i] + 1) > x);
+  }
+
+  // Sanity checks for small values of Li(x)
+  for (x = 0; x < 300000; x++)
+  {
+    uint64_t lix = Li(x);
+    double logx = std::log(max((double) x, 2.0));
+
+    if ((x >= 11 && lix < x / logx) ||
+        (x >= 2  && lix > x * logx))
+    {
+      std::cout << "Li(" << x << ") = " << lix << "   ERROR" << std::endl;
+      std::exit(1);
+    }
+  }
+
+  // Sanity checks for small values of Li_inverse(x)
+  for (x = 2; x < 30000; x++)
+  {
+    uint64_t res = Li_inverse(x);
+    double logx = std::log((double) x);
+
+    if (res < x ||
+        (x >= 4 && res > x * logx * logx))
+    {
+      std::cout << "Li_inverse(" << x << ") = " << res << "   ERROR" << std::endl;
+      std::exit(1);
+    }
+  }
+
+  std::cout << std::endl;
+  std::cout << "All tests passed successfully!" << std::endl;
+
+  return 0;
+}
diff --git a/test/Ri.cpp b/test/Ri.cpp
new file mode 100644
index 000000000..ad71c321a
--- /dev/null
+++ b/test/Ri.cpp
@@ -0,0 +1,127 @@
+///
+/// @file   Ri.cpp
+/// @brief  Test the Riemann R function.
+///
+/// Copyright (C) 2024 Kim Walisch, <kim.walisch@gmail.com>
+///
+/// This file is distributed under the BSD License. See the COPYING
+/// file in the top level directory.
+///
+
+#include <primesieve/nthPrimeApprox.hpp>
+
+#include <stdint.h>
+#include <iostream>
+#include <cstdlib>
+#include <cmath>
+#include <vector>
+
+using std::max;
+using std::size_t;
+using namespace primesieve;
+
+std::vector<uint64_t> Ri_table =
+{
+                4, // Ri(10^1)
+               25, // Ri(10^2)
+              168, // Ri(10^3)
+             1226, // Ri(10^4)
+             9587, // Ri(10^5)
+            78527, // Ri(10^6)
+           664667, // Ri(10^7)
+          5761551, // Ri(10^8)
+         50847455, // Ri(10^9)
+        455050683, // Ri(10^10)
+     4118052494ll, // Ri(10^11)
+    37607910542ll, // Ri(10^12)
+   346065531065ll, // Ri(10^13)
+  3204941731601ll  // Ri(10^14)
+};
+
+void check(bool OK)
+{
+  std::cout << "   " << (OK ? "OK" : "ERROR") << "\n";
+  if (!OK)
+    std::exit(1);
+}
+
+int main()
+{
+  uint64_t x = 1;
+  for (size_t i = 0; i < Ri_table.size(); i++)
+  {
+    x *= 10;
+    std::cout << "Ri(" << x << ") = " << Ri(x);
+    check(Ri(x) == Ri_table[i]);
+  }
+
+  x = 1;
+  for (size_t i = 0; i < Ri_table.size(); i++)
+  {
+    x *= 10;
+    std::cout << "Ri_inverse(" << Ri_table[i] << ") = " << Ri_inverse(Ri_table[i]);
+    check(Ri_inverse(Ri_table[i]) < x &&
+          Ri_inverse(Ri_table[i] + 1) >= x);
+  }
+
+  // Sanity checks for tiny values of Ri(x)
+  for (x = 0; x < 10000; x++)
+  {
+    uint64_t rix = Ri(x);
+    double logx = std::log(max((double) x, 2.0));
+
+    if ((x >= 20 && rix < x / logx) ||
+        (x >= 2  && rix > x * logx))
+    {
+      std::cout << "Ri(" << x << ") = " << rix << "   ERROR" << std::endl;
+      std::exit(1);
+    }
+  }
+
+  // Sanity checks for small values of Ri(x)
+  for (; x < 100000; x += 101)
+  {
+    uint64_t rix = Ri(x);
+    double logx = std::log(max((double) x, 2.0));
+
+    if ((x >= 20 && rix < x / logx) ||
+        (x >= 2  && rix > x * logx))
+    {
+      std::cout << "Ri(" << x << ") = " << rix << "   ERROR" << std::endl;
+      std::exit(1);
+    }
+  }
+
+  // Sanity checks for tiny values of Ri_inverse(x)
+  for (x = 2; x < 1000; x++)
+  {
+    uint64_t res = Ri_inverse(x);
+    double logx = std::log((double) x);
+
+    if (res < x ||
+        (x >= 5 && res > x * logx * logx))
+    {
+      std::cout << "Ri_inverse(" << x << ") = " << res << "   ERROR" << std::endl;
+      std::exit(1);
+    }
+  }
+
+  // Sanity checks for small values of Ri_inverse(x)
+  for (; x < 100000; x += 101)
+  {
+    uint64_t res = Ri_inverse(x);
+    double logx = std::log((double) x);
+
+    if (res < x ||
+        (x >= 5 && res > x * logx * logx))
+    {
+      std::cout << "Ri_inverse(" << x << ") = " << res << "   ERROR" << std::endl;
+      std::exit(1);
+    }
+  }
+
+  std::cout << std::endl;
+  std::cout << "All tests passed successfully!" << std::endl;
+
+  return 0;
+}
diff --git a/test/moebius.cpp b/test/moebius.cpp
new file mode 100644
index 000000000..579ac0d97
--- /dev/null
+++ b/test/moebius.cpp
@@ -0,0 +1,99 @@
+///
+/// @file   moebius.cpp
+/// @brief  Test the generate_moebius(n) function
+/// @link   https://en.wikipedia.org/wiki/Moebius_function
+///
+/// Copyright (C) 2024 Kim Walisch, <kim.walisch@gmail.com>
+///
+/// This file is distributed under the BSD License. See the COPYING
+/// file in the top level directory.
+///
+
+#include <primesieve/nthPrimeApprox.hpp>
+
+#include <stdint.h>
+#include <iostream>
+#include <cstdlib>
+#include <vector>
+
+using std::size_t;
+using namespace primesieve;
+
+/// Values of mu(n) for the first 1000 integers
+/// https://oeis.org/A008683/b008683.txt
+std::vector<int> moebius =
+{
+  0, 1, -1, -1, 0, -1, 1, -1, 0, 0, 1, -1, 0, -1, 1, 1, 0, -1, 0, -1,
+  0, 1, 1, -1, 0, 0, 1, 0, 0, -1, -1, -1, 0, 1, 1, 1, 0, -1, 1, 1,
+  0, -1, -1, -1, 0, 0, 1, -1, 0, 0, 0, 1, 0, -1, 0, 1, 0, 1, 1, -1,
+  0, -1, 1, 0, 0, 1, -1, -1, 0, 1, -1, -1, 0, -1, 1, 0, 0, 1, -1, -1,
+  0, 0, 1, -1, 0, 1, 1, 1, 0, -1, 0, 1, 0, 1, 1, 1, 0, -1, 0, 0,
+  0, -1, -1, -1, 0, -1, 1, -1, 0, -1, -1, 1, 0, -1, -1, 1, 0, 0, 1, 1,
+  0, 0, 1, 1, 0, 0, 0, -1, 0, 1, -1, -1, 0, 1, 1, 0, 0, -1, -1, -1,
+  0, 1, 1, 1, 0, 1, 1, 0, 0, -1, 0, -1, 0, 0, -1, 1, 0, -1, 1, 1,
+  0, 1, 0, -1, 0, -1, 1, -1, 0, 0, -1, 0, 0, -1, -1, 0, 0, 1, 1, -1,
+  0, -1, -1, 1, 0, 1, -1, 1, 0, 0, -1, -1, 0, -1, 1, -1, 0, -1, 0, -1,
+  0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, -1, 0, 1, 1, 1, 0, 1, 1, 1,
+  0, 1, -1, -1, 0, 0, 1, -1, 0, -1, -1, -1, 0, -1, 0, 1, 0, 1, -1, -1,
+  0, -1, 0, 0, 0, 0, -1, 1, 0, 1, 0, -1, 0, 1, 1, -1, 0, -1, -1, 1,
+  0, 0, 1, -1, 0, 1, -1, 1, 0, -1, 0, -1, 0, -1, 1, 0, 0, -1, 1, 0,
+  0, -1, -1, -1, 0, -1, -1, 1, 0, 0, -1, 1, 0, -1, 0, 1, 0, 0, 1, 1,
+  0, 1, 1, 1, 0, 1, 0, -1, 0, 1, -1, -1, 0, -1, 1, 0, 0, -1, -1, 1,
+  0, 1, -1, 1, 0, 0, 1, 1, 0, 1, 1, -1, 0, 0, 1, 1, 0, -1, 0, 1,
+  0, 1, 0, 0, 0, -1, 1, -1, 0, -1, 0, 0, 0, -1, -1, 1, 0, -1, 1, -1,
+  0, 0, 1, 0, 0, 1, -1, -1, 0, 0, -1, 1, 0, -1, -1, 0, 0, 1, 0, -1,
+  0, 1, 1, -1, 0, -1, 1, 0, 0, -1, 1, 1, 0, 1, 1, 1, 0, -1, 1, -1,
+  0, -1, -1, 1, 0, 0, -1, 1, 0, -1, -1, 1, 0, 1, 0, 1, 0, 1, -1, -1,
+  0, -1, 1, 0, 0, 0, -1, 1, 0, -1, -1, -1, 0, -1, -1, -1, 0, 1, -1, -1,
+  0, 0, -1, -1, 0, 1, 1, 1, 0, -1, 0, 1, 0, 1, 1, -1, 0, -1, 1, 0,
+  0, -1, 1, -1, 0, -1, 1, -1, 0, 1, -1, 1, 0, 1, -1, 0, 0, 0, 1, -1,
+  0, 1, 1, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, -1, 0, 0, 1, -1, -1,
+  0, 1, 1, -1, 0, 1, -1, 0, 0, -1, 1, 1, 0, 0, 1, 1, 0, 1, -1, 1,
+  0, -1, 0, -1, 0, 0, 1, 1, 0, 0, -1, 0, 0, 1, -1, 1, 0, 1, 1, 0,
+  0, -1, 1, 1, 0, 1, 1, -1, 0, 0, 0, 1, 0, 1, 1, -1, 0, -1, 0, 1,
+  0, -1, 1, -1, 0, 1, 1, 0, 0, -1, 1, -1, 0, 1, -1, 0, 0, -1, 0, 1,
+  0, 1, -1, 1, 0, 0, 1, -1, 0, 1, -1, 1, 0, -1, 0, -1, 0, 1, -1, -1,
+  0, -1, -1, 0, 0, 0, -1, -1, 0, -1, -1, 1, 0, -1, 1, -1, 0, -1, -1, -1,
+  0, 0, 1, 1, 0, 0, 1, -1, 0, 1, 0, -1, 0, 1, 1, 1, 0, 0, -1, 0,
+  0, -1, -1, -1, 0, -1, -1, -1, 0, 1, 0, -1, 0, -1, -1, 1, 0, 0, -1, -1,
+  0, -1, 1, -1, 0, -1, 0, 1, 0, 1, -1, 1, 0, -1, 1, 0, 0, -1, -1, 1,
+  0, 1, -1, -1, 0, 1, 0, 1, 0, 1, 1, -1, 0, 0, 1, 1, 0, 1, 1, 1,
+  0, -1, 0, 1, 0, -1, 1, 1, 0, -1, -1, 0, 0, 1, 1, -1, 0, 1, 1, -1,
+  0, 1, 0, 1, 0, 0, 0, -1, 0, 0, -1, 1, 0, -1, 1, 0, 0, 1, 0, -1,
+  0, -1, -1, -1, 0, 1, 1, 0, 0, 1, 0, -1, 0, 1, -1, 1, 0, -1, 1, -1,
+  0, -1, -1, 1, 0, 0, 1, 1, 0, -1, 1, 1, 0, -1, 0, 0, 0, -1, 1, 1,
+  0, 1, -1, 0, 0, 1, -1, -1, 0, 1, -1, 1, 0, 1, 1, -1, 0, -1, 1, 1,
+  0, 0, 1, 1, 0, -1, -1, 1, 0, -1, 0, -1, 0, 1, -1, 1, 0, 1, 1, 0,
+  0, -1, -1, -1, 0, 0, -1, -1, 0, -1, -1, 1, 0, 0, -1, 1, 0, 0, 1, -1,
+  0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, -1, -1, 0, 0, -1, 1, -1,
+  0, -1, 1, -1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, -1, 0, 0, -1, 1, 1,
+  0, -1, 0, -1, 0, -1, 1, -1, 0, 1, -1, 0, 0, 1, -1, 1, 0, -1, 1, 1,
+  0, 1, -1, -1, 0, 1, -1, -1, 0, 0, 1, -1, 0, 1, 1, -1, 0, 1, 0, -1,
+  0, 1, 1, 1, 0, 0, 1, 0, 0, -1, 1, 0, 0, 1, 1, -1, 0, -1, -1, 1,
+  0, -1, -1, 1, 0, 0, -1, -1, 0, 1, 0, 1, 0, -1, 0, 1, 0, -1, 1, 1,
+  0, 0, -1, 0, 0, 1, 1, -1, 0, -1, -1, -1, 0, 1, 1, 0, 0, -1, -1, 1,
+  0, 0, 1, -1, 0, 1, -1, -1, 0, 1, 0, -1, 0, 1, -1, 1, 0, -1, 1, 0
+};
+
+void check(bool OK)
+{
+  std::cout << "   " << (OK ? "OK" : "ERROR") << "\n";
+  if (!OK)
+    std::exit(1);
+}
+
+int main()
+{
+  auto mu = generate_moebius(moebius.size() - 1);
+
+  for (size_t i = 1; i < mu.size(); i++)
+  {
+    std::cout << "mu(" << i << ") = " << mu[i];
+    check(mu[i] == moebius[i]);
+  }
+
+  std::cout << std::endl;
+  std::cout << "All tests passed successfully!" << std::endl;
+
+  return 0;
+}

From ac63cd42b66a310fa1b87c562e1ad65f53ffde9c Mon Sep 17 00:00:00 2001
From: Kim Walisch <kim.walisch@gmail.com>
Date: Mon, 8 Jan 2024 10:13:59 +0100
Subject: [PATCH 19/22] Update ChangeLog

---
 ChangeLog | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/ChangeLog b/ChangeLog
index 2fb20de57..ca8b49e00 100644
--- a/ChangeLog
+++ b/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.

From 896c7027953865c16422f8b213300248c1115be7 Mon Sep 17 00:00:00 2001
From: Kim Walisch <kim.walisch@gmail.com>
Date: Mon, 8 Jan 2024 10:24:20 +0100
Subject: [PATCH 20/22] Refactor

---
 src/nthPrime.cpp | 19 ++++++++++---------
 1 file changed, 10 insertions(+), 9 deletions(-)

diff --git a/src/nthPrime.cpp b/src/nthPrime.cpp
index 4c1b68a81..1d9548974 100644
--- a/src/nthPrime.cpp
+++ b/src/nthPrime.cpp
@@ -112,6 +112,7 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
   return prime;
 }
 
+/// Used for n < 0
 uint64_t PrimeSieve::negativeNthPrime(int64_t n, uint64_t start)
 {
   ASSERT(n < 0);
@@ -147,7 +148,15 @@ uint64_t PrimeSieve::negativeNthPrime(int64_t n, uint64_t start)
 
   // 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 (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();
+  }
+  else // if (countApprox < n)
   {
     start = checkedSub(start, 1);
     uint64_t dist = (countApprox - n) * avgPrimeGap(start);
@@ -160,14 +169,6 @@ uint64_t PrimeSieve::negativeNthPrime(int64_t n, uint64_t start)
         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<double> seconds = t2 - t1;

From 64a1365c054977ebe530885b3bc5996a14211de3 Mon Sep 17 00:00:00 2001
From: Kim Walisch <kim.walisch@gmail.com>
Date: Mon, 8 Jan 2024 12:21:20 +0100
Subject: [PATCH 21/22] Fix bounds

---
 src/nthPrime.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/nthPrime.cpp b/src/nthPrime.cpp
index 1d9548974..fb18268af 100644
--- a/src/nthPrime.cpp
+++ b/src/nthPrime.cpp
@@ -150,7 +150,7 @@ uint64_t PrimeSieve::negativeNthPrime(int64_t n, uint64_t start)
   // we simply iterate over the primes until we find it.
   if (countApprox >= n)
   {
-    uint64_t dist = (n - countApprox) * avgPrimeGap(start);
+    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--)
@@ -159,7 +159,7 @@ uint64_t PrimeSieve::negativeNthPrime(int64_t n, uint64_t start)
   else // if (countApprox < n)
   {
     start = checkedSub(start, 1);
-    uint64_t dist = (countApprox - n) * avgPrimeGap(start);
+    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++)

From 2b45426f4433c06feda13a2ba73339df0b80ec1d Mon Sep 17 00:00:00 2001
From: kimwalisch <kim.walisch@gmail.com>
Date: Mon, 8 Jan 2024 12:33:13 +0100
Subject: [PATCH 22/22] Improve error message

---
 src/nthPrime.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/nthPrime.cpp b/src/nthPrime.cpp
index fb18268af..e524220dc 100644
--- a/src/nthPrime.cpp
+++ b/src/nthPrime.cpp
@@ -101,7 +101,7 @@ uint64_t PrimeSieve::nthPrime(int64_t n, uint64_t start)
     {
       prime = iter.prev_prime();
       if_unlikely(prime == 0)
-        throw primesieve_error("nth_prime(n): n is too small, nth_prime(n) < 2!");
+        throw primesieve_error("nth_prime(n): invalid n, nth prime < 2 is impossible!");
     }
   }
 
@@ -166,7 +166,7 @@ uint64_t PrimeSieve::negativeNthPrime(int64_t n, uint64_t start)
     {
       prime = iter.prev_prime();
       if_unlikely(prime == 0)
-        throw primesieve_error("nth_prime(n): n is too small, nth_prime(n) < 2!");
+        throw primesieve_error("nth_prime(n): invalid n, nth prime < 2 is impossible!");
     }
   }