lambertw, all branches, Cas code

[Albert Chan]

Python W, version 4, translated for HP Prime Cas, with some differences.

1. HP Primc Cas use 48 bits, truncated float, instead of 53 bits, round-to-nearest.
   
   1/e ≈ 0x0.bc5ab1b1677ap-1 - 0x0.83951484e1f6p-50      // both terms 48 bits
         ≈ 0.3678794411714428 - 4.565179635228518e-16    // d = ceil(b*log10(2))+1
    
2. lyuka branch solve with e^W(a) = y = r+r*x, for x
   see https://www.hpmuseum.org/forum/thread-19...#pid183435
    
3. Cas code ignored edge cases phase angle, same as WA for Wk(0), Wk(±∞)
   
   Wk(±∞) = ∞
   Wk(0) = k ? -∞ : 0
   
   This is because Cas arg(±inf) = undef
   undef is similar to nan, which polluted everywhere: inf + undef*i → undef
   
   To get phase angle for infinite cases, customized arg() is needed, for little gain.
   

Code:

#cas
lambertw(a) :=
BEGIN
LOCAL h,s,x,y,z,cln,A,T,bad,(k:=0);
IF dim(a)!=1 THEN a,k := a END
a := float(a);
h := a + 0.3678794411714428 - 4.565179635228518e-16;
s := (im(a)<0)*2 - 1;
IF abs(h)<.25 and (k==0 or k==s) THEN
  x := sqrt(2*(h*=e)) * (1-2*k*k);
  x *= sqrt(1 + x/3);
  REPEAT
    y := 1+x;
    z := ln(y) - (y-1-x)/y; // = lnp1(x);
    x := (x-z+h)/z;
    z := 1+x;
  UNTIL abs((z-y)/z) < 1e-8;
  RETURN (h-1)/z;
END;
A, T := abs(a), (2.*k*pi+arg(a))*i;
IF A-A==undef THEN RETURN inf END;
IF A==0 THEN RETURN k ? -inf : 0 END;

bad := k ? -sign(k) : s;
cln := x -> BEGIN x:=ln(x); sign(im(x))==bad ? conj(x) : x END;
x := T;         // rough guess
IF k==0 THEN x := (A+100)/(A+50); x := lnp1(a*x)/x END;
IF k==s THEN x := cln(-a) + (A<.5 ? 0 : T/2) END;
T += log(A);    // = ln_k(a)
REPEAT
  h := x/(x+1) * (x + cln(x) - T);
  x -= h;
UNTIL abs(h/x) < 1e-6;
h := (x-a*exp(-x))/(x+1);   // final touch
RETURN h-h==undef ? x : x-h;
END;
#end

Cas> lambertw(-0.3678)      → −0.979360714958
Cas> lambertw(-0.3678,-1)  → −1.02092723941
Cas> lambertw(3+4i, 5)       → −1.81700589185+30.713334137*i
Cas> lambertw(3+4i, -5)      → −1.75476362795−28.8571010474*i