Use a faster series for calculating the Riemann R function

[nipzu]

Use the Gram series
$$R(x) = 1 + \sum_{k=1}^\infty \frac{(\ln x)^k}{\zeta(k + 1) \cdot k \cdot k! }$$
to calculate the Riemann R function. I used a precalculated table for $\zeta(k+1)$. Since this is (almost?) as fast as the series for $\text{li}(x)$, I took the liberty to just use $R(x)$ everywhere.

I also updated the Newton iteration to use $R'(x)$ instead of the approximation $\frac1{\ln x}$. This guarantees quadratic convergence and reduces the number of iterations needed, so it should outweigh the cost of calculating $R'(x)$. The derivative $R'(x)$ can be calculated as
$$R'(x) = \frac1{x \ln x} \sum_{k=1}^\infty \frac{(\ln x)^k}{\zeta(k + 1) \cdot k!}.$$
For the initial guess, I used the Cesaro formula (see https://en.wikipedia.org/wiki/Prime_number_theorem#Approximations_for_the_nth_prime_number) with some modifications that make it work better with smaller values of $x$ (see https://www.desmos.com/calculator/4jm9ulxsvh).

[kimwalisch]

Thanks, I am definitely interested in a faster R(x)!

In fact R(x) is even more important in my other primecount project. I see you have a zetaInv lookup table with a "precision of 128 bits". Does this mean R(x) can be computed fairly precisely (similar to my current implementation) up to 2^128? (I am mainly interested in understanding the upper limit of your implementation)

It will take me some time to study your implementation and throughly test it before I can eventually merge it. Note that I will take care to make the required changes in my other primecount project.

[nipzu]

In short, the value of $R(x)$ can be computed as accurately as the floating point type allows (as long as the type has at most 128 mantissa bits) for any value of $x$.

The values of zetaInv have "a precision of 128" bits in the sense that if some platform were to have a long double with a mantissa of 128 bits, the constants have enough decimal digits to specify the exact value of this 128-bit long double that is closest to $1/\zeta(k)$. This may be a bit overkill since long double has usually less than 128 bits of precision.

Let's see how precisely we can calculate $R(x)$ with these values. Consider an approximation
$$R_\text{approx}(x) = 1 + \sum_{k=1}^\infty \frac{(\ln x)^k}{k \cdot k!} \cdot \text{zetaInv}[k]$$
of $R(x)$ using the values of zetaInv rounded to those floating point representations. Note that "zetaInv[k] = 1" if k >= 128. Since $0.5 < 1/\zeta(k+1) \le 1$, the rounding error is bounded by $|\text{zetaInv}[k] - 1/\zeta(k+1)| \le 2^{-129}$.
Then we can calculate a bound on the absolute error:
$$|R_\text{approx}(x) - R(x)| = \left|\sum_{k=1}^\infty \frac{(\ln x)^k}{k \cdot k!} \cdot (\text{zetaInv}[k] - 1/\zeta(k+1))\right|\le \sum_{k=1}^\infty \frac{(\ln x)^k}{k \cdot k!} \cdot |\text{zetaInv}[k] - 1/\zeta(k+1)|,$$
which can be furher bounded by
$$\sum_{k=1}^\infty \frac{(\ln x)^k}{k \cdot k!} \cdot |\text{zetaInv}[k] - 1/\zeta(k+1)| \le \sum_{k=1}^\infty \frac{(\ln x)^k}{k!} \cdot 2^{-129} = 2^{-129} \exp(\ln(x)) = x / 2^{129}.$$
Thus, the absolute error $|R_\text{approx}(x) - R(x)|$ would be less than ${x}/{2^{129}}$ for any value of $x$.

Of course, the above ignores any other floating point rounding that would happen in practice. The series would need to be truncated, but it turns out that summing the first $\max(n, 2 e \ln x )$ terms (in practice probably fewer) is enough to reach a maximum error of $2^{-n}$.

Since $\zeta(k) \approx 1 + 1/2^k$ for large $k$, it turns out that if we want to calculate $1/\zeta(k+1)$ to $n$ bits of precision, we only need to precompute about the first $n$ values. The rest can be safely set to $1$.

In conclusion, if you wanted to calculate $R(x)$ with a maximum absolute error of $x/2^{n}$, you would need:

• a floating point type with at least $n$ mantissa bits,
• a precomputed table of about the first $n$ values of $1/\zeta(k+1)$ computed to $n$ bits,
• to evaluate the sum up to $\max(n, 2 e \ln x )$ terms.

[nipzu]

You can see a derivation of the Gram series here: https://math.stackexchange.com/questions/1506980/could-someone-explain-how-the-gram-series-relates-to-riemanns-function

[kimwalisch]

@nipzu I have now merged your code into the master branch #145. Your code is really of excellent quality 😊

I made a few changes on top of your code:

• I renamed nthPrimeApprox.cpp to RiemannR.cpp.
• I added -R and --R-inverse command-line options to the primesieve application, so that it is easier to play around with those new functions.
• I simplified your code e.g. instead of a zetaInv lookup with 1/Zeta[k+1] values I now have a zeta lookup table with Zeta[k] values. In my tests this did not deteriorate the performance and accuracy. I also merged your 2 for loops (k < 128 & k >= 128) into a single loop.
• I have added more documentation.