[HP41] Lambert function RPN; question

[Pekis]

Here's a simple algorithm derived from the one I found on math.stackexchange for both branches of the Lambert function:

The iteration
Using Newton's method to solve x * exp(x) = p yields the following iterative step for finding x = W(p):
x' = (p * exp(−x) + x * x) / (x + 1)

Initial values x0
For the W0 branch,
When p is in [−1/e, 10], use x0 = 0
When p>10, use x0 = ln(p)−ln(ln(p)), ie following asymptotic behavior

For the W-1 branch,
When p is in [−1/e,−0.1], use x0 = −2
When p is in ]−0.1,0], use x0 = ln(−p)−ln(−ln(−p)), ie following asymptotic behavior