Lambert W Function (hp-42s)

[Albert Chan]

(10-05-2020 03:20 PM)Albert Chan Wrote:  Tracking LASTX is hard, it might be a good idea to use it up ASAP, like this ...

(12-29-2022 04:17 AM)Albert Chan Wrote:  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.

This updated version fixed both issues at the same time.
Termination test is now: y' == y' * (1 + relative_error^2), independent of hardware.

Code:

00 { 56-Byte Prgm }
01▸LBL "eW"
02 3            @   X        Y        Z        T
03 1/X
04 -1
05 E^X @        @   r       1/3       x
06 RCL+ ST Z    @  r+x
07 STO× ST Y    @       (r+x)/3
08 LASTX        @   r       r+x    (r+x)/3     x
09 STO+ ST Z
10 STO+ ST X    @  2r       r+x    (r+x)/3+r   x
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↑
25 RCL+ ST Z
26 X<>Y         @ ln(y)+1   y+x     y       x
27 ÷
28 -
29 LASTX        @  y'       dy      x       x
30 ÷
31 STO× ST X
32 LASTX        @  y'      err^2    x       x
33 STO× ST Y
34 STO+ ST Y    @  y'  y*(1+err^2)  x       x
35 GTO 01
36 END