# Computing a slowly-converging limit

Let $$f(x)=-\log\log x+\sum_{2\le n\le x}\frac{1}{n\log n}.$$

How can I efficiently compute $$f(\infty)=:\lim_{x\to\infty}f(x)?$$

Brute force suffices to find 0.7946786454... but I would like several hundred digits.

It seems that I should be able to use numerical integration, since $$\frac{1}{n\log n}-\log\log n+\log\log(n-1)$$ is smooth (and appears to be monotonic). (In fact, it even has a closed-form integral in li.) Alternately, various forms of series acceleration may apply.

(N.B. I have no real experience with numerical analysis outside an undergraduate class a few years back.)

-
The short almost-answer: Look at the Euler–Maclaurin formula and apply it to $\int 1/(x\log x)\,dx=\log\log x$. Use it to get good approximations to tails of the series. – Harald Hanche-Olsen Nov 9 '12 at 17:08

This can be evaluated numerically using the Euler-Maclaurin formula. Let $t(n) = (x\log x)^{-1}$ be the $n$-th term, pick a threshold $a$ and write the sum as $$\sum_{2\leq k < a}t(n) - \log\log a+\frac12t(a) + \sum_{k\geq 1}R_k,$$ $$R_k = \frac{B_{2k}}{(2k)!}\partial_n^{2k-1}\big|_{n=a}t(n).$$ The sum of $R_k$ is asymptotic in nature and diverges as $k\to\infty$, so it is necessary to sum $R_k$ only until the point $k$ when $|R_k|$ starts increasing.

Here is the value to what should be about $200$ digits: $$\begin{array}{rrrrr} 0.7946786454 & 5289940220 & 3897962065 & 1495140649 & 9959088280 \\ 4968901512 & 0950148178 & 5896068756 & 6696614777 & 6273344714 \\ 3933647273 & 2836173295 & 5585193860 & 6587583982 & 7169312934 \\ 2858106252 & 2997064914 & 8599822651 & 4777778609 & 7884785692 \end{array}$$ This was produced with the following code:

from mpmath import *

@memoize
@monitor
def fd(n, x, a=1, b=1):
# n-th derivative of 1/(x**a*log(x)**b)
if n < 0: raise ValueError('fd must have n >= 0')

dp = [[1/(x**(i+a)*log(x)**(j+b)) for j in range(n+1)] for i in range(n+1)]
for k in range(1,n+1):
for i in range(n+1-k):
for j in range(n+1-k):
dp[i][j] = -(b+j)*dp[i+1][j+1]-(a+i)*dp[i+1][j]
return dp[0][0]

@monitor
def qq(a = None):
if a is None: a = int(2 + ceil(mp.dps / 2.7))
ax = mpf(a)
def term(x):
x = mpf(x)
return 1/(x*log(x))
s = sum(term(x) for x in range(2, a)) - log(log(ax)) + term(a)/2

b, bPrev, bCount = None, None, 0
for k in range(1, 1000):
b = bernoulli(2*k)/gamma(2*k+1) * fd(2*k-1, ax)
if (bPrev is not None and abs(b) > abs(bPrev)) or abs(b) < mpf(10)**(-mp.dps):
break
bPrev = b
s -= b
bCount += 1
print('Precision %.3f (%.3f) from %d terms at a=%d' % (log10(abs(bPrev)), log10(abs(b)), bCount, a))

return s

-