Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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

share|cite|improve this question
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
up vote 1 down vote accepted

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 *

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]

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):
    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
share|cite|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.