# Numerical inverse of logarithmic integral

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=?$

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

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

• Care to explain how you get these results? – lhf Jul 1 '14 at 12:43
• @lhf: Do you mean the starting values? – gammatester Jul 1 '14 at 12:45
• 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$. – lhf Jul 1 '14 at 13:01
• @lhf: Added some information about starting values and iteration progress. – gammatester Jul 1 '14 at 14:02
• Very nice, thanks. – lhf Jul 1 '14 at 14:22