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);
                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;
}