Compute the Lambert W function

Question

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

y = W(x) if x = y*(e^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.

Answer 1

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.

Answer 2

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

Answer 3

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 = 9e⁹.

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 = 10⁸.

#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!

Answer 4

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!