Efficient calculation of Khinchin's constant

In the previous post we (me, Sander & Sander) set upon ourselves the challenge to beat Xavier Gour­don’s 1997 record of hundred and ten thousand digits of Khinchin’s constant. This quickly diverted into a challenge of calculating the Riemann Zeta function (again, see wikipedia and mathworld for more info).

Khinchin’s constant, K_0 and the Zeta function ζ are related by the following identity:

\ln K_0 = \frac{1}{\ln{2}} \sum_{[1,\infty]}^{n} \frac{\zeta(2n) -1}{n} \sum_{[1, 2n-1]}^{k} \frac{(-1)^{k+1}}{k}

Analysing bounds and limits

The inner series is an alternating harmonic series and has the following limit:

\lim_{n \to \infty} \sum_{[1, 2n-1]}^{k} \frac{(-1)^{k+1}}{k} = \ln{2}

Since we are taking two terms at a time this series it is strictly increasing and \ln{2} is also an upper bound.

An upper bound for ζ can be easily seen from the series definition, \zeta(2n) - 1 < \frac{2}{2^{2n}}. Combing these results gives an upper bound for the terms in the outer series:

\frac{\zeta(2n) -1}{n} \sum_{[1, 2n-1]}^{k} \frac{(-1)^{k+1}}{k} < \frac{\ln{4}}{n 4^n}

This expression also allows us to calculate an upper bound for the cutoff error, i.e. the error we introduce by not summing to n=\infty but stopping at n=N:

\sum_{[N, \infty]} \frac{\ln{4}}{n 4^n} < \ln{4} \sum_{[N, \infty]} \frac{1}{4^n} = \frac{\ln{16}}{3 \cdot 4^N}

If we now take the binary logarithm of this number we know how many fractional bits we get at least correct when summing up to N. In fact, since \ln{2} \ln K_0 \approx 1, we know the total number of bits we get correct:

-\log_2 \left( \frac{\ln{16}}{3 \cdot 4^N} \right) = 2 N + \log_2 3 - \log_2 \ln{16}

Or, put differently, to obtain B bits of K_0 we must sum up to (but not including):

N =\frac{1}{2}\left(B - \log_2 3 + \log_2 \ln{16}\right) < \left\lceil\frac{B}{2}\right\rceil + 1

The uncertainty of the sum is the sum of the uncertainties of the terms. So if we need an accuracy of 2^{-B} for the sum, we can divide this accuracy evenly over the terms, giving each an accuracy of \frac{1}{N} 2^{-B}. So the precision (in bits) of the terms is:

B_{term} = \left\lceil B + \log_2 N \right\rceil

The ζ can be less accurate, since it will be multiplied by \ln 2 < 1 and divided by n, this amounts to an required accuracy of:

B_{zeta}(n) = \left\lceil B_{term} - \log_2 n -\log_2 \ln{2} \right\rceil

The million digit goal

To calculate K_0 accurate up to a million digits we require the following:

The problem can now be restated as: How does one calculate ζ(2) … ζ(2N) fast and accurate up to B_{zeta}(n) bits?