Accurate computation of Gauss-Kuzmin entropy

Question

The Gauss-Kuzmin distribution gives the probability of an integer appearing as a partial denominator in the continued fraction of a real number $x$ as $$ P(a_k = k) = -\log_2\left(1 - \frac{1}{(k+1)^2} \right) $$

The entropy (in bits) of this distribution is \begin{align*} H = \sum_{k=1}^{\infty} \log_2\left(1 - \frac{1}{(k+1)^2} \right)\log_{2}\left(\log_2\left(1 - \frac{1}{(k+1)^2} \right) \right) \end{align*} The sequence is logarithmically convergent with only positive terms. I used 1000000 terms at 2000 digits of precision using MPFR, and was only able to recover 4 decimal digits. Others have indicated that this constant is very difficult to compute.

In the spirit of the SIAM 100 digit challenge, are there sequence acceleration tools that could make the sequence converge faster?

(Hope you have fun with this one!)

Answer 1

It's fairly easy to evaluate, to do this expand the logs in Taylor series in $x=(k+1)^{-2}$: $$ \log_2(1-x) = \frac{-1}{\log 2}\sum_{m\geq1}\frac{x^{m}}{m}$$ $$ \log_2(-\log_2(1-x)) = \frac{\log x}{\log 2} - \frac{\log\log 2}{\log 2} + \sum_{n\geq 1}a_n x^n, $$ where $a_n$ are Taylor series coefficients of the l.h.s. after the log-singularity is subtracted. These can be easily calculated directly using numerical quadrature (code below).

Using the identities (here, as above, $x=(k+1)^{-2}$) $$ \sum_{k\geq 1}x^s = \zeta(2s)-1, $$ $$ \sum_{k\geq 1}x^s\log x = 2\zeta'(2s), $$ we can rewrite the goal sum as a sum of three terms: $$ \frac{-1}{(\log 2)^2}\sum_{k,m\geq 1}\frac{x^m\log x}{m} = \frac{-1}{(\log 2)^2} \sum_{m\geq 1}\frac{2\zeta'(2s)}{m} = 4.067782509260337209451548799089911685899767925643908, $$ the second $$ \frac{\log\log 2}{(\log 2)^2}\sum_{k,m\geq 1}\frac{x^m}{m} = \frac{\log\log 2}{\log 2} = -0.52876637294489761424749777977881481518723706368332, $$ and a triple sum with separately calculated coefficients: $$ \frac{-1}{\log 2}\sum_{k,m,n\geq1} \frac{a_n}{m}x^{m+n} = \frac{-1}{\log 2}\sum_{m,n\geq1}\frac{a_n}{m}(\zeta(2m+2n)-1) = -0.1064886215397004957702703848332558106111776391548. $$ Because $\zeta(a)-1 = O(2^{-a})$, these sums converge extremely quickly, and can be evaluated directly without numerical extrapolation techniques, yielding the result $$ 3.4325275147757390994337806344778410601013532228057895 $$ It looks like it matches the number your link gives ($3.432527514776$) to at least the digits given there. I used fifty digits, and the most time-consuming portion of this is calculating the derivatives.

from mpmath import mp
from mpmath import *
import mpmath
mp.dps = 50

lg = lambda x: mp.log(x, b=2)

@mpmath.memoize
def A(n):
    if n == 0:
        return 0
    return (mp.diff(lambda x: lg(-lg(1-x))-(log(x)-log(log(2)))/log(2), 0, n, method="quad")/mp.gamma(n+1)).real

s1 = -log(2)**(-2)*nsum(lambda m: 2*zeta(2*m, derivative=1)/m, (1, mp.inf))
s2 = log(log(2))/log(2)**2 * nsum(lambda m: (zeta(2*m)-1)/m, (1, mp.inf))
s3 = -1/log(2) * nsum(monitor(lambda m, n: A(n)/m * (zeta(2*m+2*n)-1)), (1, mp.inf), (1, mp.inf), method="direct")

print(s1 + s2 + s3)