lambertw, all branches, Cas code
02-05-2024, 12:46 AM (This post was last modified: 02-05-2024 09:30 PM by Albert Chan.)
Post: #2
 Albert Chan Senior Member Posts: 2,707 Joined: Jul 2018
RE: lambertw, all branches, Cas code
(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)
 « 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: