Even faster ζ(2n) calculation!

Remco Bloemen

2009-11-10, last updated 2014-02-27

In the previous post I benchmarked several algorithms for computing ζ(2n)\zeta (2n). All of them where, at some point, an improvement over the implementation in mpfr. At very high nn 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 nn, 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 nn 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 BnB_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:

|Bn|=ζ(n)2n!(2π)n<(1+22n)2enlnnn(2π)n \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 BnB_n is

log2(1+22n)+1+(nlnnn)log2enlog22π \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+(nlnnn)log2enlog22π 42 + (n \ln n - n) \log_2 \mathrm{e} - n \log_2 2 \pi

Less bits in zeta

In the last half of the nn the ζ(2n)1\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 122n\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

Khinchin.cpp

Benchmark

Line Method
Red line mpfr_zeta_ui
Blue line mpfr_zeta_ui with our von Staudt-Clausen
Purple line Darvid Harvey’s bernoulli
Green line Series summation without recycling
Black line Series summation with power table
Orange line Power table + von Staudt-Clausen
Benchmark for ζ(2n) with 34 kb accuracy
Benchmark for ζ(2n) with 34 kb accuracy
Bench­mark for ζ(2n) with 340 kb ac­cu­ra­cy
Bench­mark for ζ(2n) with 340 kb ac­cu­ra­cy