Calculating Khinchin’s constant

Remco Bloemen

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

Today I and two friends decided that Xavier Gourdon’s 1997 result of 110 000 digits of Khinchin’s constant, K0K_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 K0K_0 on the usual sources, wikipedia and mathworld.

The defining series used to calculate K0K_0 is:

lnK0=1ln2[1,]nζ(2n)1n([1,2n1]k(1)k+1k) \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 S1S_1 and the second summation S2(n)S_2(n) such that:

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

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

Trivial implementation:

Now, one could calculate K0K_0 using something like:

Optimisation

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

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

S2(n)S2(n1)=([1,2n1]k[1,2(n1)1]k)(1)k+1k=[2n2,2n1]k(1)k+1k=12(1n)(2n1) \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:

We eliminated a loop!

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

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 ζ(2n)\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);
                mpfr_add(S2, S2, temp, 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);
                mpfr_add(S1, S1, temp, 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;
}