1# principal branch of the Lambert W function and its minus one branch
2
3SetHelp ("LambertW", "functions", "Principal branch of the Lambert W function for real values greater than or equal to -1/e");
4function LambertW(x) = (
5  if not IsReal(x) or x < -1/e then (error ("LambertW: only defined for real x >= -1/e");bailout);
6  if (x < 0.367) then
7    y := x-x^2+(3/2.0)*x^3
8  else if (x < 24) then
9    (lnxm1 := ln(x)-1; y := 1+lnxm1/2.0+(lnxm1^2)/16.0)
10  else
11    (l1:=ln(x);l2:=ln(l1);y:=l1-l2+l2/l1);
12  y := y-(y*e^y-x)/(e^y+y*e^y);
13  y := y-(y*e^y-x)/(e^y+y*e^y);
14  y := y-(y*e^y-x)/(e^y+y*e^y);
15  y := y-(y*e^y-x)/(e^y+y*e^y);
16  for n=1 to 100 do (
17    ylast := y;
18    y := y-(y*e^y-x)/(e^y+y*e^y);
19    if y == ylast then return y;
20    ylast2 := y;
21    y := y-(y*e^y-x)/(e^y+y*e^y);
22    if (y == ylast2 or y == ylast) then return y
23  );
24  y
25)
26
27
28SetHelp ("LambertWm1", "functions", "The minus-one branch of the Lambert W function for real values between -1/e and 0");
29function LambertWm1(x) = (
30  if not IsReal(x) or x < -1/e or x >= 0.0 then (error ("LambertWm1: only defined for real x with -1/e <= x < 0");bailout);
31  if x < -0.116 then y := -8.22*x-4.23 else (l1:=ln(-x);l2:=ln(-l1);y:=l1-l2+l2/l1);
32  for n=1 to 10 do (
33    ylast := y;
34    y := y-(y*e^y-x)/(e^y+y*e^y);
35    if (y == ylast) then return y
36  );
37  y := y-(y*e^y-x)/(e^y+y*e^y);
38  y := y-(y*e^y-x)/(e^y+y*e^y);
39  y := y-(y*e^y-x)/(e^y+y*e^y);
40  y := y-(y*e^y-x)/(e^y+y*e^y);
41  for n=1 to 100 do (
42    ylast := y;
43    y := y-(y*e^y-x)/(e^y+y*e^y);
44    if y == ylast then return y;
45    ylast2 := y;
46    y := y-(y*e^y-x)/(e^y+y*e^y);
47    if (y == ylast2 or y == ylast) then return y
48  );
49  y
50)
51
52