Lambert W Function (hp-42s)

[Albert Chan]

(09-30-2020 09:12 AM)Werner Wrote:  That is a direct consequence of using y'=y'+(y'-y)^2 as a stopping criterion..
You should do 1 extra Newton-Rhapson step

If y' = y' + (y'-y)^2, it implied y has at least half-precision.
Assuming Newton's step doubled its precision, y' has about full precision.

However, this does not apply on singularity. For y = e^W(x), the newton formula is:

y ← y - (y*ln(y) - x) / (ln(y) + 1)

Or, equivalent version (simpler, but slightly less accurate): y ← (y + x) / (ln(y) + 1)

As y approach 1/e, slope (denominator) goes to zero.
With limited precision, as soon as y reached half precision, ln(y) + 1 will lose half precison too.
First Newton formula suffered the same issue, with questionable correction term.

Both form we are hit with 0/0 issue. However, limit do exist:

XCas> x0 := -1/e
XCas> limit([y - (y*ln(y)-x0)/(ln(y)+1), (y+x0)/(ln(y)+1)], y=1/e)       → [1/e, 1/e]

Let's try e^W(10^-34 - 1/e), using lyuka "eW":

-.3678794411714423215955237701614608       x = 10^-34 - 1/e = 10^-34 - r
 .3678794411714423301731626197685289       g = r + sqrt(2*r*(x+r)) + 0.3*(x+r)
 .3678794411714423295023829145484696       y = newton(s) w/ guess g
 .3678794411714423286399438752146243055... Mathematica for e^W(x)

Newton's step(s) gives back only 17 good digits (as expected)

Had we apply 1 more newton step (with last y)

-.9999999999999999785069284676275602       ln(y)
 .0000000000000000214930715323724398       ln(y) + 1
-.3678794411714423215955237701614608       y ln(y)

We have only half-precision slope. Also, y ln(y) - x = 0. Newton steps does nothing ...