Lambert W Function (hp-42s)

[Albert Chan]

(10-02-2020 03:02 PM)Werner Wrote:  eg. y = 1/e + 1e-17
then   LN(y) + 1 = 2.71828182845904521e-17
LN1P(e*(y-1/e)) = 2.718281828459045198415006976699411e-17

Above amazing accurate results are misleading.
You have y-1/e = 1e-17 (exactly). In other words, y has infinite number of digits ...

For rounded 34 digits of y, Free42 will evaluate to the same result.
However, both answers only matched half precision. (log1p version does not help here)

-1 [EXP] 1e-17 [+]       // 3.678794411714423315955237701614609e-1 = y

To get an accurate slope, log(y)+1 = log1p(ε = e*y-1) = log1p(ε = e*(y-1/e))
But, this shift the problem to get accurate ε, which required more precise 1/e

1e-17         // y - 1/e
-3.255418886896823216549216319830254E-35             // (more precise 1/e) - 1/e
−             // 1.000000000000000003255418886896823e-17 = y - (more precise 1/e)
1 [EXP] ×     // 2.718281828459045244209433475626668e-17 = ε
[LN1+X]       // 2.718281828459045207264152980973417e-17 = log1p(ε)

Mathematica gives 2.718281828459045207264152980973418e-17

Quote:But the cancellation happens also in y + x, and there's nothing we can do. Or at least, nothing *I* can do ;-)

Since y ≈ -x, (y + x) is exact, without loss of precision
(both x, y are inputs, thus considered exact)