11
\$\begingroup\$

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.

enter image description here

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.

Improve this question
\$\endgroup\$
8
  • 1
    \$\begingroup\$ 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. \$\endgroup\$ Commented Nov 23, 2017 at 5:58
  • 1
    \$\begingroup\$ 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. \$\endgroup\$ Commented Nov 23, 2017 at 8:04
  • 1
    \$\begingroup\$ If one input should give two values, this is not a function \$\endgroup\$ Commented Nov 23, 2017 at 10:43
  • \$\begingroup\$ @Mr.Xcoder Sure, though the non-real output should be reasonable. \$\endgroup\$ Commented Nov 23, 2017 at 12:50
  • 2
    \$\begingroup\$ 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? \$\endgroup\$ Commented Nov 23, 2017 at 16:56

9 Answers 9

Reset to default
3
\$\begingroup\$

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
Improve this answer
\$\endgroup\$
2
\$\begingroup\$

05AB1E, 98 86 75 74 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н, { { } { } { { } { } {
Improve this answer
\$\endgroup\$
1
  • \$\begingroup\$ 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. \$\endgroup\$ Commented Apr 20, 2020 at 12:10
1
\$\begingroup\$

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.

Improve this answer
\$\endgroup\$
4
  • 3
    \$\begingroup\$ Why non-competing? \$\endgroup\$ Commented Nov 23, 2017 at 3:38
  • 4
    \$\begingroup\$ 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. \$\endgroup\$ Commented Nov 23, 2017 at 4:51
  • \$\begingroup\$ 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). \$\endgroup\$ Commented Nov 23, 2017 at 9:31
  • 2
    \$\begingroup\$ @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 \$\endgroup\$ Commented Nov 23, 2017 at 12:49
1
\$\begingroup\$

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!

Improve this answer
\$\endgroup\$
1
\$\begingroup\$

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.

Improve this answer
\$\endgroup\$
1
\$\begingroup\$

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 ;
Improve this answer
\$\endgroup\$
0
\$\begingroup\$

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!

Improve this answer
\$\endgroup\$
0
\$\begingroup\$

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

Improve this answer
\$\endgroup\$
0
\$\begingroup\$

SageMath, 39 bytes

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

Try it on SageCell!

Improve this answer
\$\endgroup\$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.