[HP41] Lambert function RPN; question

[Albert Chan]

When I try to understand lyuka's eW formula (previous post), I found a better one.

e^W(a) = y = (y+a) / (ln(y)+1)

Let r = 1/e
If a = -r + h (tiny h), we have y = r + z (tiny z)

e^W(a) = (r+z) = (h+z) / ln(1+e*z)

Let H = e*h (known), Z = e*z (unknown), we have:

(1+Z) = (H+Z) / ln(1+Z)

H = (1+Z) * ln(1+Z) - Z = Z^2/2 - Z^3/6 + Z^4/12 - Z^5/20 + ...

2H ≈ Z^2 / (1+Z/3)
Z^2 - 2*(H/3)*Z - 2H ≈ 0
Z ≈ H/3 ± √((H/3)^2 + 2H) ≈ H/3 ± √(2H)

Divide both side by e, we have: y = r + z ≈ r + h/3 ± √(2*r*h)

We could estimate Z another way, by first ignoring O(Z^3) terms
Note: +sign for branch 0, -sign for branch -1

2H ≈ Z^2                      → Z = ±√(2H)

Then, we solve again, using above rough estimated Z

2H ≈ Z^2 / (1+Z/3)      → Z = ±√(2H*(1+Z/3))

Redo previous examples, with this 2 step Z estimate, then convert back to y = r + r*Z

e^W₀ (-0.367879) ≈ 0.368449 321311
e^W_-1(-0.367879) ≈ 0.367309 855146
e^W₀ (-0.3678) ≈ 0.375551 106237
e^W_-1(-0.3678) ≈ 0.360260 737204