Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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|improve this question
3  
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

1 Answer 1

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 *

@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
share|improve this answer

Your Answer

 
discard

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.