# Calculating Khinchin's constant

Today I and two friends decided that Xavier Gourdon’s 1997 result of 110 000 digits of Khinchin’s constant, $K_0$ could be improved. Khinchin’s constant is the result of an interesting theorem on the properties of continued fraction expansion, you can read more about $K_0$ on the usual sources, wikipedia and mathworld.

The defining series used to calculate $K_0$ is:

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

Let’s call the first summation $S_1$ and the second summation $S_2(n)$ such that:

$$S_2(n) = \sum_{[1, 2n-1]}^{k} \frac{(-1)^{k+1}}{k}$$

$$S_1(N) = \sum_{[1, N]}^{n} \frac{\zeta(2n) -1}{n} S_2(n)$$

## Trivial implementation:

Now, one could calculate $K_0$ using something like:

double K0, S1, S2;
S1 = 0.0;
for (int n = 1; n <= N; n++) {
S2 = 0.0;
for(int k = 1; k <= 2n-1; k++) {
S2 += (-1)^(k+1)/k;
}
S1 += S2 · (ζ(2n)-1)/n;
}
K0 = exp(S1 / ln(2));


## Optimisation

Right now the implementation recalculates $S_2(n)$ entirely for every n. It would be an improvement to re-use the previous value, $S_2(n-1)$ to calculate the new one, $S_2(n)$.

To implement this we calculate the difference between two successive series $S_2(n)$:

\begin{aligned} S_2(n) - S_2(n-1) &= \left( \sum_{[1, 2n-1]}^{k} - \sum_{[1, 2(n-1)-1]}^{k}\right) \frac{(-1)^{k+1}}{k} \\ &= \sum_{[2n-2, 2n-1]}^{k} \frac{(-1)^{k+1}}{k} \\ &= \frac{1}{2(1-n)(2n-1)} \end{aligned}

Which results in the following algorithm:

double K0, S1, S2;
S1 = ζ(2) - 1;
S2 = 0.0;
for (int n=2; n <= N; n++) {
S2 += 1/(2(1-n)(2n-1));
S1 += S2 · (ζ(2n)-1)/n;
}
K0 = exp(S1 / ln(2));


We eliminated a loop!

But how big should N be to reach convergence? One solution is to loop until $S_1$ no longer changes:

double K0, S1, S2, error, S1_old;
S1 = ζ(2) - 1;
S2 = 0.0;
unsigned int n = 2;
do
{
S1_old = S1;
S2 += 1/(2(1-n)(2n-1));
S1 += S2 · (ζ(2n)-1)/n;
difference = S1 - S1_old;
n++;
} while(difference > 0);
K0 = exp(S1 / ln(2));


## Implementation and analysis

My first thought was to use C++ and GMP, but GMP lacked support for calculating zeta functions. The documentation points to MPFR for floating point extensions. The MPFR library uses GMP for calculations and contains a ton of useful functions besides zeta.

After compiling (with the latest GCC of course) the program produced thousand digits in less than a second and two thousand digits in 7 seconds, we will never get to a hundred thousand at that rate.

So we analysed this program using callgrind and KCachegrind which revealed that 99% of the time is spent calculating the $\zeta(2n)$'s. This will be the target of optimisation in following posts.

## Complete source

// g++ -O3 -lmpfr

#include <stdio.h>
#include <stdarg.h>
#include <mpfr.h>

int one_loop_K0(mpfr_t K0)
{
mp_prec_t precision = mpfr_get_prec(K0);

mpfr_t zeta, S1, S1_old, error, S2, temp;
mpfr_init2(zeta, precision);
mpfr_init2(S1, precision);
mpfr_init2(S2, precision);
mpfr_init2(S1_old, precision);
mpfr_init2(temp, precision);
mpfr_init2(error, precision);

// n = 1
// S1 = ζ(2n) - 1
// S2 = 1
int n = 1;
mpfr_zeta_ui(S1, 2*n, GMP_RNDN);
mpfr_sub_ui(S1, S1, 1, GMP_RNDN);
mpfr_set_d(S2, 1.0, GMP_RNDN);

do {
n++;

// S1_old = S1
mpfr_set(S1_old, S1, GMP_RNDN);

// S2 += -1/(2(n-1)(2n-1))
mpfr_set_d(temp, -0.5, GMP_RNDN);
mpfr_div_ui(temp, temp, (n -1) * (2*n-1), GMP_RNDN);

// zeta = ζ(2n)
mpfr_zeta_ui(zeta, 2*n, GMP_RNDN);

// S1 += (ζ(2n) - 1) * S2 / n;
mpfr_sub_ui(temp, zeta, 1, GMP_RNDN);
mpfr_div_ui(temp, temp, n, GMP_RNDN);
mpfr_mul(temp, temp, S2, GMP_RNDN);

// error = S1 - S1_old
mpfr_sub(error, S1, S1_old, GMP_RNDN);

if(n % 100 == 0) {
int bits = -mpfr_get_exp(error);
printf("Term %d using series summation has revealed %d bits.n", n, bits);
}
} while (mpfr_cmp_ui(error, 0));

// temp =  1 / ln 2
mpfr_const_log2(temp, GMP_RNDN);
mpfr_ui_div(temp, 1, temp,  GMP_RNDN);

// K = exp(temp * S1)
mpfr_mul(K0, temp, S1,  GMP_RNDN);
mpfr_exp(K0, K0, GMP_RNDN);

mpfr_clear(zeta);
mpfr_clear(S1);
mpfr_clear(S1_old);
mpfr_clear(S2);
mpfr_clear(temp);
mpfr_clear(error);

return n;
}

int main()
{
mpfr_t K0;
mpfr_prec_t precision = 2 * 3400;
mpfr_init2(K0, precision);
int n = one_loop_K0(K0);
mpfr_printf("N = %d t %RNfn", n, K0);
mpfr_clear(K0);
mpfr_free_cache();
return 0;
} Remco Bloemen
Math & Engineering
https://2π.com