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