lambertw, all branches, Cas code

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