Faster and simpler RiemannR_inverse(x)

ChangeLog

@@ -1,8 +1,9 @@
-Changes in version 12.2, 18/03/2024
+Changes in version 12.2, 19/03/2024
 ===================================
 
 * RiemannR.cpp: Fix infinite loop on Linux i386,
   see https://github.com/kimwalisch/primecount/issues/66.
+* RiemannR.cpp: Faster and simpler RiemannR_inverse(x).
 * test/Riemann_R.cpp: Add more tests.
 
 Changes in version 12.1, 09/03/2024

src/RiemannR.cpp

@@ -173,22 +173,21 @@ const primesieve::Array<long double, 128> zeta =
 template <typename T>
 T initialNthPrimeApprox(T x)
 {
-  if (x < 2)
+  if (x < 1)
     return 0;
+  else if (x >= 1 && x < 2)
+    return 2;
+  else if (x >= 2 && x < 3)
+    return 3;
 
   T logx = std::log(x);
-  T t = logx;
+  T loglogx = std::log(logx);
+  T t = logx + 0.5 * loglogx;
 
-  if (x > /* e = */ 2.719)
-  {
-    T loglogx = std::log(logx);
-    t += 0.5 * loglogx;
-
-    if (x > 1600)
-      t += 0.5 * loglogx - 1.0 + (loglogx - 2.0) / logx;
-    if (x > 1200000)
-      t -= (loglogx * loglogx - 6.0 * loglogx + 11.0) / (2.0 * logx * logx);
-  }
+  if (x > 1600)
+    t += 0.5 * loglogx - 1.0 + (loglogx - 2.0) / logx;
+  if (x > 1200000)
+    t -= (loglogx * loglogx - 6.0 * loglogx + 11.0) / (2.0 * logx * logx);
 
   return x * t;
 }
@@ -232,81 +231,31 @@ T RiemannR(T x)
   return sum;
 }
 
-/// Calculate the derivative of the Riemann R function.
-/// RiemannR'(x) = 1/x * \sum_{k=1}^{∞} ln(x)^(k-1) / (zeta(k + 1) * k!)
-///
-template <typename T>
-T RiemannR_prime(T x)
-{
-  if (x < 0.1)
-    return 0;
-
-  T epsilon = std::numeric_limits<T>::epsilon();
-
-  // RiemannR_prime(1) = NaN.
-  // Hence we return RiemannR_prime(1.0000000000000001).
-  // Required because: sum / log(1) = 0 / 0.
-  if (std::abs(x - 1.0) <= epsilon)
-    return (T) 0.60792710185402643042L;
-
-  T sum = 0;
-  T term = 1;
-  T logx = std::log(x);
-
-  // The condition k < ITERS is required in case the computation
-  // does not converge. This happened on Linux i386 where
-  // the precision of the libc math functions is very limited.
-  for (unsigned k = 1; k < 1000; k++)
-  {
-    term *= logx / k;
-    T old_sum = sum;
-
-    if (k + 1 < zeta.size())
-      sum += term / T(zeta[k + 1]);
-    else
-      // For k >= 127, approximate zeta(k + 1) by 1
-      sum += term;
-
-    // Not converging anymore
-    if (std::abs(sum - old_sum) <= epsilon)
-      break;
-  }
-
-  return sum / (x * logx);
-}
-
 /// Calculate the inverse Riemann R function which is a very
 /// accurate approximation of the nth prime.
-/// This implementation computes RiemannR^-1(x) as the zero of the
-/// function f(z) = RiemannR(z) - x using the Newton–Raphson method.
-/// https://math.stackexchange.com/a/853192
-///
-/// Newton–Raphson method:
-/// zn+1 = zn - (f(zn) / f'(zn)).
-/// zn+1 = zn - (RiemannR(zn) - x) / RiemannR'(zn)
+/// This implementation computes RiemannR^-1(x) = t as the zero of the
+/// function f(t) = RiemannR(t) - x using the Newton–Raphson method.
+/// https://en.wikipedia.org/wiki/Newton%27s_method
 ///
 template <typename T>
 T RiemannR_inverse(T x)
 {
-  if (x < 2)
-    return 0;
-
   T t = initialNthPrimeApprox(x);
   T old_term = std::numeric_limits<T>::infinity();
 
+  if (x < 3)
+    return t;
+
   // The condition i < ITERS is required in case the computation
   // does not converge. This happened on Linux i386 where
   // the precision of the libc math functions is very limited.
   for (int i = 0; i < 100; i++)
   {
-    T term;
-
-    if (x < 1e10)
-      // Converges faster for small x
-      term = (RiemannR(t) - x) / RiemannR_prime(t);
-    else
-      // Converges faster for large x
-      term = (RiemannR(t) - x) * std::log(t);
+    // term = f(t) / f'(t)
+    // f(t) = RiemannR(t) - x
+    // RiemannR(t) ~ li(t), with li'(t) = 1 / log(t)
+    // term = (RiemannR(t) - x) / li'(t) = (RiemannR(t) - x) * log(t)
+    T term = (RiemannR(t) - x) * std::log(t);
 
     // Not converging anymore
     if (std::abs(term) >= std::abs(old_term))