hp41programs

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<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