Skip to content

Commit

Permalink
Rename Ri(x) to RiemannR(x)
Browse files Browse the repository at this point in the history
  • Loading branch information
kimwalisch committed Feb 9, 2024
1 parent aae96b3 commit deb5196
Show file tree
Hide file tree
Showing 5 changed files with 155 additions and 155 deletions.
4 changes: 2 additions & 2 deletions include/primesieve/nthPrimeApprox.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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);

Expand Down
4 changes: 2 additions & 2 deletions src/app/main.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}
Expand Down
26 changes: 13 additions & 13 deletions src/nthPrimeApprox.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -169,7 +169,7 @@ const primesieve::Array<long double, 128> 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;
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand All @@ -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))
Expand All @@ -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<uint64_t>::max())
Expand All @@ -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
Expand All @@ -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
138 changes: 0 additions & 138 deletions test/Ri.cpp

This file was deleted.

138 changes: 138 additions & 0 deletions test/Riemann_R.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,138 @@
///
/// @file Riemann_R.cpp
/// @brief Test the Riemann R function.
///
/// Copyright (C) 2024 Kim Walisch, <[email protected]>
///
/// 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 <limits>
#include <vector>

using std::max;
using std::size_t;
using namespace primesieve;

std::vector<uint64_t> 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<uint64_t>::max() / 10;
uint64_t res = RiemannR_inverse(x);
if (res != std::numeric_limits<uint64_t>::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;
}

0 comments on commit deb5196

Please sign in to comment.