# hp41programs

Complex Complex Functions for the HP-41

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  = ( ez - e-z )/2                           ArcSinh z  = Ln [ z + ( z2 + 1 )1/2 ]
Cosh z = ( ez + e-z )/2                          ArcCosh z = Ln [ z + ( z + 1 )1/2 ( z - 1 )1/2 ]
Tanh z = ( e2z - 1 )/( e2z + 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:  /

 01  LBL "SINZ"    02  X<>Y   03  CHS   04  XEQ 01   05  CHS   06  X<>Y   07  RTN   08  LBL "TANZ"   09  CHS   10  X<>Y   11  XEQ 03   12  X<>Y   13  CHS   14  RTN   15  LBL "ASINZ"   16  LBL 04   17  X<>Y   18  CHS   19  XEQ 05   20  CHS   21  X<>Y   22  RTN   23  LBL "ATANZ"   24  CHS   25  X<>Y   26  XEQ 06   27  X<>Y   28  CHS   29  RTN   30  LBL "SHZ"   31  LBL 01   32  XEQ 00   33  XEQ 01   34  GTO 10   35  LBL "COSZ"   36  X<>Y   37  CHS 38  LBL "CHZ"   39  XEQ 00   40  XEQ 02   41  GTO 10   42  LBL "ATHZ"     43  LBL 06   44  RCL X   45  1   46  ST- Z   47  +   48  R^   49  X<>Y   50  XEQ 12   51  R^   52  CHS   53  R^   54  CHS   55  XEQ 12   56  XEQ 01   57  LBL 10   58  2   59  ST/ Z   60  /   61  RTN   62  LBL "THZ"   63  LBL 03   64  2   65  ST* Z   66  *   67  XEQ 08   68  RCL X   69  1   70  ST- Z   71  +   72  R^   73  X<>Y   74  LBL "Z/Z" 75  R-P   76  R^   77  R^   78  R-P   79  R^   80  ST- Z   81  X<> T   82  /   83  P-R   84  RTN   85  LBL 00   86  RCL Y   87  RCL Y   88  XEQ 08   89  R^   90  CHS   91  R^   92  CHS   93  GTO 08   94  LBL "ACOSZ"   95  XEQ 04   96  X<>Y   97  CHS   98  PI   99  2 100  / 101  R^ 102  - 103  RTN 104  LBL "ACHZ" 105  STO M 106  X<>Y 107  STO N 108  X<>Y 109  1 110  + 111  XEQ 07 112  RCL M 113  1 114  - 115  RCL N 116  X<>Y 117  XEQ 07 118  XEQ 03 119  X<>Y 120  0 121  X<> N 122  + 123  X<>Y 124  0 125  X<> M 126  + 127  GTO 12 128  LBL "ASHZ"   129  LBL 05 130  RCL Y 131  RCL Y 132  XEQ 05 133  SIGN 134  ST* X 135  ST+ L 136  X<> L 137  XEQ 07 138  XEQ 02 139  LBL "LNZ" 140  LBL 12 141  RAD 142  R-P 143  LN 144  DEG 145  RTN 146  LBL "Z^Z" 147  R^ 148  R^ 149  XEQ 12 150  XEQ 03 151  LBL "E^Z" 152  LBL 08 153  RAD 154  E^X 155  P-R 156  DEG 157  RTN 158  LBL "Z^X" 159  RDN 160  R-P 161  R^ 162  ST* Z 163  Y^X 164  P-R 165  RTN 166  LBL "Z-Z"      167  LBL 01 168  CHS 169  X<>Y 170  CHS 171  X<>Y 172  LBL "Z+Z" 173  LBL 02 174  X<>Y 175  ST+ T 176  RDN 177  + 178  RTN 179  LBL "Z*Z" 180  LBL 03 181  R-P 182  R^ 183  R^ 184  R-P 185  X<>Y 186  ST+ T 187  RDN 188  * 189  P-R 190  RTN 191  LBL "1/Z" 192  R-P 193  1/X 194  X<>Y 195  CHS 196  X<>Y 197  P-R 198  RTN 199  LBL "Z^2" 200  LBL 05 201  R-P 202  X^2 203  X<>Y 204  ST+ X 205  X<>Y 206  P-R 207  RTN 208  LBL "SQRTZ" 209  LBL 07 210  R-P 211  SQRT 212  SIGN 213  ST+ X 214  ST/ Y 215  X<> L 216  P-R 217  END

( 448 bytes / SIZE 000 )

1°) First case:  Function of one complex number

 STACK INPUTS OUTPUTS Y y v X x u

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

Example:       z = 2+3.i

-Calculate  z2 ;  z1/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)2   =   -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

 STACK INPUTS OUTPUTS Z y / Y x v X r u

where   z = x+i.y  ;  zr = 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.

 STACK INPUTS OUTPUTS T y / Z x / Y y' v X x' u

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' ; zz'

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    zz' = 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  = ( ez - e-z )/2                           ArcSinh z  = Ln [ z + ( z2 + 1 )1/2 ]
Cosh z = ( ez + e-z )/2                          ArcCosh z = Ln [ z + ( z2 - 1 )1/2 ]
Tanh z = ( e2z - 1 )/( e2z + 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:  /

 01  LBL "SINZ"   02  X<>Y   03  CHS   04  XEQ 01   05  CHS   06  X<>Y   07  RTN   08  LBL "TANZ"   09  CHS   10  X<>Y   11  XEQ 03   12  X<>Y   13  CHS   14  RTN   15  LBL "ASINZ"   16  X<>Y   17  CHS   18  XEQ 05   19  CHS   20  X<>Y   21  RTN   22  LBL "ATANZ"   23  CHS   24  X<>Y   25  XEQ 06   26  X<>Y   27  CHS   28  RTN   29  LBL "SHZ"   30  LBL 01   31  XEQ 00   32  XEQ 01   33  GTO 10   34  LBL "COSZ"   35  X<>Y   36  CHS   37  LBL "CHZ"   38  XEQ 00   39  XEQ 02 40  GTO 10   41  LBL "ATHZ"     42  LBL 06   43  RCL X   44  1   45  ST+ Z   46  X<>Y   47  -   48  R^   49  CHS   50  X<>Y   51  XEQ 07   52  XEQ 12   53  LBL 10   54  2   55  ST/ Z   56  /   57  RTN   58  LBL "THZ"   59  LBL 03   60  2   61  ST* Z   62  *   63  XEQ 08   64  RCL X   65  1   66  ST- Z   67  +   68  R^   69  X<>Y   70  LBL "Z/Z"   71  LBL 07   72  R-P   73  R^   74  R^   75  R-P   76  R^   77  ST- Z   78  X<> T 79  /   80  P-R   81  RTN   82  LBL 00   83  RCL Y   84  RCL Y   85  XEQ 08   86  R^   87  CHS   88  R^   89  CHS   90  GTO 08   91  LBL "ACOSZ"   92  XEQ 04   93  CHS   94  X<>Y   95  RTN   96  LBL "ACHZ"   97  LBL 04   98  XEQ 13   99  CHS 100  GTO 09 101  LBL "ASHZ" 102  LBL 05 103  XEQ 13 104  LBL 09 105  ST+ L 106  X<> L 107  XEQ 11 108  XEQ 02 109  LBL "LNZ" 110  LBL 12 111  RAD 112  R-P 113  LN 114  DEG 115  RTN 116  LBL 13 117  RCL Y 118  RCL Y 119  XEQ 05 120  SIGN 121  ST* X 122  RTN 123  LBL "Z^Z" 124  R^ 125  R^ 126  XEQ 12 127  XEQ 03 128  LBL "E^Z" 129  LBL 08 130  RAD 131  E^X 132  P-R 133  DEG 134  RTN 135  LBL "Z^X"      136  RDN 137  R-P 138  R^ 139  ST* Z 140  Y^X 141  P-R 142  RTN 143  LBL "Z-Z" 144  LBL 01 145  CHS 146  X<>Y 147  CHS 148  X<>Y 149  LBL "Z+Z" 150  LBL 02 151  X<>Y 152  ST+ T 153  RDN 154  + 155  RTN 156  LBL "Z*Z" 157  LBL 03 158  R-P 159  R^ 160  R^ 161  R-P 162  X<>Y 163  ST+ T 164  RDN 165  * 166  P-R 167  RTN 168  LBL "1/Z" 169  R-P 170  1/X 171  X<>Y 172  CHS 173  X<>Y 174  P-R 175  RTN 176  LBL "Z^2" 177  LBL 05 178  R-P 179  X^2 180  X<>Y 181  ST+ X 182  X<>Y 183  P-R 184  RTN 185  LBL "SQRTZ" 186  LBL 11 187  R-P 188  SQRT 189  SIGN 190  ST+ X 191  ST/ Y 192  X<> L 193  P-R 194  END

( 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  z0 = x0+i.y0 ;  z1 = x1+i.y1 ,  it calculates  zk+1 = zk - f(zk).(zk - zk-1) / [ f(zk) - f(zk-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

 01  LBL "ROOTZ" 02  STO 01 03  X<>Y 04  STO 02 05  R^ 06  STO 06 07  R^ 08  STO 05 09  XEQ IND 00 10  STO 03 11  X<>Y 12  STO 04 13  LBL 01 14  VIEW 01 15  RCL 02 16  RCL 01 17  XEQ IND 00  18  RCL 04 19  RCL 03 20  R^ 21  STO 04 22  ST- Z 23  R^ 24  STO 03 25  ST- Z 26  R-P 27  R^ 28  R^ 29  R-P 30  X<>Y 31  ST- T 32  RDN 33  X#0? 34  / 35  RCL 02           36  ENTER^ 37  X<> 06 38  - 39  RCL 01           40  STO L 41  X<> 05 42  ST- L 43  X<> L 44  R-P 45  X<>Y  46  ST+ T 47  RDN 48  * 49  P-R 50  ST+ 01           51  X<>Y 52  ST+ 02 53   E-8 54  LASTX 55  X>Y? 56  GTO 01  57  RCL 02 58  RCL 01 59  END

( 86 bytes / SIZE 007 )

 STACK INPUTS OUTPUTS T y0 / Z x0 / Y y1 y X x1 x

Example:   Find a root of      sinh z + z2 + 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   z0 = 0  ;  z1 = 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  zk+1 =  zk  or when   f(zk+1) =  f(zk)  Therefore, check f(z) in registers R03 & R04

b) Quadratic Interpolation

-Starting with 3 initial guesses:    z0 = x0+i.y0 ;  z1 = x1+i.y1 ;  z2 = x2+i.y2 , this program calculates

A = ( z1-z0 ).f(z2) + ( z0-z2 ).f(z1) + ( z2-z1 ).f(z0)
B = ( z1-z0 ).( 2.z2-z1-z0 ).f(z2) - ( z0-z2 )2.f(z1) + ( z2-z1 )2.f(z0)
C = ( z2-z1 ).( z1-z0 ).( z2-z0 ).f(z2)

and   z3 = z2 - 1 / [ B/(2.C) + E.( B2/(4C2) - A/C )1/2 ]      where   E = +1 or -1   ( the sign is choosen to produce the larger denominator )

-The process is repeated until  | zk+1-zk |  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

 01  LBL "RTZ2"   02  STO 01   03  RDN   04  STO 02   05  X<>Y   06  STO 05   07  STO 06   08  STO 09   09  STO 10   10  ST- T   11  -   12  R^   13  XEQ IND 00    14  STO 07   15  X<>Y   16  STO 08   17  RCL 02   18  RCL 01   19  RCL 05   20  ST+ X   21  ST- Z   22  -   23  XEQ IND 00   24  STO 11   25  X<>Y   26  STO 12   27  LBL 01   28  VIEW 01   29  RCL 02   30  RCL 01   31  XEQ IND 00   32  STO 03   33  X<>Y   34  STO 04   35  X<>Y 36  RCL 10   37  RCL 09   38  XEQ "Z*Z"   39  STO 13   40  X<>Y   41  STO 14   42  RCL 06   43  RCL 10   44  +   45  STO 20   46  RCL 05   47  RCL 09   48  +   49  STO 19   50  RCL 14   51  RCL 13   52  XEQ "Z*Z"    53  RCL 06   54  RCL 05   55  XEQ "Z*Z"       56  R-P   57  X=0?   58  GTO 02   59  CHS   60  STO 15   61  X<>Y   62  STO 16   63  RCL 06   64  RCL 20   65  +   66  RCL 05   67  RCL 19   68  +   69  RCL 14   70  RCL 13 71  XEQ "Z*Z"   72  STO 17   73  X<>Y   74  STO 18   75  RCL 20   76  RCL 19   77  RCL 08   78  RCL 07   79  XEQ "Z*Z"   80  ST- 13   81  X<>Y   82  ST- 14   83  X<>Y   84  RCL 20   85  RCL 19   86  XEQ "Z*Z"       87  ST- 17   88  X<>Y   89  ST- 18   90  RCL 06   91  RCL 05   92  RCL 12   93  RCL 11   94  XEQ "Z*Z"   95  ST+ 13   96  X<>Y   97  ST+ 14   98  X<>Y   99  RCL 06 100  RCL 05 101  XEQ "Z*Z" 102  ST+ 17 103  X<>Y 104  ST+ 18 105  RCL 14 106  RCL 13 107  R-P 108  STO 13 109  X<>Y 110  STO 14 111  RCL 18 112  RCL 17 113  R-P 114  RCL 16 115  ST- 14 116  ST- Z 117  X<> 15 118  ST/ 13 119  ST+ X 120  / 121  P-R 122  STO 15 123  X<>Y 124  STO 16 125  X<>Y 126  XEQ "Z^2" 127  RCL 14 128  RCL 13 129  P-R 130  XEQ "Z+Z" 131  XEQ "SQRTZ" 132  RCL 16 133  RCL Z 134  ST- 16 135  + 136  RCL 15 137  RCL Z 138  ST- 15 139  + 140  R-P 141  X<>Y 142  RCL 16 143  RCL 15 144  R-P 145  R^ 146  XY 153  P-R 154  ST+ 01 155  X<> 05 156  STO 09  157  X<>Y 158  ST+ 02 159  X<> 06 160  STO 10  161  RCL 03 162  X<> 07 163  STO 11 164  RCL 04 165  X<> 08           166  STO 12 167   E-8 168  LASTX 169  X>Y? 170  GTO 01  171  LBL 02 172  RCL 02 173  RCL 01 174  END

( 281 bytes / SIZE 021 )

 STACK INPUTS OUTPUTS Z h / Y y0 y X x0 x

( The 3 initial guesses are actually   z0 = x0+i.y0 ;   z0-(1+i).h ; z0-2.(1+i).h  )

Example:     With the same equation     sinh z + z2 + 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   z0 = 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  z0 = x0+i.y0  is needed but we have to compute the first derivative f '(z)

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

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

 01  LBL "NWTZ" 02  STO 01 03  X<>Y 04  STO 02 05  LBL 01 06  VIEW 01 07  RCL 02 08  RCL 01 09  XEQ IND 03 10  R-P 11  STO 04 12  X<>Y 13  STO 05 14  RCL 02 15  RCL 01 16  XEQ IND 00 17  R-P 18  ENTER^ 19  X<> 04 20  X#0? 21  / 22  X<>Y 23  RCL 05         24  - 25  X<>Y 26  P-R 27  ST- 01 28  X<>Y 29  ST- 02 30   E-8 31  LASTX         32  X>Y? 33  GTO 01  34  RCL 04         35  RCL 02 36  RCL 01 37  END

( 55 bytes / SIZE 006 )

 STACK INPUTS OUTPUTS Z / | f(z) | Y y0 y X x0 x

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

Example:       sinh z + z2 + 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  z0 = 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   z0 = x0+i.y0  , we now use the first and second derivatives   f '(z)  and  f "(z)  ,  the recursion formula is:
zk+1 = zk -  1 / [ f '(zk)/f(zk) + E.( ( f(zk)/(2.f '(zk)) )2 - f "(zk) /(2.f (zk)) )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

 01  LBL "RTZ4" 02  STO 01 03  X<>Y 04  STO 02 05  LBL 01 06  VIEW 01 07  RCL 02 08  RCL 01 09  XEQ IND 03 10  R-P 11  STO 06 12  X<>Y 13  STO 07 14  RCL 02 15  RCL 01 16  XEQ IND 04 17  R-P 18  STO 08 19  X<>Y 20  STO 09 21  RCL 02 22  RCL 01 23  XEQ IND 00 24  R-P 25  STO 05 26  X=0? 27  GTO 02 28  ST+ X 29  ST/ 06 30  ST/ 08 31  X<>Y 32  ST- 07 33  ST- 09 34  RCL 07 35  ST+ X 36  RCL 06 37  X^2 38  P-R 39  RCL 09 40  RCL 08         41  P-R 42  ST- Z 43  RDN 44  ST- Z 45  RDN 46  R-P 47  SQRT 48  X<>Y 49  2 50  / 51  X<>Y 52  P-R 53  RCL 07 54  RCL 06 55  P-R 56  STO 06 57  RDN 58  STO 07 59  RCL Z 60  ST- 07 61  + 62  RCL 06         63  RCL Z 64  ST- 06 65  + 66  R-P 67  X<>Y 68  RCL 07 69  RCL 06  70  R-P 71  R^ 72  XY 79  P-R 80  ST- 01 81  X<>Y 82  ST- 02 83   E-8 84  LASTX 85  X>Y? 86  GTO 01         87  LBL 02 88  RCL 05  89  RCL 02 90  RCL 01 91  END

( 123 bytes / SIZE 010 )

 STACK INPUTS OUTPUTS Z / | f(z) | Y y0 y X x0 x

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

Example:       sinh z + z2 + 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  z0 = 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.