File: lambert.gel

package info (click to toggle)
genius 1.0.24-2
  • links: PTS, VCS
  • area: main
  • in suites: buster
  • size: 27,952 kB
  • sloc: ansic: 105,597; xml: 67,672; sh: 4,537; makefile: 2,089; lex: 499; yacc: 298; perl: 54; python: 22
file content (51 lines) | stat: -rw-r--r-- 1,664 bytes parent folder | download | duplicates (5)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
# principal branch of the Lambert W function and its minus one branch

SetHelp ("LambertW", "functions", "Principal branch of the Lambert W function for real values greater than or equal to -1/e");
function LambertW(x) = (
  if not IsReal(x) or x < -1/e then (error ("LambertW: only defined for real x >= -1/e");bailout);
  if (x < 0.367) then
    y := x-x^2+(3/2.0)*x^3
  else if (x < 24) then 
    (lnxm1 := ln(x)-1; y := 1+lnxm1/2.0+(lnxm1^2)/16.0)
  else
    (l1:=ln(x);l2:=ln(l1);y:=l1-l2+l2/l1);
  y := y-(y*e^y-x)/(e^y+y*e^y);
  y := y-(y*e^y-x)/(e^y+y*e^y);
  y := y-(y*e^y-x)/(e^y+y*e^y);
  y := y-(y*e^y-x)/(e^y+y*e^y);
  for n=1 to 100 do (
    ylast := y;
    y := y-(y*e^y-x)/(e^y+y*e^y);
    if y == ylast then return y;
    ylast2 := y;
    y := y-(y*e^y-x)/(e^y+y*e^y);
    if (y == ylast2 or y == ylast) then return y 
  );
  y
)


SetHelp ("LambertWm1", "functions", "The minus-one branch of the Lambert W function for real values between -1/e and 0");
function LambertWm1(x) = (
  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);
  if x < -0.116 then y := -8.22*x-4.23 else (l1:=ln(-x);l2:=ln(-l1);y:=l1-l2+l2/l1);
  for n=1 to 10 do (
    ylast := y;
    y := y-(y*e^y-x)/(e^y+y*e^y);
    if (y == ylast) then return y 
  );
  y := y-(y*e^y-x)/(e^y+y*e^y);
  y := y-(y*e^y-x)/(e^y+y*e^y);
  y := y-(y*e^y-x)/(e^y+y*e^y);
  y := y-(y*e^y-x)/(e^y+y*e^y);
  for n=1 to 100 do (
    ylast := y;
    y := y-(y*e^y-x)/(e^y+y*e^y);
    if y == ylast then return y;
    ylast2 := y;
    y := y-(y*e^y-x)/(e^y+y*e^y);
    if (y == ylast2 or y == ylast) then return y 
  );
  y
)