Lambert W

# 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

( 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

( 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