Faster Convergence for the Smaller Values of the Riemann Zeta Function

Question

I have a C++ program that uses the equation

$$\zeta(s)=\sum_{n=1}^\infty\frac{1}{n^s}$$

to calculate the Riemann zeta function.

This equation converges fast for larger values, like 183, but converges much slower for smaller values, like 2. For example, calculating the value of $\zeta(2)$ took an hour to be accurate to 5 digits, but one second for $\zeta(183)$ to be accurate to 100 digits.

Are there any equations for calculating the Riemann zeta function that are faster for calculating smaller values?

Because I am coding in C++, I cannot use $\int$ (without implementing external libraries, which is not really an option here).

Answer 1

There is a derivation of $\displaystyle(1-2^{1-s})\zeta(s) = \eta(s) = \sum_{n=0}^\infty \frac{1}{2^{n+1}} \sum_{k=0}^n {n \choose k} \frac {(-1)^{k}}{(k+1)^s}$ converging in $\mathcal{O}(2^{-N})$ for every $s \in \mathbb{C}$

• $\displaystyle\sum_{k=0}^\infty x^k = \frac{1}{1-x} =\frac{1/2}{1-\frac{1+x}{2}}= \sum_{n=0}^\infty \frac{(1+x)^k}{2^{n+1}}$ $\displaystyle = \sum_{n=0}^\infty \frac{1}{2^{n+1}}\sum_{k=0}^n {n \choose k}x^k = \sum_{k=0}^\infty x^k \sum_{n=k}^\infty {n \choose k} \frac{1}{2^{n+1}}$

  identifying the coefficients : $\displaystyle \sum_{n=k}^\infty \frac{1}{2^{n+1}}{n \choose k}=1$

• For $\text{Re}(s) > 1$ where everything converges absolutely : $$\sum_{n=0}^\infty \frac{1}{2^{n+1}} \sum_{k=0}^n {n \choose k} \frac {(-1)^{k}}{(k+1)^s} =\sum_{k=0}^\infty \frac{(-1)^{k}}{(k+1)^s} \sum_{n=k}^\infty \frac{1}{2^{n+1}} {n \choose k} = \sum_{k=0}^\infty \frac{(-1)^k}{(k+1)^s} = \eta(s)$$

• The hard part is to show that for any compact $K\subset \mathbb{C}$, there is a constant $\alpha$ such that for every $n$ and every $s\in K$ : $\ \ |\sum_{k=0}^n {n \choose k} \frac{(-1)^{k}}{(k+1)^s}|< \alpha$, by noting this is the $n$th forward difference $\Delta^n(0)$ of the sequence $\left\{ \frac{(-1)^k}{(k+1)^s}\right\}_{k \in \mathbb{N}}$, so that the series converges in $\mathcal{O}(2^{-N})$ for every $s \in \mathbb{C}$ and is entire, i.e. there is a continuous function $\alpha: \mathbb{C} \to \mathbb{R}^+$ such that $$\forall s \in \mathbb{C}, \qquad \left|\eta(s)- \sum_{n=0}^N \frac{1}{2^{n+1}} \sum_{k=0}^n (-1)^{k} {n \choose k} \frac {1}{(k+1)^s}\right| < \alpha(s) 2^{-N}$$

Answer 2

As far as the speed of the algorithm and its implementation go, you probably won't find much better than implemented by the GNU scientific library. See the file zeta.c, line 776. The source is sometimes interspersed by relevant references to Abramowicz and Stegun and other sources. Update: See user1952009's comment below for the maths that's used there!

Albeit I understand that you prefer writing your own implementation, please only do so if that's a necessary exercise. It's a lot of work to do things effectively and keep all cases covered at the same time. It's easier to take the file from the GSL if you don't want the whole library (linking against which would lead to the necessity of including BLAS as well), it's designed to be separable into independent modules. You can also start exploring by modifying it to better suit your needs, in a situation similar to what you describe I once managed to shorten the code for some elliptic integrals to 5 lines or so of branchless code, executable on a GPU.

If you need some special additions like complex arguments or arbitrary precision, Pari/GP (ζ source) is similarly the place to go. But that's harder to understand and heavily relies on the stack model used in that framework. I think they use Euler's identity with precomputed primes and use a smart trick utilizing Horner's scheme to reduce the number of operations. I once studied the internals of that but just forgot, I'll update my answer if I decrypt the formula from the code.

Update: I think Pari uses a literal implementation of the algorithm described here. But they have a precomputed list of both primes and Bernoulli numbers so that makes it easier.

First, the optimum bounds $N$ and $M$ are computed from the required precision. Then $S$ (a partial sum of $k^{-s}$) is evaluated using $$\sum_{k=1}^N \frac1{k^s} = \sum_{k=1\atop k\ \mathrm{odd}}^N \frac1{k^s} + \frac1{2^s} \sum_{k=1}^{\lfloor N/2\rfloor} \frac1{k^s} = \sum_{k=1\atop k\ \mathrm{odd}}^N \frac1{k^s} + \frac1{2^s} \sum_{k=1\atop k\ \mathrm{odd}}^{\lfloor N/2\rfloor} \frac1{k^s} + \frac1{(2^s)^2} \sum_{k=1}^{\lfloor N/4\rfloor} \frac1{k^s} = \cdots$$ until the sum becomes trivial ($O(\log N)$ steps) using Horner scheme. To evaluate the finite sum over odd $k$s, they use prime factorization and a precomputed table of $\frac1{p^{ms}}$ for $p$ prime. This is viable as long as $N$ is moderately small. One can replace that by a finite version of the sum taken exactly from your former approach, just terminating much earlier.

$I$ and $T$ (the tail for $n > N$) are calculated together, evaluating the factorials and powers iteratively and using known $B_{2n}$. The error term $R$ is dropped.

Obviously choosing a good $(M,N)$ is crucial. I spent way too much time on the above to look up what they are, but I'm sure there's some freedom in how you pick them. You can choose $M$ small so you don't need very many $B_n$, and then hardcode those.