Efficient ways to calculate Lambert $W_{-1}(z)$ with an arbitrary precision.

Question

I proved and tested that for any number $z\in [-\frac{1}{e}, 0]$ and a number $\epsilon$ the following algorithm returns a number $x$ such that $|W_{-1}(z)-x| < \epsilon$ in time $O\left(\left\lceil\log_{4|W_{-1}(z)|}\left(\frac{1}{\epsilon}\right)\right\rceil\right)$ i.e. it would take $O(1)$ for values of $z$ close to $0$ and it would take $O\left(\log\left(\frac{1}{\epsilon}\right)\right)$ for values of $z$ close to $-\frac{1}{e}$.

W_{-1}(z, epsilon)
{
    upper = log(z / log(-z));
    lower = upper + log(1 - (1.0 / e));
    x     = log(2.0 * z / (lower + upper));
    while (upper - lower > 2.0 * epsilon)
    {
        if (2.0 * x < lower + upper)
        {
            upper = x;
            lower = log(z / lower);
        }
        else
        {
            upper = log(z / upper);
            lower = x;
        }
        x = log(2.0 * z / (lower + upper));
    }

    return x;
}

And I wanted to know whether or not my algorithm was useful or if maybe I should publish it or something. So I tried to search what other algorithms there are to calculate the $-1$ branch of the lambert W function with an $\epsilon$ precision but I didn't find anything (maybe I'm just not very good at searching). What efficient ways there are to calculate the $W_{-1}$ with an arbitrary precision?

Answer 1

This is just fixed-point iteration, which only converges linearly. The bounds are also unnecessary because the iterations from either side will converge in an expected manner, except for $z$ near $-1/e$ and $-0$. It would be better to use something such as Newton's method or Halley's method, which converge much faster.

Newton's method gives the iterations:

$$w_{n+1}=w_n-\frac{w_n-ze^{-w_n}}{w_n+1}\tag0$$

and with the initial estimate

$$w_0=\ln\left(\frac z{\alpha\ln(-z)}\right),\quad\alpha=\frac2{2-\ln(2)}\tag1$$

you can guarantee iterations monotonically converge to the Lambert W function.

To get faster convergence near $-1/e$ you may wish to use instead the initial estimate

$$w_0=-1-\sqrt{2(1+ez)}\tag2$$

which provides an upper bound. Taking one step further with fixed-point iteration, one may instead use

$$w_0=\ln\left(\frac{-z}{1+\sqrt{2(1+ez)}}\right)\tag3$$

which is accurate both as $z\to-1/e$ and $z\to-0$. To obtain another asymptotic term as $z\to-0$, another iteration of fixed-point iteration gives

$$w_0=\ln\left(\frac z{\ln\left(\frac{-z}{1+\sqrt{2(1+ez)}}\right)}\right)\tag4$$

If your concern is having both upper and lower bounds to guarantee an error bound, you may wish to try using the larger step-size

$$w_{n+1}=\min\left\{-1,w_n-\beta\frac{w_n-ze^{-w_n}}{w_n+1}\right\},\quad\beta>1\tag5$$

every so often to see if the solution becomes bounded. Initial bounds may be obtained by noting that $(4)$ is an upper bound while applying $(5)$ to $(4)$ with $\beta=-w_n$ gives a lower bound.

I would suggest, since one can be sure the above converges, to use $(5)$ only once $|w_{n+1}-w_n|$ is sufficiently below the desired precision. The singular initial estimate $(4)$ should also suffice to give rapid convergence.