I am looking for a simple, inexpensive and very accurate approximation of the Lambert function ($W_0$ branch) ($-1/e < x < 0$).


6 Answers 6


Here is a post that I made on sci.math a while ago, regarding a method I have used.

Analysis of $we^w$

For $w\gt0$, $we^w$ increases monotonically from $0$ to $\infty$. When $w\lt0$, $we^w$ is negative.

Thus, for $x\gt0$, $\mathrm{W}(x)$ is positive and well-defined and increases monotonically.

For $w\lt0$, $we^w$ reaches a minimum of $-1/e$ at $w=-1$. On $(-1,0)$, $w e^w$ increases monotonically from $-1/e$ to $0$. On $(-\infty,-1)$, $w e^w$ decreases monotonically from $0$ to $-1/e$. Thus, on $(-1/e,0)$, $\mathrm{W}(x)$ can have one of two values, one in $(-1,0)$ and another in $(-\infty,-1)$. The value in $(-1,0)$ is called the principal value.

The iteration

Using Newton's method to solve $we^w$ yields the following iterative step for finding $\mathrm{W}(x)$: $$ w_{\text{new}}=\frac{xe^{-w}+w^2}{w+1} $$

Initial values of $w$

For the principal value, when $-1/e\le x\lt0$, and when $0\le x\le10$, use $w=0$. When $x\gt10$, use $w=\log(x)-\log(\log(x))$.

For the non-principal value, when $x$ is in $[-1/e,-.1]$, use $w=-2$; and if $x$ is in $(-.1,0)$, use $w=\log(-x)-\log(-\log(-x))$.

This says that for the branch you want, use the iteration with an initial value of $w=0$.


Corless, et al. 1996 is a good reference for implementing the Lambert $W$. One key to being efficient is a good initial guess for the iterative refinement. This is where series approximations are useful. For $(-1/\text{e}, 0)$ it's helpful to use two different expansions: about $0$ and about the branch point $-1/\text{e}$. Otherwise you can potentially waste a lot of time iterating.

Here's some unvectorized (for clarity) Matlab code for $W_0$. I've indicated the equation numbers from Corless, et al. on which various lines are based.

function w = fastW0(x)
% Lambert W function, upper branch, W_0
% Only valid for real-valued scalar x in [-1/e, 0]

if x > -exp(-sqrt(2))
    % Series about 0, Eq. 3.1
    w = ((((125*x-64)*x+36)*x/24-1)*x+1)*x;
    % Series about branch point, -1/e, Eq. 4.22
    p = sqrt(2*(exp(1)*x+1));
    w = p*(1-(1/3)*p*(1-(11/24)*p*(1-(86/165)*p*(1-(769/1376)*p))))-1;

tol = eps(class(x));
maxiter = 4;
for i = 1:maxiter
    ew = exp(w);
    we = w*ew;
    r = we-x; % Residual, Eq. 5.2
    if abs(r) < tol

    % Halley's method, for Lambert W from Corless, et al. 1996, Eq. 5.9
    w = w-r/(we+ew-0.5*r*(w+2)/(w+1));

Here's a plot (for double precision, tolerance of machine epsilon) of number of the iterations until convergence if using no initial series guess, the series about $0$, the series about the branch point, or both series:

iterations vs. x

Of course the evaluation of the series approximation has a cost relative to performing an iteration. One can adjust how many terms are used in order to tune this.


An iterative procedure, given in wikipedia, q&d-translated to Pari/GP, which suits my needs well:

 LW(x, prec=1E-80, maxiters=200) = local(w, we, w1e); 
           w = w-(we-x)/(w1e-(w+2)*(we-x)/(2*w+2)) ;
    print("W doesn't converge fast enough for ",x,"( we=",we);


   x = -0.99999/exp(1)
 %382 = -0.367875762377

 %383 = -0.995534517079

 %384 = 
[2.49762470622 E-163]
  • $\begingroup$ The maximum number of iterations is : 1 ! $\endgroup$ Jun 14, 2013 at 11:00
  • $\begingroup$ @CLaude: ? I just tried x=-0.2 and got the result... $\endgroup$ Jun 14, 2013 at 11:17
  • $\begingroup$ Hi Gottfried ! Try at x=-0.999/e $\endgroup$ Jun 14, 2013 at 11:19
  • $\begingroup$ @Claude: perhaps you use too low precision? I'm using 200 digits per default, so I'm possibly overlooking that problem? I've appended one run for checking in my answer. $\endgroup$ Jun 14, 2013 at 11:32
  • $\begingroup$ The question arises: is this an approximation or a calculation of $LambertW$? Maple produces $evalf(LambertW(- \frac 1 {\exp(1)}),20)= 1.$. $\endgroup$
    – user64494
    Jun 14, 2013 at 11:49

Maple produces $$\operatorname{Lambert W}(x)=x-{x}^{2}+{\frac {3}{2}}{x}^{3}-{\frac {8}{3}}{x}^{4}+{\frac {125}{24 }}{x}^{5}+O \left( {x}^{6} \right), x\to 0.$$ Let us compare $\operatorname{Lambert W}(-.2)=-.2591711018$ with $-0.2-(-0.2)^{2}+3/2(-0.2)^{3}-8/3(-0.2)^{4}+{\frac {125}{24}}(-0.2)^{5}= -.2579333334$.

  • $\begingroup$ Maple also produces its Pade approximation $$LambertW(x)= \left( {\frac {451}{340}}\,{x}^{3}+{\frac {228}{85}}\,{x}^{2}+x \right) \left( 1+{\frac {313}{85}}\,x+{\frac {1193}{340}}\,{x}^{2}+{ \frac {133}{204}}\,{x}^{3} \right) ^{-1} $$ which gives the value $LambertW(-0.2)=-.2591579777$. $\endgroup$
    – user64494
    Jun 14, 2013 at 9:25
  • $\begingroup$ Thanks for this info. The problem is that I need the approximation to be almost as good close to zero as close to -1/e. $\endgroup$ Jun 14, 2013 at 10:58
  • $\begingroup$ The Pade approximation $[3,3]$ is that you want:see a pdf file. $\endgroup$
    – user64494
    Jun 14, 2013 at 11:32

For the range $-\frac 1e \leq x \leq 0$, using $y=\sqrt{2(1+ex)}$ (as in this post) and building Padé approximants $P_n$ around $y=0$, what is obtained is $$P_1=\frac{-1+\frac{2 y}{3} } {1+\frac{y}{3} }$$ $$P_2=\frac{-1+\frac{14 y}{45}+\frac{301 y^2}{1080} } {1+\frac{31 y}{45}+\frac{83 y^2}{1080} }$$ $$P_3=\frac{-1-\frac{974 y}{22659}+\frac{21865 y^2}{51792}+\frac{89131 y^3}{1165320} } {1+\frac{23633 y}{22659}+\frac{104225 y^2}{362544}+\frac{3167 y^3}{196560} }$$ $$P_4=\frac{-1-\frac{11637254 y}{29330279}+\frac{463636649 y^2}{1055890044}+\frac{23930361857 y^3}{110868454620}+\frac{192684057311 y^4}{10643371643520} } {1+\frac{40967533 y}{29330279}+\frac{659231191 y^2}{1055890044}+\frac{1928737771 y^3}{20157900840}+\frac{34384971553 y^4}{10643371643520} }$$ The maximum errors occur for $x=0$; for the above formulae, they are $\Delta_1=-3.88\times 10^{-2}$, $\Delta_2=-1.23\times 10^{-3}$, $\Delta_3=-3.73\times 10^{-5}$ and $\Delta_4=-1.12\times 10^{-6}$.

Concerning $$S_i=\int_{-\frac 1e}^0 \left(W(x)-P_i\right)^2\,dx$$ they are $S_1=3.93\times 10^{-4}$, $S_2=2.57\times 10^{-7}$, $S_3=1.76\times 10^{-10}$ and $S_4=1.26\times 10^{-13}$.


For a "sanity" check, I also considered all the possible $[k,6-k]$ Padé approximants and obtained the following results $$\left( \begin{array}{ccc} k & \Delta_k & S_k \\ 0 & -0.108806 & 1.84 \times 10^{-3} \\ 1 & -0.007522 & 6.16 \times 10^{-6} \\ 2 & -0.000171 & 3.49 \times 10^{-9}\\ \color{red} { 3} & \color{red} {-0.000037} & \color{red} {1.76 \times 10^{-10}}\\ 4 & -0.000082 & 8.17 \times 10^{-10}\\ 5 & -0.000966 & 1.05 \times 10^{-7}\\ 6 & -0.094985 & 9.04 \times 10^{-4} \end{array} \right)$$


To compensate for massive cancellations when $x≈-1/e$

Let $r = 1/e - ε$
Let $h = (x+r)+ε$

$z_0 = \sqrt{2\,r\,h} \;≈\; e^{W(x)} - r$

Note: negative square root can get $W_{-1}$. Thus, I write $W$, instead of $W_0$
We multiply by a factor for improvement: click here for explanation

$\displaystyle z = z_0 \;\sqrt{1+\frac{z_0}{3r}}$

$\displaystyle → W(x) = \frac{x}{e^{W(x)}} \;≈\; \frac{x}{r+z}$

For $-1/e < x < 0$, formula is exact for both end points.
Max error $=9·10^{-5}$, at $x=-0.17$

We can apply Newton's method for much improvment.
Max error $=3·10^{-8}$, at $x=-0.12$

$\displaystyle → W(x) \;≈\; \frac{x\,\ln(1+e\,z)}{h+z}$


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .