Lambert W Function (hp-42s)

[Albert Chan]

I discovered a bug, for termination condition.

We used y' == y' + dy^2 termination test, but it is wrong
Assuming Newton's method have quadratic convergence, it should be:

1 == 1 + (relative error)^2
1 == 1 + (dy/y')^2
y'^2 == y'^2 + dy^2

If y' is big, this does not matter, our test implied relative error test too.
It may terminate later than needed, but that's OK.

Problem is if we use this for e^W_-1(-ε), which results in y' even smaller than ε
Now, the flawed termination test quit early.

Example, for ε = 1E-99:

-1e-99 XEQ "eW"    → W₀(-ε) = 1.
- - + * R/S        → 2.201582623716450630677604534034768e-102

This quit too early (not even 1 correct digit!), we use this guess for next iteration:

CLX + R/S          → 4.281027759337267571527542476904921e-102

We now have 3 correct digits. It will take a few more iterations for convergence.

We use this formula, which is very accurate close to x = -1/e

e^W(x) ≈ 1/e + (x+1/e)/3 + sqrt ((2/e)*(x+1/e))

I am too lazy to code relative error test, and just hard coded right eps for the job.
If y and y' matched 17 digits, this implied y' converged to 34 digits.
Worse case, we have 17+ correct digits.

Code:

00 { 59-Byte Prgm }
01▸LBL "eW"
02 3
03 1/X
04 -1
05 E↑X
06 RCL+ ST Z
07 STO× ST Y
08 STO+ ST X
09 LASTX
10 STO+ ST Z
11 ×
12 SQRT
13 +
14 X<>Y
15 +/-
16 X<>Y     @   y    -x    x    x

17▸LBL 01
18 X=Y?
19 RTN
20 STO ST Y @   y     y     x    x
21 LN
22 1
23 +
24 R↑       @   x  LN(y)+1  y    x
25 RCL+ ST Z
26 X<>Y
27 ÷        @   y'    y     x    x
28 -
29 LASTX
30 X<>Y     @   dy   y'     x    x
31 1ᴇ-17
32 ×
33 X<>Y     @   y'   eps    x    x
34 +
35 LASTX    @   y'  y'+eps  x    x
36 GTO 01
37 END

With this updated eW, we get (all digits correct)

e^W_-1(-1E-99) = 4.284330166764374654650589938202028e-102
→ W_-1(-1E-99) = -233.4087152660372939791972288026527