Compute the Lambert W function

Your challenge is to compute the Lambert W function. $$\W(x)\$$ is defined to be the real value(s) $$\y\$$ such that

$$y = W(x) \text{ if } x = ye^y$$

where $$\e = 2.718281828...\$$ is Euler's number.

Sometimes, $$\y\$$ may not be real.

Examples

W(-1) = non-real
W(-0.1) = -0.11183, -3.57715
W(1) = 0.56714
W(2) = 0.85261


Here's a quick graph of what this function looks like.

Rules

Your goal is to take an input and output either nothing, 1 solution, or 2 solutions, out to 5 significant figs. You should expect float inputs within the reasonable range of -100..100.

This is code-golf, so shortest code wins.

• A note about Lambert W function: there is no solution for x < -1/e, two solutions for -1/e < x < 0, and one solution for x == -1/e or x >= 0. Commented Nov 23, 2017 at 5:58
• 5 places after the decimal, or 5 significant figures? Since W^-1(x) tends to negative infinity when x->0-, it is hard to get 5 places after the decimal. Commented Nov 23, 2017 at 8:04
• If one input should give two values, this is not a function Commented Nov 23, 2017 at 10:43
• @Mr.Xcoder Sure, though the non-real output should be reasonable. Commented Nov 23, 2017 at 12:50
• I think a nice way to handle multiple solutions is the output being a list of reals, containing 0, 1, or 2 values. Would this be allowed?
– xnor
Commented Nov 23, 2017 at 16:56

Octave, 36 34 bytes

@(x)(y=lambertw(-1:0,x))(~imag(y))


Try it online!

I solved this in 10 different ways without the builtin, but I couldn't get the solution for both the -1 branch and the 0-branch. I wasn't able to solve this without resorting to the builtin lambertw function. Even though it uses a builtin, it's far from straightforward.

lambertw will only give out one value by default, namely the W0(x), unless otherwise specified. The problem arises when we specifies that we want two values, both W0(x) and W-1(x). If this is the case, and x>=0, it will return one real and one complex result.

We therefore have to specify that we only want the results where the imaginary part of the result is zero.

Breakdown:

@(x)                              % Anonymous function that takes x as input
(y=lambertw(-1:0,x))          % The lambertw function, with both the -1 and 0 branches
(~imag(y))  % Index, where the elements with imag(y)=0 is true


05AB1E, 98867574 73 bytes

-11 bytes thanks to @KevinCruijssen

žr(zV"Ðžrsm*I-d"UY‹iõë∞.ΔX.V}®s‚ΔDÅAX.Vǝ}н,0‹i∞.Δ(X.V}(®‚ΔDÅAX.V≠ǝ}IYÊiн,


Try it online!

Explanation:

(I will let f(x)=x*e^x)

žr(zV"Ðžrsm*I-d"U #Let Y=-1/e,X=str(ÐYsm*I-d) to save characters in the future

Y‹iõ              #if input<-1/e: print empty string

ë∞.ΔX.V}          #else: find least int k>=1 such that f(k)>=0 (evaluates string X as though it were a function)

®s‚               #make list [-1,k]

ΔDÅAX.Vǝ}         #find root of f(x)-input in the interval [-1,k] using bisection method

#actually this outputs a list with 2 numbers which will be very close to each other, which represents the interval in which the root lies

#e.g. the list [-3.00001,-3] would represent the interval [-3.00001,-3]

н,                #print the first item of the list; stack is currently empty

0‹i               #if input<0:

∞.Δ(X.V}(         #find greatest int k<=-1 such that f(k)>=0

®‚                #make list [k,-1]

ΔDÅAX.V≠ǝ}        #find root of f(x)-input in the interval [k,-1] using bisection method

#again, this outputs a list rather than a single number

IYÊiн,            #if input!=-1/e, print that root (done so that "-1" is not printed twice)


Visualization of the program's indentation:

žr(zV"Ðžrsm*I-d"UY‹iõë∞.ΔX.V}®s‚ΔDÅAX.Vǝ}н,0‹i∞.Δ(X.V}(®‚ΔDÅAX.V≠ǝ}IYÊiн,
{    {   }   {       }    {  {    }   {        }   {

• I'm sure this can be golfed some more, but here is 75 bytes by replacing both 1( with ® (relevant 05AB1E tip), and putting ÐYsm*I-d inside a string variable at the start, so we can use it four times with .V. Commented Apr 20, 2020 at 12:10

Ruby, 131 bytes

->x{y,z,e,l=1,-1,Math::E,->t{Math::log t};100.times{z=x<0&&x>-1/e ?l[x/z]:p;y=x<-1/e ?"non-real":x>e ?l[x/y]:x*e**-y};p y;z&&(p z)}


Try it online!

Basically applies fixed-point iteration 100 times, as partially explained in Lambert W function aproximation.

• Why non-competing? Commented Nov 23, 2017 at 3:38
• Marking a submission non-competing does not make it valid. If your submission is invalid, please edit it to make it valid or delete it altogether. Relevant meta post. Commented Nov 23, 2017 at 4:51
• Actually the only non-competing reason I have seen is language postdate the challenge. In other cases I often see they are deleted or competing. / Well, apply the fixed-point iteration is a valid way to compute the function, but only if you can prove that it is good enough (to your own requirement). Commented Nov 23, 2017 at 9:31
• @JungHwanMin My submission is valid, I just marked it non-competing because, well, I didn't want to compete. I suppose that in itself is not a valid reason :P Commented Nov 23, 2017 at 12:49

C, 101 + 2 (-lm) bytes

Returns how many real solutions, and store them in an array.

This one is very slow, since it uses a brute-force approach to find roots. It handles up to x = 9e9.

f(a,b,p,t,c,x)float*b,a,x;{for(p=a<0,c=0,x=-50;x<9;x+=2e-6,p=t)(t=x*exp(x)>a)^p?b[c++]=x:0;return c;}


Try it online!

C, 148 146 + 2 (-lm) bytes

This function returns how many real solutions, and store them in the given pointers.

Use binary search to find a solution. It handles up to about x = 108.

#define h(x)99;for(l=-1;x=(l+r)/2,fabs(l-r)>1e-6;)x*exp(x)>a?r=x:(l=x);
f(a,b,c,l,r)float*b,*c,a,l,r;{r=h(*b)r=-h(*c)return-1/exp(1)>a?0:a<0?2:1;}


Try it online!

Casio BASIC, 192 Bytes

A port of My 05AB1E answer

?→X
""
If X≥-e^-1
Then 0→R
While Re^R<X
Isz R
WhileEnd
-1→L
While Re^R>X
.5(L+R→M
If Me^M<X
Then M→L
Else M→R
IfEnd
WhileEnd
If X<0
Then M◢
While Le^L<X
Dsz L
WhileEnd
While Le^L≠X
.5(L+R→M
If Me^M<X
Then M→R
Else M→L
IfEnd
WhileEnd
If M≠-1
Then M


Note that indents are actually not supported; I added them for legibility.

Forth (gforth), 144 bytes

: N 99 for fover fover fnegate fexp f* fover 2e f** f+ fswap 1e f+ f/ next f. ;
: f fdup 0e -1e fexp f- f> if 0e N fdup f0< if -2e N then then ;


Try it online!

Being stuck with the floating-point stack hurts a lot (with all the f-prefixes), but I guess it's pretty good score for a non-built-in after all.

Uses 100 iterations of Newton's method described here: given a previous estimate of $$\w=W(x)\$$, the next estimate $$\w'\$$ is calculated using the formula

$$w' = \frac{xe^{-w}+w^2}{w+1}$$

I simplified the starting points for the golf: starting at 0 for primary branch (if $$\x \ge -1/e\$$) and -2 (inflection point) for secondary branch (if $$\-1/e \le x < 0\$$) should suffice theoretically.

: N ( x w -- x )
99 for                        \ ( x w ) run Newton's method 100 times...
fover fover fnegate fexp f* \ ( x w x*e^-w )
fover 2e f** f+             \ ( x w x*e^-w+ww )
fswap 1e f+ f/              \ ( x (x*e^-w+ww)/(w+1) )
next
f.     \ Print the resulting w
;
: f ( x -- )
fdup 0e -1e fexp f-  \ ( x x -1/e )
f> if                \ ( x ) if less than -1/e, no solution
0e N               \ ( x ) print primary solution
fdup f0< if        \ ( x ) if less than 0, two solutions
-2e N f.         \ ( x ) print secondary solution
then
then
;


Wolfram Language (Mathematica), 28 bytes

x/.NSolve[#==E^x x,x,Reals]&


Try it online!

Returns x for no-solution.

This function should give all solutions because NSolve algorithm uses inverse functions, and there is a built-in inverse function of x = y*(e^y) called ProductLog.

Explanation

#==E^x x


Equation <input> = e^x * x

NSolve[ ... ,x,Reals]


Find solutions in the domain of real numbers.

x/.


Replace x with the solutions.

Alternative version, 35 bytes

N@ProductLog[{-1,0},#]~Cases~_Real&


Try it online!

Python 3, 258 bytes

e,f,X,R,L=2.7182818,lambda x:x*e**x,float(input()),0,-1
if X>=-1/e:
while f(R)<X:R+=1
while f(R)>X:
M=(L+R)/2
if f(M)<X:L=M
else:R=M
print(M)
if X<0:
while f(L)<X:L-=1
while f(L)!=X:
M=(L+R)/2
if f(M)<X:R=M
else:L=M
if M!=-1:print(M)


Try it online!

A port of my 05AB1E answer

SageMath, 39 bytes

lambda z:[lambert_w(l,z)for l in[0,-1]]


Try it on SageCell!