Update to latest libprimesieve

lib/primesieve/CMakeLists.txt

@@ -1,7 +1,7 @@
 cmake_minimum_required(VERSION 3.4...3.27)
 project(primesieve CXX)
-set(PRIMESIEVE_VERSION "12.1")
-set(PRIMESIEVE_SOVERSION "12.1.0")
+set(PRIMESIEVE_VERSION "12.2")
+set(PRIMESIEVE_SOVERSION "12.2.0")
 
 # Build options ######################################################
 

lib/primesieve/ChangeLog

@@ -1,3 +1,9 @@
+Changes in version 12.2, 18/03/2024
+===================================
+
+* RiemannR.cpp: Fix infinite loop on Linux i386,
+  see https://github.com/kimwalisch/primecount/issues/66.
+
 Changes in version 12.1, 09/03/2024
 ===================================
 

lib/primesieve/include/primesieve.h

@@ -15,9 +15,9 @@
 #ifndef PRIMESIEVE_H
 #define PRIMESIEVE_H
 
-#define PRIMESIEVE_VERSION "12.1"
+#define PRIMESIEVE_VERSION "12.2"
 #define PRIMESIEVE_VERSION_MAJOR 12
-#define PRIMESIEVE_VERSION_MINOR 1
+#define PRIMESIEVE_VERSION_MINOR 2
 
 #include <primesieve/iterator.h>
 

lib/primesieve/include/primesieve.hpp

@@ -13,9 +13,9 @@
 #ifndef PRIMESIEVE_HPP
 #define PRIMESIEVE_HPP
 
-#define PRIMESIEVE_VERSION "12.1"
+#define PRIMESIEVE_VERSION "12.2"
 #define PRIMESIEVE_VERSION_MAJOR 12
-#define PRIMESIEVE_VERSION_MINOR 1
+#define PRIMESIEVE_VERSION_MINOR 2
 
 #include <primesieve/iterator.hpp>
 #include <primesieve/primesieve_error.hpp>

lib/primesieve/src/RiemannR.cpp

@@ -170,17 +170,17 @@ 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
 ///
-long double initialNthPrimeApprox(long double x)
+double initialNthPrimeApprox(double x)
 {
   if (x < 2)
     return 0;
 
-  long double logx = std::log(x);
-  long double t = logx;
+  double logx = std::log(x);
+  double t = logx;
 
   if (x > /* e = */ 2.719)
   {
-    long double loglogx = std::log(logx);
+    double loglogx = std::log(logx);
     t += 0.5 * loglogx;
 
     if (x > 1600)
@@ -195,70 +195,80 @@ long double initialNthPrimeApprox(long double x)
 /// Calculate the derivative of the Riemann R function.
 /// RiemannR'(x) = 1/x * \sum_{k=1}^{∞} ln(x)^(k-1) / (zeta(k + 1) * k!)
 ///
-long double RiemannR_prime(long double x)
+template <typename T>
+T RiemannR_prime(T x)
 {
   if (x < 0.1)
     return 0;
 
-  long double epsilon = std::numeric_limits<long double>::epsilon();
+  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 0.60792710185402643042L;
+  if (std::abs(x - 1.0) <= epsilon)
+    return (T) 0.60792710185402643042L;
 
-  long double sum = 0;
-  long double old_sum = -1;
-  long double term = 1;
-  long double logx = std::log(x);
+  T sum = 0;
+  T term = 1;
+  T logx = std::log(x);
 
-  for (unsigned k = 1; std::abs(old_sum - sum) >= epsilon; k++)
+  // 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;
-    old_sum = sum;
+    T old_sum = sum;
 
     if (k + 1 < zeta.size())
-      sum += term / zeta[k + 1];
+      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);
 }
 
-} // namespace
-
-namespace primesieve {
-
 /// 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!)
 ///
-long double RiemannR(long double x)
+template <typename T>
+T RiemannR(T x)
 {
   if (x < 0.1)
     return 0;
 
-  long double epsilon = std::numeric_limits<long double>::epsilon();
-  long double sum = 1;
-  long double old_sum = -1;
-  long double term = 1;
-  long double logx = std::log(x);
+  T epsilon = std::numeric_limits<T>::epsilon();
+  T sum = 1;
+  T term = 1;
+  T logx = std::log(x);
 
-  for (unsigned k = 1; std::abs(old_sum - sum) >= epsilon; k++)
+  // 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;
-    old_sum = sum;
+    T old_sum = sum;
 
     if (k + 1 < zeta.size())
-      sum += term / (zeta[k + 1] * k);
+      sum += term / (T(zeta[k + 1]) * k);
     else
       // For k >= 127, approximate zeta(k + 1) by 1
       sum += term / k;
+
+    // Not converging anymore
+    if (std::abs(sum - old_sum) <= epsilon)
+      break;
   }
 
   return sum;
@@ -274,17 +284,21 @@ long double RiemannR(long double x)
 /// zn+1 = zn - (f(zn) / f'(zn)).
 /// zn+1 = zn - (RiemannR(zn) - x) / RiemannR'(zn)
 ///
-long double RiemannR_inverse(long double x)
+template <typename T>
+T RiemannR_inverse(T x)
 {
   if (x < 2)
     return 0;
 
-  long double t = initialNthPrimeApprox(x);
-  long double old_term = std::numeric_limits<long double>::infinity();
+  T t = (T) initialNthPrimeApprox((double) x);
+  T old_term = std::numeric_limits<T>::infinity();
 
-  while (true)
+  // 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++)
   {
-    long double term = (RiemannR(t) - x) / RiemannR_prime(t);
+    T term = (RiemannR(t) - x) / RiemannR_prime(t);
 
     // Not converging anymore
     if (std::abs(term) >= std::abs(old_term))
@@ -297,6 +311,26 @@ long double RiemannR_inverse(long double x)
   return t;
 }
 
+} // namespace
+
+namespace primesieve {
+
+long double RiemannR(long double x)
+{
+  if (x > 1e8)
+    return ::RiemannR(x);
+  else
+    return ::RiemannR((double) x);
+}
+
+long double RiemannR_inverse(long double x)
+{
+  if (x > 1e8)
+    return ::RiemannR_inverse(x);
+  else
+    return ::RiemannR_inverse((double) x);
+}
+
 /// primePiApprox(x) is a very accurate approximation of PrimePi(x)
 /// with |PrimePi(x) - primePiApprox(x)| < sqrt(x).
 /// primePiApprox(x) is currently only used in nthPrime.cpp where