diff --git a/src/RiemannR.cpp b/src/RiemannR.cpp
index 83e93aebe..4370abb14 100644
--- a/src/RiemannR.cpp
+++ b/src/RiemannR.cpp
@@ -170,17 +170,18 @@ const primesieve::Array<long double, 128> zeta =
 /// Rendus Hebdomadaires des Séances de l'Académie des Sciences. 119: 848–849.
 /// https://en.wikipedia.org/wiki/Prime_number_theorem#Approximations_for_the_nth_prime_number
 ///
-double initialNthPrimeApprox(double x)
+template <typename T>
+T initialNthPrimeApprox(T x)
 {
   if (x < 2)
     return 0;
 
-  double logx = std::log(x);
-  double t = logx;
+  T logx = std::log(x);
+  T t = logx;
 
   if (x > /* e = */ 2.719)
   {
-    double loglogx = std::log(logx);
+    T loglogx = std::log(logx);
     t += 0.5 * loglogx;
 
     if (x > 1600)
@@ -192,24 +193,20 @@ double initialNthPrimeApprox(double x)
   return x * t;
 }
 
-/// Calculate the derivative of the Riemann R function.
-/// RiemannR'(x) = 1/x * \sum_{k=1}^{∞} ln(x)^(k-1) / (zeta(k + 1) * k!)
+/// Calculate the Riemann R function which is a very accurate
+/// approximation of the number of primes below x.
+/// http://mathworld.wolfram.com/RiemannPrimeCountingFunction.html
+/// The calculation is done with the Gram series:
+/// RiemannR(x) = 1 + \sum_{k=1}^{∞} ln(x)^k / (zeta(k + 1) * k * k!)
 ///
 template <typename T>
-T RiemannR_prime(T x)
+T RiemannR(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 sum = 1;
   T term = 1;
   T logx = std::log(x);
 
@@ -222,33 +219,37 @@ T RiemannR_prime(T x)
     T old_sum = sum;
 
     if (k + 1 < zeta.size())
-      sum += term / T(zeta[k + 1]);
+      sum += term / (T(zeta[k + 1]) * k);
     else
       // For k >= 127, approximate zeta(k + 1) by 1
-      sum += term;
+      sum += term / k;
 
     // Not converging anymore
     if (std::abs(sum - old_sum) <= epsilon)
       break;
   }
 
-  return sum / (x * logx);
+  return sum;
 }
 
-/// Calculate the Riemann R function which is a very accurate
-/// approximation of the number of primes below x.
-/// http://mathworld.wolfram.com/RiemannPrimeCountingFunction.html
-/// The calculation is done with the Gram series:
-/// RiemannR(x) = 1 + \sum_{k=1}^{∞} ln(x)^k / (zeta(k + 1) * k * k!)
+/// 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(T x)
+T RiemannR_prime(T x)
 {
   if (x < 0.1)
     return 0;
 
   T epsilon = std::numeric_limits<T>::epsilon();
-  T sum = 1;
+
+  // 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);
 
@@ -261,17 +262,17 @@ T RiemannR(T x)
     T old_sum = sum;
 
     if (k + 1 < zeta.size())
-      sum += term / (T(zeta[k + 1]) * k);
+      sum += term / T(zeta[k + 1]);
     else
       // For k >= 127, approximate zeta(k + 1) by 1
-      sum += term / k;
+      sum += term;
 
     // Not converging anymore
     if (std::abs(sum - old_sum) <= epsilon)
       break;
   }
 
-  return sum;
+  return sum / (x * logx);
 }
 
 /// Calculate the inverse Riemann R function which is a very
@@ -290,7 +291,7 @@ T RiemannR_inverse(T x)
   if (x < 2)
     return 0;
 
-  T t = (T) initialNthPrimeApprox((double) x);
+  T t = initialNthPrimeApprox(x);
   T old_term = std::numeric_limits<T>::infinity();
 
   // The condition i < ITERS is required in case the computation