lambertw, all branches, Cas code
02-04-2024, 07:34 PM (This post was last modified: 02-04-2024 10:55 PM by Albert Chan.)
Post: #1
 Albert Chan Senior Member Posts: 2,677 Joined: Jul 2018
lambertw, all branches, Cas code
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

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
 « Next Oldest | Next Newest »

 Messages In This Thread lambertw, all branches, Cas code - Albert Chan - 02-04-2024 07:34 PM RE: lambertw, all branches, Cas code - Albert Chan - 02-05-2024, 12:46 AM RE: lambertw, all branches, Cas code - Albert Chan - 02-06-2024, 01:44 PM

User(s) browsing this thread: 1 Guest(s)