lambertw, all branches, Cas code - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Software Libraries (/forum-10.html) +--- Forum: HP Prime Software Library (/forum-15.html) +--- Thread: lambertw, all branches, Cas code (/thread-21240.html) lambertw, all branches, Cas code - Albert Chan - 02-04-2024 07:34 PM Python W, version 4, translated for HP Prime Cas, with some differences. 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   lyuka branch solve with e^W(a) = y = r+r*x, for x see https://www.hpmuseum.org/forum/thread-19689-post-183435.html#pid183435   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 RE: lambertw, all branches, Cas code - Albert Chan - 02-05-2024 12:46 AM (02-04-2024 07:34 PM)Albert Chan Wrote:  z := ln(y) - (y-1-x)/y; // = lnp1(x); Technically lnp1(x) == ln(1+x) + correction lnp1(x) might even be implemented this way, see latest Free42 function C.LN1+X However, HP Prime have trouble getting accurate result, if x is tiny. lnp1(x) = 2*atanh(y) = 2*(y + y^3/3 + y^5/5 + ...), where y=x/(x+2) Cas> x := -4e-16+5e-8i Cas> 2*(y := x/(x+2)) 8.5e−16+0.00000005*i Cas> 2*atanh(y) 13.3226762955e−16+0.00000005*i Cas> lnp1(x) 12.5e−16+0.00000005*i Cas> y:=x+1 Cas> ln(y) - (y-1-x)/y 9.3226762955e−16+0.00000005*i atanh/lnp1/ln are all bad, but ln(1+x) + correction version seems better (at least for this case) Hopefully, these numerical bugs would get fixed. For comparison, this is HP71B (can't test LOGP1, because it does not support complex type) >complex x,y >x = (-4e-16,5e-8) >y = x/(2+x) >2*y  (8.5E-16,.00000005) >2*ATANH(y)  (8.5E-16,.00000005) >y = 1+x >LN(y) - (y-1-x)/y  (8.5E-16,.00000005) >LN(y)*x/(y-1)                  ! Kahan's log1p way, return x if y=1  (8.5E-16,.00000005) RE: lambertw, all branches, Cas code - Albert Chan - 02-06-2024 01:44 PM I experimented ln / lnp1 for real agrument. First, some reference numbers, using lua lua> p48 = fn'x: local hi,lo = bit.split(x); bit.join(hi, lo-lo%32)' -- Lua to Prime float lua> m48 = fn'x: trunc(frexp(x) * 2^48)' -- HP Prime mantissa bits lua> x = p48(.1) lua> y = p48(1+x) lua> m48(x) 225179981368524 lua> m48(log1p(x)) 214619445125685 lua> m48(log(y)) 214619445125674 lua> m48(p48(log(y)) - p48((y-1-x)/y)) -- simulate HP Prime float operation. 214619445125684 Cas> m48(x) := trunc(x * 2^(48-frexp(x)[2])) Cas> x := .1 Cas> y := 1+x Cas> m48(x) 225179981368524 Cas> m48(lnp1(x)) 214619445125630 Cas> m48(ln(y)) 214619445125674 Cas> m48(ln(y)-(y-1-x)/y) 214619445125684 Strange, lnp1(x) is worse than ln(1+x), even without correction. Update: I tried this on XCas 1.9.0-31 (mingw win32) I was expecting same result, since both use 48-bits truncated float. XCas> log1p(.1) - ln(1+.1)      → 4.88498130835e-15 XCas> ans() * 2^51                → 11. XCas log1p(.1) gives exactly lua reference, 74 + 11 = 85. lnp1(.1) bug is only for HP Prime, error = 85 - 30 = 55 ULP I discovered 2 more HP Prime bugs, during experimentation. 1.) round(x) likely defined as floor(x+0.5), even if x is an integer. I was expecting symmetry, round(x) = sign(x)*round(|x|), like C++ round Not matching others is not a bug, but round an integer to something else is. Cas> m := float(m48(x)) Cas>  round( m)       → 225179981368524 Cas> -round(-m)      → 225179981368523 2.) frexp(x) may lose its mantissa. This is why m48 not defined as trunc(frexp(x)[1] * 2^48) Cas> frexp(1e-12)      → [0., -39]