Lambert W Function (hp-42s)

[Werner]

(09-30-2020 05:28 PM)Albert Chan Wrote:  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.

It's rather the y+x that is the culprit - or y*ln(y)-x, which is zero as you pointed out.
LN(y)+1 for y close to 1/e can be calculated precisely as (LN1P being the LN1+X function)
LN(y) + 1 = LN((1/e)*(1+e*(y-1/e)) + 1 = -1 + LN1P(e*(y-1/e)) + 1
= LN1P(e*(y-1/e))
eg. y = 1/e + 1e-17
then LN(y) + 1 = 2.71828182845904521e-17
LN1P(e*(y-1/e)) = 2.718281828459045198415006976699411e-17

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

Cheers, Werner