From deb5196909af4308714c70279c71c8a8c18e3162 Mon Sep 17 00:00:00 2001 From: kimwalisch Date: Fri, 9 Feb 2024 21:38:29 +0100 Subject: [PATCH] Rename Ri(x) to RiemannR(x) --- include/primesieve/nthPrimeApprox.hpp | 4 +- src/app/main.cpp | 4 +- src/nthPrimeApprox.cpp | 26 ++--- test/Ri.cpp | 138 -------------------------- test/Riemann_R.cpp | 138 ++++++++++++++++++++++++++ 5 files changed, 155 insertions(+), 155 deletions(-) delete mode 100644 test/Ri.cpp create mode 100644 test/Riemann_R.cpp diff --git a/include/primesieve/nthPrimeApprox.hpp b/include/primesieve/nthPrimeApprox.hpp index a8488e67d..8c612494d 100644 --- a/include/primesieve/nthPrimeApprox.hpp +++ b/include/primesieve/nthPrimeApprox.hpp @@ -12,8 +12,8 @@ namespace primesieve { -uint64_t Ri(uint64_t x); -uint64_t Ri_inverse(uint64_t x); +uint64_t RiemannR(uint64_t x); +uint64_t RiemannR_inverse(uint64_t x); uint64_t primePiApprox(uint64_t x); uint64_t nthPrimeApprox(uint64_t n); diff --git a/src/app/main.cpp b/src/app/main.cpp index 70d48d937..8bf59f9d6 100644 --- a/src/app/main.cpp +++ b/src/app/main.cpp @@ -138,9 +138,9 @@ int main(int argc, char* argv[]) if (opt.nthPrime) nthPrime(opt); else if (opt.RiemannR) - std::cout << Ri(opt.numbers[0]) << std::endl; + std::cout << RiemannR(opt.numbers[0]) << std::endl; else if (opt.RiemannR_inverse) - std::cout << Ri_inverse(opt.numbers[0]) << std::endl; + std::cout << RiemannR_inverse(opt.numbers[0]) << std::endl; else sieve(opt); } diff --git a/src/nthPrimeApprox.cpp b/src/nthPrimeApprox.cpp index efc529292..f85b2cf63 100644 --- a/src/nthPrimeApprox.cpp +++ b/src/nthPrimeApprox.cpp @@ -169,7 +169,7 @@ const primesieve::Array zetaInv = /// The calculation is done with the Gram series: /// RiemannR(x) = 1 + \sum_{k=1}^{∞} ln(x)^k / (zeta(k + 1) * k * k!) /// -long double Ri(long double x) +long double RiemannR(long double x) { if (x < 0.1) return 0; @@ -203,7 +203,7 @@ long double Ri(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 Ri_prime(long double x) +long double RiemannR_prime(long double x) { if (x < 0.1) return 0; @@ -236,15 +236,15 @@ long double Ri_prime(long double x) /// 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. +/// 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 - (Ri(zn) - x) / Ri'(zn) +/// zn+1 = zn - (RiemannR(zn) - x) / RiemannR'(zn) /// -long double Ri_inverse(long double x) +long double RiemannR_inverse(long double x) { if (x < 2) return 0; @@ -266,7 +266,7 @@ long double Ri_inverse(long double x) while (true) { - long double term = (Ri(t) - x) / Ri_prime(t); + long double term = (RiemannR(t) - x) / RiemannR_prime(t); // Not converging anymore if (std::abs(term) >= std::abs(old_term)) @@ -283,14 +283,14 @@ long double Ri_inverse(long double x) namespace primesieve { -uint64_t Ri(uint64_t x) +uint64_t RiemannR(uint64_t x) { - return (uint64_t) ::Ri((long double) x); + return (uint64_t) ::RiemannR((long double) x); } -uint64_t Ri_inverse(uint64_t x) +uint64_t RiemannR_inverse(uint64_t x) { - auto res = ::Ri_inverse((long double) x); + auto res = ::RiemannR_inverse((long double) x); // Prevent 64-bit integer overflow if (res > (long double) std::numeric_limits::max()) @@ -308,7 +308,7 @@ uint64_t Ri_inverse(uint64_t x) /// uint64_t primePiApprox(uint64_t x) { - return Ri(x); + return RiemannR(x); } /// nthPrimeApprox(n) is a very accurate approximation of the nth @@ -318,7 +318,7 @@ uint64_t primePiApprox(uint64_t x) /// uint64_t nthPrimeApprox(uint64_t n) { - return Ri_inverse(n); + return RiemannR_inverse(n); } } // namespace diff --git a/test/Ri.cpp b/test/Ri.cpp deleted file mode 100644 index 142189303..000000000 --- a/test/Ri.cpp +++ /dev/null @@ -1,138 +0,0 @@ -/// -/// @file Ri.cpp -/// @brief Test the Riemann R function. -/// -/// Copyright (C) 2024 Kim Walisch, -/// -/// This file is distributed under the BSD License. See the COPYING -/// file in the top level directory. -/// - -#include - -#include -#include -#include -#include -#include -#include - -using std::max; -using std::size_t; -using namespace primesieve; - -std::vector 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); - } - } - - { - uint64_t x = std::numeric_limits::max() / 10; - uint64_t res = Ri_inverse(x); - if (res != std::numeric_limits::max()) - { - std::cout << "Ri_inverse(" << x << ") != UINT64_MAX, failed to prevent integer overflow!" << std::endl; - std::exit(1); - } - } - - std::cout << std::endl; - std::cout << "All tests passed successfully!" << std::endl; - - return 0; -} diff --git a/test/Riemann_R.cpp b/test/Riemann_R.cpp new file mode 100644 index 000000000..910852747 --- /dev/null +++ b/test/Riemann_R.cpp @@ -0,0 +1,138 @@ +/// +/// @file Riemann_R.cpp +/// @brief Test the Riemann R function. +/// +/// Copyright (C) 2024 Kim Walisch, +/// +/// This file is distributed under the BSD License. See the COPYING +/// file in the top level directory. +/// + +#include + +#include +#include +#include +#include +#include +#include + +using std::max; +using std::size_t; +using namespace primesieve; + +std::vector RiemannR_table = +{ + 4, // RiemannR(10^1) + 25, // RiemannR(10^2) + 168, // RiemannR(10^3) + 1226, // RiemannR(10^4) + 9587, // RiemannR(10^5) + 78527, // RiemannR(10^6) + 664667, // RiemannR(10^7) + 5761551, // RiemannR(10^8) + 50847455, // RiemannR(10^9) + 455050683, // RiemannR(10^10) + 4118052494ll, // RiemannR(10^11) + 37607910542ll, // RiemannR(10^12) + 346065531065ll, // RiemannR(10^13) + 3204941731601ll // RiemannR(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 < RiemannR_table.size(); i++) + { + x *= 10; + std::cout << "RiemannR(" << x << ") = " << RiemannR(x); + check(RiemannR(x) == RiemannR_table[i]); + } + + x = 1; + for (size_t i = 0; i < RiemannR_table.size(); i++) + { + x *= 10; + std::cout << "RiemannR_inverse(" << RiemannR_table[i] << ") = " << RiemannR_inverse(RiemannR_table[i]); + check(RiemannR_inverse(RiemannR_table[i]) < x && + RiemannR_inverse(RiemannR_table[i] + 1) >= x); + } + + // Sanity checks for tiny values of RiemannR(x) + for (x = 0; x < 10000; x++) + { + uint64_t rix = RiemannR(x); + double logx = std::log(max((double) x, 2.0)); + + if ((x >= 20 && rix < x / logx) || + (x >= 2 && rix > x * logx)) + { + std::cout << "RiemannR(" << x << ") = " << rix << " ERROR" << std::endl; + std::exit(1); + } + } + + // Sanity checks for small values of RiemannR(x) + for (; x < 100000; x += 101) + { + uint64_t rix = RiemannR(x); + double logx = std::log(max((double) x, 2.0)); + + if ((x >= 20 && rix < x / logx) || + (x >= 2 && rix > x * logx)) + { + std::cout << "RiemannR(" << x << ") = " << rix << " ERROR" << std::endl; + std::exit(1); + } + } + + // Sanity checks for tiny values of RiemannR_inverse(x) + for (x = 2; x < 1000; x++) + { + uint64_t res = RiemannR_inverse(x); + double logx = std::log((double) x); + + if (res < x || + (x >= 5 && res > x * logx * logx)) + { + std::cout << "RiemannR_inverse(" << x << ") = " << res << " ERROR" << std::endl; + std::exit(1); + } + } + + // Sanity checks for small values of RiemannR_inverse(x) + for (; x < 100000; x += 101) + { + uint64_t res = RiemannR_inverse(x); + double logx = std::log((double) x); + + if (res < x || + (x >= 5 && res > x * logx * logx)) + { + std::cout << "RiemannR_inverse(" << x << ") = " << res << " ERROR" << std::endl; + std::exit(1); + } + } + + { + uint64_t x = std::numeric_limits::max() / 10; + uint64_t res = RiemannR_inverse(x); + if (res != std::numeric_limits::max()) + { + std::cout << "RiemannR_inverse(" << x << ") != UINT64_MAX, failed to prevent integer overflow!" << std::endl; + std::exit(1); + } + } + + std::cout << std::endl; + std::cout << "All tests passed successfully!" << std::endl; + + return 0; +}