(28/48/50) Lambert W Function

[Albert Chan]

The issue is how to accurately calculate Newton's iteration, around branch point, -1/e
I prefer solving for y = e^W(a). Its Newton's method is very simple.

x*e^x = ln(y)*y = a
f = y*ln(y) - a

y = y - f/f' = y - (y*ln(y) - a) / (ln(y)+1)
y = (y+a) / (ln(y)+1)

Around a ≈ -1/e, (y, a) are considered input, thus "exact".
They are about the same size, with opposite sign.
By Sterbenz lemma, (y+a) is exact.

e^W(-1/e) = 1/e      // a ≈ -1/e --> y ≈ +1/e

(12-30-2022 11:48 PM)Albert Chan Wrote:  ln(y)+1 = ln(e*y) = log1p(e*(y-1/e))

For accurate denominator, we need "doubled" precision for 1/e
For 12 digits, 1/e ≈ 0.367879441171 + 0.442321595524E-12

To reduce iterations, we also need good starting estimate.
Lua e^W code (0, -1 branch): https://www.hpmuseum.org/forum/thread-19...#pid167919