Lambert W Function for the HP-41
Overview
1°) Real Arguments
2°) Complex Arguments
3°) Quaternionic Arguments
-We have to solve x = w(x) exp [ w(x) ]
-The following programs use Newton's method.
1°) Real Arguments
Data Register: R00 = x
Flags: /
Subroutines: /
-Line 17 may be replaced by X#Y?? to avoid an infinite loop
( cf "A Few M-Code Routines for the HP-41" )
| 01 LBL "LAMB"
02 STO 00 03 LN1+X 04 LBL 01 05 ENTER^ 06 ENTER^ 07 CHS 08 E^X 09 RCL 00 10 * 11 - 12 X<>Y 13 1 14 + 15 / 16 - 17 X#Y? 18 GTO 01 19 END |
( 29 bytes / SIZE 001 )
| STACK | INPUT | OUTPUT |
| X | x | w(x) |
Examples:
-1/e XEQ "LAMB" >>>> w(-1/e) = -1
-0.1 R/S
>>>> w(-0.1) = -0.111832559
0
R/S
>>>> w(0) = 0
e
R/S
>>>> w(e) = 1
100 R/S
>>>> w(100) = 3.385630140
E99 R/S
>>>> w(1099) = 222.55007690
Notes:
-The initial guess has been changed to produce the standard definition
when -1/e < x < 0
and to avoid divergence for small positive arguments.
-In fact, "LAMB" now returns w(-1/e) = -1.000004459 instead
of -1
-So, you could add 1 CHS E^X
CHS X<>Y X#Y? GTO 00
LASTX RTN LBL 00 after line 02
-If x < -1/e use the routine hereunder, because the result
is a complex number.
2°) Complex Arguments
Data Registers: R01 R02 = z ,
R03 R04 = w , R00-R05: temp
Flags: /
Subroutines: "LNZ" "E^Z" "Z*Z"
"Z/Z" ( cf "Complex Functions for the HP-41" ) or similar M-Code
routines.
-The successive u-approximations are displayed. Otherwise, delete line
54.
| 01 LBL "LAMBZ"
02 STO 01 03 X<>Y 04 STO 02 05 X<>Y 06 1 07 STO 03 08 STO 04 09 + 10 X=0? 11 X#Y? 12 FS?30 |
13 GTO 01
14 XEQ "LNZ" 15 STO 03 16 X<>Y 17 STO 04 18 LBL 01 19 RCL 04 20 RCL 03 21 XEQ "E^Z" 22 STO 05 23 X<>Y 24 STO 00 |
25 X<>Y
26 RCL 04 27 RCL 03 28 XEQ "Z*Z" 29 ST+ 05 30 X<>Y 31 ST+ 00 32 RCL 02 33 - 34 X<>Y 35 RCL 01 36 - |
37 RCL 00
38 RCL 05 39 XEQ "Z/Z" 40 CHS 41 RCL 03 42 + 43 STO 03 44 LASTX 45 - 46 X<>Y 47 CHS 48 RCL 04 |
49 +
50 STO 04 51 LASTX 52 - 53 R-P 54 VIEW 03 55 E-9 56 X<Y? 57 GTO 01 58 RCL 04 59 RCL 03 60 END |
( 84 bytes / SIZE 006 )
| STACK | INPUTS | OUTPUTS |
| Y | y | v |
| X | x | u |
where z = x + i y and w(z) = u + i v
Examples:
9 ENTER^
4 XEQ "LAMBZ" >>>> 1.680157300
---Execution time = 19s---
RDN 0.738465807
whence W(4+9i) = 1.680157300 + 0.738465807 i
0 ENTER^
-1 R/S
>>>> -0.318131505
---Execution time = 21s---
RDN 1.337235701
so
W(-1) = -0.318131505 + 1.337235701 i
3°) Quaternionic Arguments
Data Registers: R00 thru R15: temp , when
the program stops, the solution is in R13-R14-R15-R00
Flags: /
Subroutines: "LNQ" "E^Q" "Q*Q"
"1/Q" ( cf "Quaternions for the HP-41" )
-The M-Code routine X+1 add one to register X
-You can replace line 10 by 1 +
RCL 12 RDN
| 01 LBL "LAMBQ"
02 STO 09 03 RDN 04 STO 10 05 RDN 06 STO 11 07 RDN 08 STO 12 09 RDN 10 X+1 11 XEQ "LNQ" 12 STO 13 13 RDN 14 STO 14 15 RDN 16 STO 15 17 X<>Y 18 STO 00 19 LBL 01 |
20 RCL 00
21 CHS 22 RCL 15 23 CHS 24 RCL 14 25 CHS 26 RCL 13 27 CHS 28 XEQ "E^Q" 29 STO 05 30 RDN 31 STO 06 32 RDN 33 STO 07 34 X<>Y 35 STO 08 36 RCL 09 37 STO 01 38 RCL 10 |
39 STO 02
40 RCL 11 41 STO 03 42 RCL 12 43 STO 04 44 XEQ "Q*Q" 45 STO 05 46 RDN 47 STO 06 48 RDN 49 STO 07 50 X<>Y 51 STO 08 52 RCL 13 53 ST- 05 54 1 55 + 56 RCL 00 57 ST- 08 |
58 RCL 15
59 ST- 07 60 RCL 14 61 ST- 06 62 R^ 63 XEQ "1/Q" 64 STO 01 65 RDN 66 STO 02 67 RDN 68 STO 03 69 X<>Y 70 STO 04 71 XEQ "Q*Q" 72 ST+ 13 73 ABS 74 X<>Y 75 ST+ 14 76 ABS |
77 +
78 X<>Y 79 ST+ 15 80 ABS 81 + 82 X<>Y 83 ST+ 00 84 ABS 85 + 86 VIEW 13 87 E-8 88 X<Y? 89 GTO 01 90 RCL 00 91 RCL 15 92 RCL 14 93 RCL 13 94 END |
( 137 bytes / SIZE 016 )
| STACK | INPUTS | OUTPUTS |
| T | t | t' |
| Z | z | z' |
| Y | y | y' |
| X | x | x' |
with W(x+y.i+z.j+t.k) = x'+y'.i+z'.j+t'.k
Example: q = 1 + 2 i + 3 j + 4 k
4 ENTER^
3 ENTER^
2 ENTER^
1 XEQ "LAMBQ" >>>>> 1.281459811
= R13
---Execution time = 64s---
RDN 0.304051227 = R14
RDN 0.456076840 = R15
RDN 0.608102453 = R00
Whence, W(1 + 2 i + 3 j + 4 k) = 1.281459811 + 0.304051227 i + 0.456076840 j + 0.608102453 k
Notes:
-This routine doesn't work if q = -1 ( cf paragraph 2 )
-Though the multiplication of 2 quaternions is not commutative, we
have: w exp(w) = exp(w) w