hp41programs

Overview
 

 1°)  Complex Operations
 2°)  Complex Equations

   a)  Secant Method
   b)  Quadratic Interpolation
   c)  Newton's Method
   d)  A Better Method
 

1°) Complex Operations
 

-The global labels  "Z+Z"  "Z-Z"  "Z*Z"  "Z/Z"   mean addition, subtraction, multiplication and division.
  "Z^2"  "SQRTZ"  "LNZ"  "E^Z"  "1/Z"   mean  square, square root, logarithm, exponential, reciprocal.
  "SINZ"  "COSZ"  "TANZ"  "ASINZ"  "ACOSZ"  "ATANZ"  correspond to sine, cosine, tangent and their inverses.
  "SHZ"  "CHZ"  "THZ"  "ASHZ"  "ACHZ"  "ATHZ"  correspond to hyperbolic sine, hyperbolic cosine, hyperbolic tangent and their inverses:
-I've used the French notations ( sh z instead of sinh z  ... etc ... ) which have the advantage of saving several bytes and keystrokes.
  "Z^X" and "Z^Z"  stand for raising a complex number to a real power and a complex power respectively.

-These programs use no data register and no external subroutine.
-However, ArcCosh z clear synthetic registers M &N - which can be replaced by R00 & R01.
-And, for instance, "TANZ"  calls  "THZ"  as a subroutine since  tan z = i.tanh (-i.z) and similarly for many other complex functions.
-In these cases, local labels have been inserted after the global labels  ( for example,  XEQ 03  LBL 03 instead of  XEQ "THZ",
  which saves 1 byte and executes faster ) and so on ...
-Lines 144 & 156 are only useful te restore the DEG mode: they may be deleted if RAD is your favorite angular mode.

-To get the magnitude of a complex number, simply press  R-P.

 Formulae:

    Sinh z  = ( e^z - e^-z )/2                           ArcSinh z  = Ln [ z + ( z² + 1 )^1/2 ]
    Cosh z = ( e^z + e^-z )/2                          ArcCosh z = Ln [ z + ( z + 1 )^1/2 ( z - 1 )^1/2 ]
    Tanh z = ( e^2z - 1 )/( e^2z + 1 )               ArcTanh z = [ Ln ( 1 + z ) - Ln ( 1 - z ) ]/2

    Sin z  = -i Sinh iz                                  ArcSin z = -i ArcSinh iz
    Cos z =  Cosh iz                                  ArcCos z = PI/2 - ArcSin z
    Tan z =  i Tanh -iz                                ArcTan z = i ArcTanh -iz
 

Data Registers: /
Flags: /
Subroutines:  /
 
 

 
    ( 448 bytes / SIZE 000 )
 

1°) First case:  Function of one complex number
 

where   z = x+i.y  ;  f(z) = u+i.v

Example:       z = 2+3.i

-Calculate  z² ;  z^1/2  ;  1/z  ;  exp(z)  ;  Ln(z)  ;  sin z  ;  cos z  ;  tan z ;  Asin z  ; Acos z ; Atan z ;  Sinh z ; Cosh z ; Tanh z ; Asinh z ; Acosh z  ; Atanh z

3  ENTER^   2   XEQ "Z^2"        >>>      -5      X<>Y      12        whence    (2+3.i)²   =   -5+12.i
3  ENTER^   2   XEQ "SQRTZ"  >>>   1.6741  X<>Y   0.8960   whence    (2+3.i)^1/2 =   1.6741+0.8960.i
3  ENTER^   2   XEQ "1/Z"         >>>   0.1538  X<>Y  -0.2308   whence   1/(2+3.i)  =    0.1538-0.2308.i
3  ENTER^   2   XEQ "E^Z"        >>>  -7.3151  X<>Y   1.0427   whence   exp(2+3.i) =  -7.3151+1.0427.i
3  ENTER^   2   XEQ "LNZ"       >>>   1.2825  X<>Y    0.9828   whence   Ln(2+3.i)  =   1.2825+0.9828.i
3  ENTER^   2   XEQ "SINZ"     >>>    9.1545  X<>Y   -4.1689  whence   Sin(2+3.i)  =   9.1545-4.1689.i
3  ENTER^   2   XEQ "COSZ"    >>>   -4.1896  X<>Y  -9.1092  whence   Cos(2+3.i) =  -4.1896-9.1092.i
3  ENTER^   2   XEQ "TANZ"    >>>   -0.0038  X<>Y   1.0032  whence    Tan(2+3.i) =  -0.0038+1.0032.i
3  ENTER^   2   XEQ "ASINZ"   >>>    0.5707  X<>Y   1.9834  whence   Asin(2+3.i) =   0.5707+1.9834.i
3  ENTER^   2   XEQ "ACOSZ"  >>>   1.0001  X<>Y  -1.9834  whence   Acos(2+3.i) =  1.0001-1.9834.i
3  ENTER^   2   XEQ "ATANZ"  >>>   1.4099  X<>Y   0.2291  whence    Atan(2+3.i) =  1.4099+0.2291.i
3  ENTER^   2   XEQ "SHZ"        >>>  -3.5906  X<>Y  0.5309  whence    Sinh(2+3.i)  =  -3.5906+0.5309.i
3  ENTER^   2   XEQ "CHZ"       >>>   -3.7245  X<>Y  0.5118  whence    Cosh(2+3.i) = -3.7245+0.5118.i
3  ENTER^   2   XEQ "THZ"       >>>    0.9654   X<>Y  -0.0099  whence   Tanh(2+3.i)  =  0.9654-0.0099.i
3  ENTER^   2   XEQ "ASHZ"    >>>    1.9686   X<>Y   0.9647  whence   Asinh(2+3.i) =  1.9686+0.9647.i
3  ENTER^   2   XEQ "ACHZ"   >>>     1.9834   X<>Y  1.0001  whence    Acosh(2+3.i) = 1.9834+1.0001.i
3  ENTER^   2   XEQ "ATHZ"    >>>    0.1469   X<>Y   1.3390  whence    Atanh(2+3.i) =  0.1469+1.3390.i

Note:    "Z^2"  "SQRTZ"  "1/Z"  "E^Z"  and  "LNZ"   save registers Z and T ( this feature is important when they are called as subroutines by the other functions ).
 

2°) Second case:  raising a complex number to a real power
 

where   z = x+i.y  ;  z^r = u+i.v

Example:   Evaluate  (2+3.i)^1/5

3  ENTER^   2  ENTER^   5   1/X   XEQ "Z^X"   >>>  1.2675   X<>Y   0.2524    whence   (2+3.i)^1/5 = 1.2675+0.2524.i

Note:   "Z^X"  saves T-register in  registers Z and T.
 

3°) Third case:  Function of two complex numbers.
 

where   z = x+i.y  ;  z' = x'+i.y' ; f(z;z') = u+i.v

Example:   z = 2+3.i  ;  z' = 4+7.i  Compute  z+z' ; z-z' ; z.z' ; z/z' ; z^z'
 

3  ENTER^   2  ENTER^   7  ENTER^   4  XEQ "Z+Z"    gives      6       X<>Y       10      whence    z+z' = 6+10.i
3  ENTER^   2  ENTER^   7  ENTER^   4  XEQ "Z-Z"     gives     -2       X<>Y       -4      whence    z-z' = -2-4.i
3  ENTER^   2  ENTER^   7  ENTER^   4  XEQ "Z*Z"    gives   -13       X<>Y       26      whence    z.z' = -13+26.i
3  ENTER^   2  ENTER^   7  ENTER^   4  XEQ "Z/Z"     gives   0.4462  X<>Y  -0.0308   whence   z/z' = 0.4462-0.0308.i
3  ENTER^   2  ENTER^   7  ENTER^   4  XEQ "Z^Z"    gives  0.1638   X<>Y   0.0583    whence    z^z' = 0.1638+0.0583.i
 

Notes:

-This program may of course be simplified if you can use the M-code functions:

    1/Z  Z*Z  Z/Z  Z^2  SQRTZ  which are listed in "A few M-code Routines for the HP-41"

-Slightly different formulae may also be used to compute  ArcCosh  ArcTanh and ArcCos:

    Sinh z  = ( e^z - e^-z )/2                           ArcSinh z  = Ln [ z + ( z² + 1 )^1/2 ]
    Cosh z = ( e^z + e^-z )/2                          ArcCosh z = Ln [ z + ( z² - 1 )^1/2 ]
    Tanh z = ( e^2z - 1 )/( e^2z + 1 )               ArcTanh z = [ Ln (( 1 + z ) / ( 1 - z )) ]/2

    Sin z  = -i Sinh iz                                  ArcSin z = -i ArcSinh iz
    Cos z =  Cosh iz                                  ArcCos z = -i ArcCosh z
    Tan z =  i Tanh -iz                                ArcTan z = i ArcTanh -iz

-These definitions are possible and easier to program since no synthetic register is necessary.
-But they don't always give the principal value of some inverse hyperbolic functions.
-Anyway, here is this simpler program:
 

Data Registers: /
Flags: /
Subroutines:  /
 
 

 
    ( 415 bytes / SIZE 000 )
 

-Same instructions and - in most cases - same results...
 

2°) Complex Equations:   f(z) = 0
 

        a) The "Secant" Method
 

-This program uses an iterative algorithm.
-Starting with 2 distinct initial guesses  z₀ = x₀+i.y₀ ;  z₁ = x₁+i.y₁ ,  it calculates  z_k+1 = z_k - f(z_k).(z_k - z_k-1) / [ f(z_k) - f(z_k-1) ]     ( k = 1 ; 2 ; 3 ; ....... )
-The process is repeated until 2 consecutive approximations are equal.

Data Registers:     • R00 = Function name   ( Global label, maximum of 6 characters. This register is to be initialized before executing "ROOTZ" )

                                R01 = x   R03 = a   R05 = x'
                                 R02 = y   R04 = b   R06 = y'   where  f(z) = f(x+i.y) = a+i.b
Flag: /
Subroutine:  A program which computes  f(z) = a+i.b  assuming  x is in X-register and y is in Y-register upon entry
                       In other words, the stack      Y = y           Y = b
                        must be changed from          X = x    to    X = a      where z = x+i.y   and   f(z) = f(x+i.y) = a+i.b
 
 

 
    ( 86 bytes / SIZE 007 )
 

 
Example:   Find a root of      sinh z + z² + Pi = 0.      First, we load a routine to compute f(z) , for instance:
 

  LBL "T"
  STO 07
  X<>Y
  STO 08
  X<>Y
  XEQ "SHZ"
  RCL 08
  RCL 07
  XEQ "Z^2"
  XEQ "Z+Z"
  PI
  +
  RTN

And if we choose   z₀ = 0  ;  z₁ = 1+i

   T  ASTO 00
   0  ENTER^   ENTER^
   1  ENTER^   XEQ "ROOTZ"  the succesive x-values are displayed

        and we get   -0.278189856                  ( execution time = 77 seconds )
            X<>Y      1.812880366

Whence     -0.278189856 + 1.812880366.i    is a solution.   ( There are many other solutions like  4.419361546 + 4.697881722.i  ... etc ... )

Note:

-The iterations stop when  z_k+1 =  z_k  or when   f(z_k+1) =  f(z_k)  Therefore, check f(z) in registers R03 & R04
 

        b) Quadratic Interpolation
 

-Starting with 3 initial guesses:    z₀ = x₀+i.y₀ ;  z₁ = x₁+i.y₁ ;  z₂ = x₂+i.y₂ , this program calculates

   A = ( z₁-z₀ ).f(z₂) + ( z₀-z₂ ).f(z₁) + ( z₂-z₁ ).f(z₀)
   B = ( z₁-z₀ ).( 2.z₂-z₁-z₀ ).f(z₂) - ( z₀-z₂ )².f(z₁) + ( z₂-z₁ )².f(z₀)
   C = ( z₂-z₁ ).( z₁-z₀ ).( z₂-z₀ ).f(z₂)

and   z₃ = z₂ - 1 / [ B/(2.C) + E.( B²/(4C²) - A/C )^1/2 ]      where   E = +1 or -1   ( the sign is choosen to produce the larger denominator )

-The process is repeated until  | z_k+1-z_k |  is smaller than  10^-8  or  C = 0   ( check  f(z) in  R03 & R04 )

Data Registers:   • R00 = Function name   ( Global label, maximum of 6 characters. This register is to be initialized before executing "RTZ2" )

                                 R01 = x   R03 = a   R05 thru R20:  temp
                                 R02 = y   R04 = b   where  f(z) = f(x+i.y) = a+i.b
Flag: /
Subroutines:  "Z+Z"  "Z*Z"  "Z^2"  "SQRTZ"

  and a program which computes  f(z) = a+i.b  assuming  x is in X-register and y is in Y-register upon entry

       The stack must       Y = y           Y = b      where     z = x+i.y
       be changed from     X = x    to    X = a        and   f(z) = f(x+i.y) = a+i.b

 ( Registers R03 R04 R13 R14 ... etc ... may be used by this subroutine )

-Line  58  is a three-byte  GTO 02
-Line 170 is a three-byte  GTO 01
 
 

 
    ( 281 bytes / SIZE 021 )
 
 

( The 3 initial guesses are actually   z₀ = x₀+i.y₀ ;   z₀-(1+i).h ; z₀-2.(1+i).h  )

Example:     With the same equation     sinh z + z² + Pi = 0.     We load a routine to compute f(z)
 

  LBL "T"
  STO 03
  X<>Y
  STO 04
  X<>Y
  XEQ "SHZ"
  RCL 04
  RCL 03
  XEQ "Z^2"
  XEQ "Z+Z"
  PI
  +
  RTN

And if we choose   z₀ = 1+i  and  h = 1

   T  ASTO 00
   1  ENTER^   ENTER^
       XEQ "RTZ2"  the succesive x-values are displayed and we get   -0.278189857                  ( execution time = 156 seconds )
                                                                                          X<>Y      1.812880366

Whence    z = -0.278189857 + 1.812880366.i

Note:     In this example,  "ROOTZ"  is faster than  "RTZ2"   but if  f is a very complicated function, the advantage is clearly with "RTZ2"
 

        c) Newton's Method
 

-Only one initial guess  z₀ = x₀+i.y₀  is needed but we have to compute the first derivative f '(z)

   Formula:      z_k+1 = z_k - f(z_k)/f '(z_k)

Data Registers:     • R00 = f name       • R03 = f ' name    ( these 2 registers are to be initialized before executing "NWTZ" )

                                 R01 = x
                                 R02 = y   R04 = | f '(z) |   R05: temp
Flag: /
Subroutines:

-A program which computes  f(z)   assuming  x is in X-register and y is in Y-register upon entry
-a program which computes f '(z)   assuming  x is in X-register and y is in Y-register upon entry
 
 
 

 
    ( 55 bytes / SIZE 006 )
 

 -The output in Z-register should be a "small" number  if  z = x+i.y  is a "true" root

Example:       sinh z + z² + Pi = 0.     We have to load 2 routines to compute f(z) and f ' (z) = cosh z + 2.z  for instance:
 

  LBL "T"
  STO 07
  X<>Y
  STO 08
  X<>Y
  XEQ "SHZ"
  RCL 08
  RCL 07
  XEQ "Z^2"
  XEQ "Z+Z"
  PI
  +
  RTN
  LBL "U"
  STO 07
  X<>Y
  STO 08
  X<>Y
  XEQ "CHZ"
  RCL 08
  ST+ X
  RCL 07
  ST+ X
  XEQ "Z+Z"
  RTN

Then    T  ASTO 00   U  ASTO 03   and with  z₀ = 1+i
           1  ENTER^    XEQ "NWTZ"  the successive  x-approximations are displayed
                                      and we get  -0.278189857
                                              RDN    1.812880366     whence  z = -0.278189857 + 1.812880366.i      ( execution time = 54 seconds )

Note:   | f(z) |  in Z-register actually concerns the next to last approximation.
 

        d) A Better method
 

-Starting with one initial guess   z₀ = x₀+i.y₀  , we now use the first and second derivatives   f '(z)  and  f "(z)  ,  the recursion formula is:
   z_k+1 = z_k -  1 / [ f '(z_k)/f(z_k) + E.( ( f(z_k)/(2.f '(z_k)) )² - f "(z_k) /(2.f (z_k)) )^1/2 ]   where   E = +1 or -1   ( the sign is choosen to produce the larger denominator )

Data Registers:   • R00 = f name       • R03 = f ' name     • R04 = f " name  ( these 3 registers are to be initialized before executing "RTZ4" )

                                 R01 = x
                                 R02 = y   R05 = | f '(z) |    R06 thru R09: temp
Flag: /
Subroutines:

-A program which computes  f(z)    assuming  x is in X-register and y is in Y-register upon entry
-A program which computes  f '(z)   assuming  x is in X-register and y is in Y-register upon entry
-A program which computes  f "(z)   assuming  x is in X-register and y is in Y-register upon entry
 
 

 
    ( 123 bytes / SIZE 010 )
 

 -The output in Z-register should be a "small" number  if  z = x+i.y  is a "true" root

Example:       sinh z + z² + Pi = 0.     We have to load 3 routines to compute f(z) ; f ' (z) = cosh z + 2.z and  f "(z) = sinh z + 2  ,  for instance:
 

   LBL "T"
   STO 10
   X<>Y
   STO 11
   X<>Y
   XEQ "SHZ"
   RCL 11
   RCL 10
   XEQ "Z^2"
   XEQ "Z+Z"
   PI
   +
   RTN
   LBL "U"
   STO 10
   X<>Y
   STO 11
   X<>Y
   XEQ "CHZ"
   RCL 11
   ST+ X
   RCL 10
   ST+ X
   XEQ "Z+Z"
   RTN
   LBL "V"
   XEQ "SHZ"
   2
   +
   RTN
 

Then    T  ASTO 00   U  ASTO 03   V ASTO 04  and with  z₀ = 1+i
           1  ENTER^    XEQ "RTZ4"    the successive  x-approximations are displayed
                                      and we get  -0.278189857
                                              RDN    1.812880366     whence  z = -0.278189857 + 1.812880366.i      ( execution time = 70 seconds )

Notes:

   | f(z) |  in Z-register actually concerns the next to last approximation.
   Among these 4 root-finding programs, "RTZ4" has the best rate of convergence.