Lambert W Function (hp-42s)

[Werner]

Ah yes, my simple program of course only works for reals.
Okay then, let's use a general purpose complex solver.
I converted the one written by Valentin Albillo for the HP-35S, to Free42/DM42.
(you can find it here on his site)
Of course, it is equally usable on the 42S if you change all LSTO statements to STO, and change the 1E-15 to 1E-5.
I changed three things in the program:
- a more accurate version of the root of the quadratic
- a simpler test to see whether the number is sufficiently close to a real (that is not possible on a 35S or Valentin would've used it)
- and just taking the real part of the current estimate if the test is positive

Usage:
ALPHA: name of function to solve
X: guess
XEQ "CSLV"

The guess may be real or complex, and the resulting root as well.

Code:

00 { 145-Byte Prgm }
01▸LBL 01
02 RCL+ "X"
03 ASTO ST L
04 GTO IND ST L
05▸LBL "CSLV"
06 0
07 LSTO "U"
08 LSTO "V"
09 ENTER
10 COMPLEX
11 +
12 LSTO "X"
13 1ᴇ-15
14 LSTO "S"
15▸LBL 02
16 CLX
17 XEQ 01
18 STO+ ST X
19 STO "U"
20 RCL "S"
21 +/-
22 XEQ 01
23 STO "V"
24 RCL "S"
25 XEQ 01
26 ENTER
27 RCL+ "V"
28 RCL- "U"
29 RCL "S"
30 X↑2
31 ÷
32 X<>Y
33 RCL- "V"
34 RCL "S"
35 STO+ ST X
36 ÷
37 ÷
38 RCL "U"
39 LASTX
40 ÷
41 STO ST Z
42 ×
43 1
44 -
45 +/-
46 SQRT
47 1
48 +
49 ÷
50 STO- "X"
51 RCL÷ "X"
52 ABS
53 RCL "S"
54 X↑2
55 X<Y?
56 GTO 02
57 RCL "X"
58 SIGN
59 COMPLEX
60 X<>Y
61 R↓
62 ABS
63 X>Y?
64 GTO 00
65 RCL "X"
66 COMPLEX
67 R↓
68 STO "X"
69▸LBL 00
70 RCL "X"
71 END

So, to find x=i^i^i..., solve f(x) = i^x - x = 0

define FX:

Code:

>LBL "FX"
 -1
 SQRT
 X<>Y
 Y^X
 LASTX
 -
 END

"FX"
0
XEQ "CSLV"
(0.43829..,0.36059..)
XEQ "FX" results in (0,-1e-34)

Implement Lambert's W:

Code:

>LBL "W"
 LSTO "W"
 "WX"
 1.4
 X<>Y
 REAL?
 +
 XEQ "CSLV"
 RTN
>LBL "WX"
 E^X
 LASTX
 x
 RCL- "W"
 END

The solution to the above problem can be written as

x = i^x

x = e^(x*ln(i))

x*ln(i) = ln(i)*e^(x*ln(i))

ln(i) = x*ln(i)*e^(-x*ln(i))

Let y = -x*ln(i), then solve

y^*e^y = -ln(i) or calculate W(-ln(i))

-1
SQRT
LN
+/-
XEQ "W"
-1
SQRT
LN
+/-
/

results in the almost the the same number as above.

However, Houston, we have a problem here. The definition of WX must not use the local variables used by CSLV, namely X, S, U and V. If we had used variable X instead of W, the routine would not have worked. To remove this dependency, here's a version of CSLV that does not use alphanumeric variables, except REGS:
(making it less usable on a real 42S because it will overwrite REGS)

00 { 115-Byte Prgm }
01▸LBL 01
02 RCL 00
03 +
04 ASTO ST L
05 GTO IND ST L
06▸LBL "CSLV"
07 2
08 ENTER
09 NEWMAT
10 ENTER
11 COMPLEX
12 +
13 LSTO "REGS"
14 1ᴇ-15
15 STO 01
16▸LBL 02
17 CLX
18 XEQ 01
19 STO+ ST X
20 STO 02
21 RCL 01
22 +/-
23 XEQ 01
24 STO 03
25 RCL 01
26 XEQ 01
27 ENTER
28 RCL+ 03
29 RCL 02
30 -
31 RCL 01
32 X↑2
33 ÷
34 X<>Y
35 RCL 03
36 -
37 RCL 01
38 STO+ ST X
39 ÷
40 ÷
41 RCL 02
42 LASTX
43 ÷
44 STO ST Z
45 ×
46 1
47 -
48 +/-
49 SQRT
50 1
51 +
52 ÷
53 STO- 00
54 RCL 00
55 ÷
56 ABS
57 RCL 01
58 ABS
59 X↑2
60 X<Y?
61 GTO 02
62 RCL 00
63 SIGN
64 COMPLEX
65 X<>Y
66 R↓
67 ABS
68 X>Y?
69 GTO 00
70 RCL 00
71 COMPLEX
72 R↓
73 STO 00
74▸LBL 00
75 RCL 00
76 END

Hope you like it, Werner