Even faster ζ(2n) calculation!

In the previous post I benchmarked several algorithms for computing $ζ(2n)$. All of them where, at some point, an improvement over the implementation in mpfr. At very high $n$ nothing seemed to come even near the series summation with recycled powers method, but this method is impossibly slow at the start. So we used Harvey’s method to bridge the gap.

I originally tried the von Staudt-Clausen method with the series summation, but it preformed no different from the version with mpfr_zeta_ui. Using von Staudt-Clausen with Harvey’s method is pointless, since it already uses it internally. The last option, von Staudt-Clausen with power table recycling series is also impossible, the number of terms and precision required would increase with $n$, so the whole table would have to be regenerated at every step.

And then the obvious struck me: the powtable method can be combined with von Staudt-Clausen if $n$ is decreasing! Evaluating the first terms of the Khinchin calculation backwards is not too hard and a quick test showed that von Staudt-Clausen has a considerable performance boost over Harvey’s method.

A nice consequence of this choice is that all the code is currently developed by ourselves, excluding the GMP, MPFR and C++ libraries, but those are mainly for simple arithmetic operations (simple as in the operations, not the implementation!) and input/output routines.

More efficient von Staudt-Clausen

The low precison $B_n$ was currently evaluated at high precision. But it only has to be accurate to a few digits after the decimal (binary?) point since we will be rounding it to integers. Let’s estimate its magnitude:

$$ \left\vert B_n \right\vert = \zeta(n)\frac{2n!}{(2 \pi)^n} < \left( 1 + \frac{2}{2^{n}}\right)\frac{2 \mathrm{e}^{n \ln n - n}}{(2 \pi)^n} $$

Now the number of bits required in $B_n$ is

$$ \log_2 \left( 1 + \frac{2}{2^{n}}\right) + 1 + (n \ln n - n) \log_2 \mathrm{e} - n \log_2 2 \pi $$

With a safety margin of 39 bits and an upper bound of 2 on the first log₂ this becomes:

$$ 42 + (n \ln n - n) \log_2 \mathrm{e} - n \log_2 2 \pi $$

Less bits in zeta

In the last half of the $n$ the $\zeta(2n)-1$ consists of a bunch of zeroes, a one, a smaller bunch of zeroes and then some interesting bits. The second set of zeroes takes up valuable computation time which can be eliminated by taking out the $\frac{1}{2^{2n}}$ term and implementing it separately using bit shifting. A nice consequence is that after a certain n there is no zeta calculation necessary anymore!

Source code



Red linempfr_zeta_ui
Blue linempfr_zeta_ui with our von Staudt-Clausen
Purple lineDarvid Harvey’s bernoulli
Green lineSeries summation without recycling
Black lineSeries summation with power table
Orange linePower table + von Staudt-Clausen

Benchmark for ζ(2n) with 34 kb accuracy

Bench­mark for ζ(2n) with 340 kb ac­cu­ra­cy

Remco Bloemen
Math & Engineering