7
$\begingroup$

What is the best way to numerically calculate the inverse of the logarithmic integral, defined by

$$ \operatorname{li}(x)=\int_{0}^{x}\dfrac{1}{\log(t)}\operatorname{d}t $$

eg $\operatorname{li}(x)=100,\ x=?$

$\endgroup$
  • $\begingroup$ Ask WA. $\endgroup$ – lhf Jul 1 '14 at 13:17
  • $\begingroup$ May I suggest to write $\operatorname{li}(x)=\int_{0}^{x}\dfrac{1}{\log(t)}\operatorname{d}t$ ? Cheers. $\endgroup$ – Claude Leibovici Jul 3 '14 at 16:07
  • $\begingroup$ @Claude Leibovici, yes - thanks $\endgroup$ – martin Jul 4 '14 at 15:03
13
+50
$\begingroup$

The answer is $x=\operatorname{li}^{-1}(100) \approx 488.871909852807532.\;$ My implementation computes $\operatorname{li}^{-1}(x)$ as the zero of the function $f(z)=\operatorname{li}(z)-x$ using Halley iterations $$z_{n+1}=z_n-\frac{f(z_n)}{f'(z_n)}\left[1 - \frac{f(z_n)}{f'(z_n)}\frac{f''(z_n)}{2f'(z_n)}\right]^{-1}\\ =z_n-\Delta_n \ln(z_n) \left[1 + \frac{\Delta_n}{2z_n}\right]^{-1}$$

with $\Delta_n = f(z_n) = \operatorname{li}(z_n)-x.$ At most 3 iterations are needed for all $x$ (using 80 bit extended floats) with the following starting values $z_0$ \begin{equation*} z_0 = z_0(x) = \begin{cases} x \ln x &x > 3.5\\ 1 + x &0.75 < x \le 3.5\\ 1.45137 + 0.37251\cdot x &-0.5 < x \le 0.75\\ 1 + e^{x-\gamma} &-43.8 < x \le -0.5 \end{cases} \end{equation*} (For $x\le -43.8$ the result is $1$ accurate to extended precision.)

Edit: The 1. and 4. starting values are from asymptotic expansions $$x\rightarrow \infty: \operatorname{li}(x) \sim \frac{x}{\ln x} \Longrightarrow \operatorname{li}(x \ln x) \sim \frac{x\ln x }{\ln \ln x+ \ln x} \sim x$$ $$x\rightarrow 0^{+}: \operatorname{li}(1+x) \sim \ln x + \gamma$$ the third from the Taylor expansion for $\operatorname{li}$ around its zero at $x_0 = 1.451369234883381\dots$ and the second is plain heuristic. For $x=100$ the Halley iterations are as follows (first is the starting value):

 460.517018598809137
 488.874968191031427
 488.871909852807528
 488.871909852807532
$\endgroup$
  • $\begingroup$ Care to explain how you get these results? $\endgroup$ – lhf Jul 1 '14 at 12:43
  • $\begingroup$ @lhf: Do you mean the starting values? $\endgroup$ – gammatester Jul 1 '14 at 12:45
  • $\begingroup$ I get $488.87190985279381\ldots$ using bisection in the interval $[2,1000]$ and adaptive Simpson's rule to compute $Li(x)+li(2)$ for $x\ge 2$, taking $li(2)=1.045163780117492784844588$. $\endgroup$ – lhf Jul 1 '14 at 13:01
  • $\begingroup$ @lhf: Added some information about starting values and iteration progress. $\endgroup$ – gammatester Jul 1 '14 at 14:02
  • $\begingroup$ Very nice, thanks. $\endgroup$ – lhf Jul 1 '14 at 14:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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