3
$\begingroup$

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.)

  • 3
    The short almost-answer: Look at the [Euler–Maclaurin formula](http://en.wikipedia.org/wiki/Euler%E2%80%93Maclaurin_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. – 2012-11-09

1 Answers 1

1

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