Calculating Lambert W Function

Question

I'm trying to evaluate Lambert W Function , I used the formula $$ W(z)e^{W(z)} = z \Rightarrow W(z) = \frac{z}{W(z)} $$ $$ W(z) \approx ln(z)-ln(ln(z)-ln(...)) $$ But the result is very bad If I used ln(z)-ln(ln(z)) when it used in solving many cases like : $$ x^x = 100 \Rightarrow x = 3.015 $$ which is very far from 100 (27.68) . It works good in big values as shown in this graph. I've searched and found this equation : $$ L_1 - L_2 + \sum_{i=0}^\infty \sum_{j=0}^\infty (-1)^i \begin{bmatrix} i+j \\ i+1 \end{bmatrix}\frac{L_2\strut^{j}}{L_1\strut^{(i+j)}(i)!} \\ $$ $$ \text{Where} \hspace{10pt}L_1:ln(x), \hspace{10pt} L_2:ln(ln(x)) $$ Tried it but worse than the previous one and wolfram can't solve more than 15 terms but it can solve W(x) with prefect accuracy although i found the last formula on wolfram reference, so my question is " Is there a better way or formula to calculate Lambert W Function with high accuracy ? " another small question : I've used Taylor Series Expansion for Lambert W Function $$ W(z) = \sum_{k=0}^\infty (-k)^{(k-1)}\frac{z^k}{k!} $$ but the radius of convergence is 1/e ,so it's not useful for real computations and wolfram formula works for x > e is there a way to get the values in between 1/e and e

Answer 1

$\require{begingroup} \begingroup$ $\def\e{\mathrm{e}}\def\W{\operatorname{W}}\def\Wp{\operatorname{W_0}}\def\Wm{\operatorname{W_{-1}}}$

Is there any need to reinvent the wheel? Just use Halley's method, which provides successive approximations to $w = \W(z)$ (so $z = w\exp(w)$) being

\begin{align} w_{j+1} &= w_j-{\frac {w_j \exp(w_j)-z} {\exp(w_j)(w_j+1)- \displaystyle\frac{(w_j+2)(w_j\exp(w_j)-z)}{2w_j+2}}} \end{align}

as an established very efficient way to calculate $\W(x)$.

You can even check out open-source code, like for example, specfunc/lambert.c, which is part of GSL - GNU Scientific Library, for the details of realization.

$\endgroup$