Numerical inverse of logarithmic integral

Question

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

Answer 1

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