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