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

[Albert Chan]

(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)

[Albert Chan]

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]