# hp41programs

Polynomials for the HP-41

# Polynomials for the HP-41

Overview

1°) Real Polynomials

b) Cubic equations
c) Quartic equations
d) Evaluating a polynomial
e) First derivative of a polynomial
f)  Polynomial roots ( provided all roots are real )
g) Quadratic factors ( Bairstow's method )
h) Multiple roots
i)  Euclidean division
j)  Multiplication
m) Composition of 2 Polynomials

2°) Complex Polynomials

b) Cubic equations
c) Quartic equations
d) Evaluating a polynomial
e) First derivative of a polynomial
f) Second derivative of a polynomial
g) Polynomial roots
h) Laguerre's method
i) Multiple roots
j) Euclidean division
k) Multiplication

3°) Three short routines

a) Extremum of the curve y =  a.x2+b.x+c
b) Extrema of the curve y =  a.x3+b.x2+c.x+d
c) Center of symmetry of the curve y = a.x3+b.x2+c.x+d

Note:

-Some of these programs use synthetic registers M N O P
-They may be replaced by any standard registers but avoid any register usage conflict.

1°) Real Polynomials

-"P2" solves the equation:   a.x2+b.x+c = 0    with a # 0

Data Registers:  /
Flags:  F00  ( indicates complex roots )
Subroutines:  /

 01  LBL "P2"  02  CF 00  03  X<> Z  04  ST/ Z  05  ST+ X  06  /  07  CHS  08  ENTER^  09  X^2  10  RCL Z  11  -  12  X<0?  13  GTO 00  14  SQRT  15  RCL Y  16  SIGN  17  *  18  +  19  X#0?  20  ST/ Y  21  RTN  22  LBL 00  23  SF 00  24  CHS  25  SQRT  26  X<>Y  27  END

( 43 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Z a # 0 c/a Y b y X c x

-If  CF 00  the 2 solutions are  x , y
-If  SF 00   ------------------  x+i.y , x-i.y

Examples:     Solve  2.x2 + 3.x - 4 = 0   and   2.x2 + 3.x + 4 = 0

2  ENTER^
3  ENTER^
-4  XEQ "P2"  >>>>   -2.350781059
X<>Y                   0.850781060      Flag  F00 is clear, therefore the 2 solutions are  -2.350781059  and  0.850781060

2  ENTER^
3  ENTER^
4     R/S         >>>>   -0.75
X<>Y                    1.198957881      Flag  F00 is set, therefore the 2 solutions are  -0.75 + 1.198957881.i  and   -0.75 - 1.198957881.i

b) Cubic Equations

-"P3" finds the 3 roots of   a.x3+b.x2+c.x+d  by the Cardan's ( or Tartaglia's ) formulae:  ( with a # 0 )
-First, the term in x2 is removed by a change of argument, leading to  x3+p.x+q = 0
-Then, x = u+v  with u.v = -p/3 leads to a quadratic equation in u3

Data Registers:  /
Flags:  F01             ( indicates complex roots )
Subroutine:  none  if d # 0 , "P2" if d = 0

-Lines 04 to 22 are useful to produce exactly  x = 0 if  d = 0
-This is important when "P3" is called as a subroutine by "P4" or some "cosmological" programs

 01  LBL "P3"   02  DEG   03  CF 01   04  X#0?    05  GTO 00    06  RDN   07  XEQ "P2"   08  FS?C 00   09  SF 01   10  0   11  FS? 01   12  RTN   13  XY   15  RDN   16  XY   18  R^   19  XY   21  RTN   22  LBL 00   23  R^ 24  ST/ Z   25  ST/ T   26  /   27  R^   28  3   29  /   30  STO M      31  ST* Z   32  X^2   33  RDN   34  -   35  2   36  /   37  X<>Y   38  3   39  /   40  X<>Y   41  R^   42  ST- Z   43  RCL M   44  *   45  -   46  X<>Y 47  3   48  Y^X   49  RCL Y   50  X^2   51  +   52  X<=0?   53  GTO 01      54  SF 01   55  SQRT   56  RCL Y   57  SIGN   58  *   59  +   60  SIGN   61  LASTX   62  ABS   63  3   64  1/X   65  Y^X   66  *   67  ST/ Y   68  STO Z   69  X<>Y 70  ST- Z   71  +   72  60   73  SIN   74  *   75  RCL Y       76  2   77  /   78  CHS   79  R^   80  R^   81  GTO 02   82  LBL 01   83  CHS   84  SQRT   85  X<>Y   86  R-P   87  3   88  ST/ Z   89  1/X   90  Y^X   91  ST+ X   92  X<>Y 93  COS   94  LASTX       95  120   96  +   97  COS   98  R^   99  ST* Z 100  * 101  ENTER^ 102  CHS 103  RCL Z 104  ST- Y 105  RCL M 106  ST- T 107  LBL 02 108  CLX 109  X<> M 110  ST- Z 111  - 112  END

( 154 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS T a # 0 / Z b z Y c y X d x

-If  CF 01  the 3 solutions are  x ; y ; z
-If  SF 01   ------------------  x ; y+i.z ; y-i.z

Example:      Solve  2.x3-5.x2-7.x+1 = 0  and   2.x3-5.x2+7.x+1 = 0

2   ENTER^
-5  ENTER^
-7  ENTER^
1   XEQ "P3"  >>>>   3.467727065   RDN   0.131206041   RDN   -1.098933107   which are the 3 real solutions since flag F01 is clear.

2  ENTER^
-5  ENTER^
7   ENTER^
1     R/S         >>>>   -0.130131632   RDN   1.315065816   RDN   1.453569820

-Flag F01 is set , therefore the 3 solutions are  -0.130131633  ;  1.315065816 - 1.453569820.i  ;  1.315065816 + 1.453569820.i

Note:  The order of the 3 outputs x ; y ; z  is important when "P3" is used as a subroutine by "P4" hereafter and "WEF2" ( Weierstrass elliptic functions )

c) Quartic Equations

-"P4" solves the equation   x4+a.x3+b.x2+c.x+d = 0   ( if the leading coefficient is not 1 , divide all the equation by this coefficient ).
-First, the term in x3 is removed by a change of argument, leading to  x4+p.x2+q.x+r = 0  (E')
-Then, the resolvant  z3+p.z2/2+(p2-4r).z/16-q2/64 = 0 is solved by "P3" and if we call  z1 , z2 , z3  the 3 roots of this equation, the zeros of (E') are

x = z11/2 sign(-q) +/- ( z21/2+z31/2 )    ;    x = -(z11/2) sign(-q) +/- ( z21/2-z31/2 )

Data Registers:  R00 thru R04: temp
Flags:   F01 F02
Subroutine:    "P3"

-Lines 03 to 12 are only useful if  d = 0  to compute exactly the solution  x = 0.

 01  LBL "P4"   02  CF 02   03  X#0?     04  GTO 00   05  SIGN   06  RDN   07  XEQ "P3"   08  0   09  R^   10  R^   11  RTN   12  LBL 00   13  STO 03   14  RDN   15  STO 02   16  RDN   17  STO 04   18  X<>Y   19  CHS   20  STO 00   21  X^2   22  3   23  *   24  8   25  /   26  -   27  STO 01 28  RCL 00     29  4   30  ST/ 00   31  SQRT   32  ST/ 01   33  /   34  ENTER^   35  X^2   36  RCL 04   37  -   38  *   39  X<> 02   40  ST- 02   41  RCL 00   42  ENTER^   43  X^2   44  3   45  *   46  RCL 04   47  -   48  *   49  -   50  RCL 00   51  STO 04   52  *   53  ST+ 03   54  RCL 01 55  ENTER^   56  X^2   57  RCL 03   58  -   59  4   60  /   61  RCL 02   62  8   63  /   64  X^2   65  CHS   66  1   67  RDN   68  XEQ "P3"   69  SQRT   70  RCL 02   71  SIGN   72  *   73  ST+ 00   74  ST- 04   75  RDN   76  FS? 01   77  GTO 01   78  X<>Y   79  X<0?   80  GTO 02   81  SQRT 82  STO 03   83  X<>Y   84  SQRT   85  ST+ 03   86  -   87  ENTER^   88  CHS   89  RCL 04     90  ST+ Z   91  +   92  RCL 00   93  RCL 03   94  ST+ 00   95  -   96  RCL 00   97  RTN   98  LBL 01   99  R-P 100  SQRT 101  2 102  ST/ Z 103  * 104  P-R 105  ENTER^ 106  CHS 107  RCL 00 108  ST+ Z 109  + 110  R^ 111  RCL 04   112  RTN 113  LBL 02 114  CHS 115  SQRT 116  STO 02 117  X<>Y 118  CHS 119  SQRT 120  STO 03 121  ST+ 02 122  - 123  X#0? 124  SF 02 125  X=0? 126  RCL 00 127  RCL 00 128  SF 01 129  RCL 02 130  RCL 04 131  END

( 164 bytes / SIZE 005 )

 STACK INPUTS OUTPUTS T a t Z b z Y c y X d x

-If  CF 01  &  CF 02  the 4 solutions are    x ; y ; z ; t
-If  SF 01  &  CF 02   ------------------    x+i.y ; x-i.y ; z ; t
-If  SF 01  &  SF 02   ------------------    x+i.y ; x-i.y ; z+i.t ; z-i.t

Example1:    Solve  x4-2.x3-35.x2+36.x+180 = 0

-2   ENTER^
-35  ENTER^
36   ENTER^
180  XEQ "P4"  >>>>   6  RDN   3   RDN   -2   RDN   -5    Flags F01 & F02 are clear , the 4 solutions are   6 ; 3 ; -2 ; -5

Example2:    Solve  x4-5.x3+11.x2-189.x+522 = 0

-5   ENTER^
11   ENTER^
-189   ENTER^
522      R/S     >>>>   -2  RDN    5   RDN   3   RDN   6  Flag F01 is set & F02 is clear , the 4 solutions are  -2+5.i ; -2-5.i ; 3 ; 6

Example3:    Solve  x4-8.x3+26.x2-168.x+1305 = 0

-8   ENTER^
26   ENTER^
-168  ENTER^
1305    R/S      >>>>   -2  RDN    5   RDN   6   RDN   3    Flags F01 & F02 are set , the 4 solutions are  -2+5.i ; -2-5.i ; 6+3.i ; 6-3.i

N.B:   The method used in this version is less simple but much more reliable than the previous program.

d) Evaluating a real Polynomial

-The following program calculates  p(x) = a0.xn+a1.xn-1+ ... + an-1.x+an  for any given x-value.
-The coefficients   a0 , a1 ,  ...... , an-1 , an  are to be stored into contiguous registers.

Data Registers:      •  Rbb = a0  •  Rbb+1 = a1 , ....... , •  Ree = an         ( These n+1 registers are to be initialized before executing "PL" )
Flags: /
Subroutines:  /

 01  LBL "PL"  02  0  03  LBL 01  04  RCL Y  05  *  06  RCL IND Z  07  +  08  ISG Z  09  GTO 01  10  X<>Y  11  SIGN  12  RDN  13  END

( 24 bytes )

 STACK INPUTS OUTPUTS Y bbb.eee / X x p(x) L / x

Example:     p(x) = 2.x3+3.x2-4.x+7     Find  p(5)

-If we store these 4 coefficients into R01 to R04  ( 2  STO 01  3  STO 02  -4  STO 03  7  STO 04 )

1.004  ENTER^
5     XEQ "PL"  >>>>  p(5) = 312

Note:   Lines 10-11-12  are only usefull to save x in L-register.

e) First derivative of a Polynomial

-"dPL" calculates  p'(x) = n.a0.xn-1+(n-1).a1.xn-2+ ... + an-1  for any given x-value.
-The coefficients   a0 , a1 ,  ...... , an-1 , an  are to be stored into contiguous registers.

Data Registers:   •  Rbb = a0 , •  Rbb+1 = a1 , ....... , •  Ree = an                        ( These n+1 registers are to be initialized before executing "dPL" )
In fact, Ree is unused.
Flags: /
Subroutines:  /

 01  LBL "dPL" 02  RCL Y 03  FRC 04   E3 05  * 06  R^ 07  INT 08  - 09  STO M 10  CLX 11  LBL 01 12  RCL Y 13  * 14  RCL IND Z 15  X<> M 16  ST* M 17  X<> M 18  + 19  ISG Z 20  ""        21  DSE M 22  GTO 01 23  X<>Y 24  SIGN 25  RDN 26  END

( 45 bytes )

 STACK INPUTS OUTPUTS Y bbb.eee / X x p'(x) L / x

Example:     p(x) = 2.x3+3.x2-4.x+7     Find  p'(5)

-If we store these 4 coefficients into R01 to R04  ( 2  STO 01  3  STO 02  -4  STO 03  7  STO 04 )

1.004  ENTER^
5     XEQ "dPL"  >>>>  p'(5) = 176

f) Real Polynomial Roots ( assuming all roots are real )

-The following program solves the equation   p(x) = a0.xn+a1.xn-1+ ... + an-1.x+an = 0  provided all roots are real.
-The coefficients   a0 , a1 ,  ...... , an-1 , an  are to be stored into contiguous registers ( Rbb thru Ree ).

-Starting with an initial guess  x0 , the Newton's formula   xk+1 = xk-p(xk)/p'(xk) is applied until  | p(xk)/p'(xk) | is smaller than 10-9  ( line 25 )
-Then, p is replaced by  p/(x-r)  (  lines 33 to 40 ) and the root is stored into Ree
-The process is repeated until all roots are found ( polynomial deflation ).
-Finally, the roots are in registers  Rbb+1 , ..... , Ree.

Data Registers:  R00 to R04: temp                     ( Registers Rbb thru Ree are to be initialized before executing "PLR" )

•  Rbb = a0 , •  Rbb+1 = a1 , ....... , •  Ree = an       bb > 04
Flags: /
Subroutines:   "PL"  "dPL"

 01  LBL "PLR"   02  STO 01 03  X<>Y 04  STO 00 05  STO 03 06  STO 04 07  ISG 04 08  LBL 01 09  VIEW 01 10  RCL 03 11  RCL 01 12  XEQ "dPL"   13  STO 02 14  RCL 03 15  RCL 01 16  XEQ "PL" 17  RCL 02 18  / 19  ST- 01 20  RCL 01 21  X=0? 22  SIGN 23  / 24  ABS 25   E-9 26  X

( 79 bytes )

 STACK INPUTS OUTPUTS Y bbb.eee / X x0 1+bbb.eee

where  bbb.eee  is the control number of the polynomial ( bbb > 004 )  and  x0 is an initial guess.

Example:   Find all the roots of   p(x) = 2.x5+3.x4-35.x3-10.x2+128.x-74

For instance   2  STO 05  3  STO 06  -35  STO 07  -10  STO 08  128 STO 09  -74 STO 10  ( the control number is 5.010 )  and if we choose  x0 = 1

5.010  ENTER^
1     XEQ "PLR"  the successive approximations are displayed and 1mn55s  later,
we get  6.010 = the control number of the solutions ( in R06 R07 R08 R09 R10 )

RCL 06  >>>>  -4.373739462
RCL 07  >>>>  -2.455070118
RCL 08  >>>>   2.984066207
RCL 09  >>>>   1.641131729
RCL 10  >>>>   0.703611645

N.B:     If you get "DATA ERROR" at line 18 ( meaning p'(x) = 0 ) , change register R01 and  XEQ 01

-"BRST" factorizes the polynomial   p(x) = a0.xn+a1.xn-1+ ... + an-1.x+an  into quadratic factors and solves  p(x) = 0      ( n > 1 )
-If deg(p) is odd, we have   p(x) =    (a0.x+b).(x2+u1.x+v1)...............(x2+um.x+vm)      where m = (n-1)/2
-If deg(p) is even                p(x) = (a0x2+u1.x+v1)(x2+u2.x+v2)..........(x2+um.x+vm)      where m  =  n/2

-The coefficients u and v are found by the Newton method for solving 2 simultaneous equations.
-Then  p is divided by (x2+u.x+v)  and  u & v  are stored into Ree-1 & Ree respectively
-The process is repeated until all quadratic factors are found, and when the program stops ( line 89 ):

If deg(p) is odd          Rbb = a0 , Rbb+1 = b , Rbb+2 = u1 , Rbb+3 = v1 , ........ , Ree-1 = um , Ree = vm
If deg(p) is even         Rbb = a0 , Rbb+1 = u1 , Rbb+2 = v1 , ............................ , Ree-1 = um , Ree = vm

-Then, press R/S to get the successive roots.

Data Registers:  R00 to R08: temp.   ( R06 = u ; R07 = v ; R08 = bbb.eee )

•  Rbb = a0  •  Rbb+1 = a1 , ....... , •  Ree = an       ( These n+1 registers are to be initialized before executing "BRST")     bb > 08
Flags:  F00
Subroutines:  "DIV"  ( see below )   and "P2"    ( only to find the roots )

-Lines 67-68 may also be replaced by RND  X#0?  and the accuracy will be controlled by the display setting.

 01  LBL "BRST"   02  STO 06   03  RDN   04  STO 07   05  X<>Y   06  STO 08   07  STO 04   08  LBL 01   09  VIEW 06   10  CLX   11  STO 00   12  STO 01   13  STO 02   14  STO 03   15  RCL 04   16  STO 05   17  LBL 02   18  RCL IND 05   19  RCL 00   20  RCL 07   21  *   22  -   23  RCL 01   24  RCL 06   25  *   26  -   27  X<> 01 28  STO 00   29  RCL 02   30  RCL 07   31  *   32  -   33  RCL 03   34  RCL 06   35  *   36  -   37  X<> 03   38  X<> 02   39  ISG 05   40  GTO 02   41  STO 05   42  RCL 01   43  *   44  RCL 00          45  RCL 02   46  ST* 01   47  *   48  -   49  RCL 03   50  RCL 05   51  *   52  RCL 02   53  X^2   54  - 55  STO 05   56  /   57  ST+ 06   58  RCL 00   59  RCL 03   60  *   61  RCL 01   62  -   63  RCL 05   64  /   65  ST+ 07   66  R-P   67   E-8      68  XY 112  / 113  CHS 114  SF 00 115  TONE 9 116  STOP 117  ISG 04 118  LBL 04 119  1 120  LBL 05 121  RCL IND 04 122  ISG 04 123  RCL IND 04 124  XEQ "P2" 125  ISG 04 126  FS? 30 127  GTO 06 128  TONE 9 129  STOP 130  GTO 04 131  LBL 06 132  BEEP 133  END

( 190 bytes )

 STACK INPUTS OUTPUT0 OUTPUTS1 ......................... OUTPUTSk Z bbb.eee / / / Y v0 / y1 yk X u0 bbb.eee x1 xk

where  bbb.eee  is the control number of the polynomial ( bbb > 008 )  and  u0 and  v0  are initial guesses.

-When the program stops ( line 89 )  Rbb thru Ree contain the factorization
-Then, the successive outputs are the different pairs of roots:     xk  ,  yk   if flag F00 is clear  or   xk+i.yk  ,   xk-i.yk  if F00 is set

However , when deg(p) is odd , the first root is always real:  F00 is set but y = 0

-The last pair of roots is indicated by a BEEP, the others by a TONE 9.

Example1:     Find all the roots of   2.x5+7.x4+20.x3+81.x2+190.x+150

-For instance,  2 STO 11   7 STO 12   20 STO 13   81 STO 14   190 STO 15   150 STO 16   ( control number = 11.016 ) and if we choose  u0 = v0 = 0

11.016  ENTER^
0         ENTER^
XEQ "BRST"  the successive uk are displayed and 2 minutes later, we hear a BEEP and X = 11.016

R11 = 2   R12 = 3     //     R13 = -2   R14 = 10     //    R15 = 4   R16 = 5            whence   p(x) = (2x+3).(x2-2x+10).(x2+4x+5)    then:

R/S  >>> ( TONE 9 )  -1.5    X<>Y   0   with  F00 set     the first root is   -1.5    ( deg(p) is odd , the 1st root is real )
R/S  >>> ( TONE 9 )    1      X<>Y   3    with  F00 set     1+3.i  and  1-3.i
R/S  >>>   ( BEEP )     -2      X<>Y   1    with  F00 set    -2+i  and  -2-i            the 5 roots are  -1.5 ; 1+3.i ; 1-3.i ; -2+i ; -2-i

Example2:     Solve      x6-6.x5+8.x4+64.x3-345.x2+590.x-312 = 0

-For instance,    1 STO 11   -6 STO 12   8 STO 13   64 STO 14   -345 STO 15   590 STO 16   -312 STO 17  ( bbb.eee = 11.017 ) and with  u0 = v0 = 0

11.017  ENTER^
0         ENTER^
XEQ "BRST"  the successive uk are displayed and 2 minutes later, we hear a BEEP and X = 11.017

R11 = 1   R12 = 1   R13 = -12    //     R14 = -4   R15 = 13     //     R16 = -3   R17 = 2            whence   p(x) = (x2+x-12).(x2-4x+13).(x2-3x+2)   then:

R/S  >>>  ( TONE 9 )  -4      X<>Y   3   with  F00 clear  the first pair of roots are -4 and  3
R/S  >>>  ( TONE 9 )   2      X<>Y   3    with  F00 set        2+3.i  and  2-3.i
R/S  >>>    ( BEEP )     2      X<>Y   1    with  F00 clear        2     and    1                 the 6 roots are   -4 ; 3 ; 2+3.i ; 2-3.i ; 2 ; 1

Notes:

-If you get "DATA ERROR"  ( line 56 ) or if the process seems to diverge, stop the program, change registers R06 & R07 and press  XEQ 01
-If you want to see the roots again, press XEQ 10
-If you need the factorization only, lines 89 to 132 may be deleted.

h) Multiple Roots

-"PLR"  may be disappointing to find multiple roots: slow convergence and low accuracy are to be expected.
-"MSR" changes multiple roots into single roots.
-More exactly, if p(x) is a polynomial,  s = p/GCD(p;p') and p have the same distinct roots but s has single roots only. ( GCD = Greatest Common Divisor )
and the following program calculates the coefficients of s(x) using the Euclidean algorithm.
-Then, "PLR" or "BRST" can be applied to s(x).

Data Registers:  R00 to R06: temp   when the program stops, R02 = R05 = the control number of  GCD(p;p')

•  Rbb = a0  •  Rbb+1 = a1 , ....... , •  Ree = an       ( These n+1 registers are to be initialized before executing "MSR" )  bb > 06

Registers  R(ee+1) ........R(bb+3n-1) are also used
Flags: /
Subroutines:  "DIV"  ( see below )

 01  LBL "MSR"  02  ENTER^ 03  STO 04 04  FRC 05   E3 06  * 07  RCL 04 08  INT 09  - 10  STO 01 11  1 12  + 13  .1 14  % 15  + 16  + 17  STO 05 18  LASTX 19  + 20   E-3 21  - 22  STO 06 23  RCL 04 24  RCL 05 25  LBL 00 26  RCL IND Y 27  STO IND Y 28  ISG Y 29  RDN 30  ISG Y 31  GTO 00 32  RCL 04 33  RCL 06 34  LBL 01 35  RCL IND Y 36  RCL 01 37  * 38  STO IND Y 39  ISG Z 40  ISG Y 41  RDN 42  DSE 01 43  GTO 01 44  LBL 02 45  RCL 05 46  RCL 06 47  XEQ "DIV" 48  SIGN 49  - 50  STO 05 51  LBL 03 52  ISG 05 53  X<0? 54  GTO 04 55  RCL IND 05 56  ABS 57   E-4 58  X>Y? 59  GTO 03 60  LBL 04 61  RCL 06 62  X<> 05 63  STO 06 64  1 65  - 66  ISG X 67  GTO 02 68  RCL 05 69  RCL IND 05 70  LBL 05 71  ST/ IND Y 72  ISG Y 73  GTO 05 74  RCL 04 75  RCL 05 76  XEQ "DIV" 77  END

( 121 bytes / SIZE>3n+8 )

 STACK INPUTS OUTPUTS X (bbb.eee)p (bbb.eee)s

( bbb > 006 )

Example1:    Find all the roots of the polynomial  p(x) = x4+2.x3+5.x2+4.x+4

-If we choose  bbb = 007    1  STO 07   2  STO 08   5  STO 09   4  STO 10  STO 11  ( control number = 7.011 )

7.011  XEQ "MSR"  >>>>  7.009  ( in 19 seconds )

RCL 07 = 1
RCL 08 = 1
RCL 09 = 2   whence  s(x) = x2+x+2      ( in fact,  p(x) = ( x2+x+2 )2 )

1  ENTER^
1  ENTER^
2  XEQ "P2"  >>>   -0.5   X<>Y  1.322875656  with flag F00 set.  Therefore the 2 solutions are    -0.5 - 1.322875656.i  and  -0.5 + 1.322875656.i

Example2:    Find all the roots of  p(x) = x10-17.x9+127.x8-549.x7+1521.x6-2823.x5+3557.x4-3007.x3+1634.x2-516.x+72

-Let's store these 11 coefficients 1 ; -17 ; 127 ; .... ; -516 ; 72   into   R07 ; R08 ; R09 ;        ; R16 ; R17 respectively  ( control number = 7.017 )

7.017  XEQ "MSR"  >>>> 7.010  ( in 49 seconds )    R07 = 1 ;  R08 = -6 ; R09 = 11 ; R10 = -6  whence  s(x) = x3-6.x2+11.x-6
( All the coefficients of p(x) are integers, so we can be sure all the coefficients of s(x) are integers too )

1   ENTER^
-6  ENTER^
11  ENTER^
-6  XEQ "P3"  >>>  3  RDN  2  RDN  1 with flag F01 clear.  Therefore p(x) = 0 has 3 solutions:  1 ; 2 ; 3   ( Actually p(x) = (x-1)5(x-2)3(x-3)2  )

-we have also:       R03 = 31.038   R31 = 1 ; R32 = -11 ; R33= 50 ; R34 = -122 ; R35 = 173 ; R36 = -143 ; R37 = 64 ; R38 = -12
whence          GCD(p;p') =  x7-11.x6+50.x5-122.x4+173.x3-143.x2+64.x-12

Notes:

-"MSR" may be applied to GCD(p;p') to find the multiplicities of the distinct roots.
-The Euclidean algorithm stops when the remainder equals zero.
-In this program, the leading coefficient of the remainder is deleted if it is smaller than 10-4 ( line 57 ).
-Another choice is sometimes better. For instance,  s(x) will be wrongly computed if p(x) = x3+0.0001.x2
but if you replace line 57 by  E-9 the good result  s(x) = x2+0.0001.x  will be obtained!

Alternative:   Replace lines 56-57-58  with  RND  X=0?  and execute "MSR" several times  in different  FIX  formats

i) Euclidean Division

-If   a(x) = a0.xn+a1.xn-1+ ... + an-1.x+an  and   b(x) = b0.xm+b1.xm-1+ ... + bm-1.x+bm  are 2 polynomials,
there are only 2 polynomials q(x) and r(x)  such that   a = b.q + r  with  deg(r) < deg(b)

-The following program calculates q and r.

Data Registers:  R00 to R03: temp                       ( Registers Rbb thru Ree & Rbb' thru Ree' are to be initialized before executing "DIV" )

•  Rbb = a0   •  Rbb+1 = a1  , ....... , •  Ree = an      bb > 03
•  Rbb' = b0  •  Rbb'+1 = b1 , ....... , •  Ree' = bm       bb' > 03
Flags: /
Subroutines: /

 01  LBL "DIV"   02  STO 02 03  FRC 04  X<>Y 05  STO 01 06  FRC 07  X<>Y 08  - 09   E3 10  * 11  RCL 01 12  INT 13  STO 03 14  - 15  RCL 02 16  INT 17  + 18  STO 00 19  ISG 00 20  LBL 01 21  RCL 01 22  RCL 02 23  RCL IND Y 24  RCL IND Y 25  / 26  STO IND Z 27  ISG Z 28  SIGN 29  ISG Y 30  X=0? 31  GTO 03 32  LBL 02 33  CLX 34  RCL IND Y 35  LASTX 36  * 37  ST- IND Z  38  ISG Z 39  CLX 40  ISG Y 41  GTO 02 42  LBL 03 43  ISG 01 44  CLX 45  DSE 00 46  GTO 01 47  RCL 01 48  RCL 03 49  RCL 01       50  INT 51  1 52  - 53   E3 54  / 55  + 56  END

( 81 bytes )

 STACK INPUTS OUTPUTS Y (bbb.eee)a (bbb.eee)r X (bbb.eee)b (bbb.eee)q

where         (bbb.eee)a = the control number of the dividend                        (bbb.eee)r  = the control number of the remainder        ( all  bbb > 003 )
(bbb.eee)b = the control number of the divisor                           (bbb.eee)q = the control number of the quotient

Example:     a(x) =  2.x5+5.x4-21.x3+23.x2+3.x+5      b(x) = 2.x2-3.x+1     Find   q(x) and r(x)

For instance:        2 STO 04   5 STO 05   -21 STO 06   23 STO 07   3 STO 08   5 STO 09     ( control number = 4.009 )
and        2  STO 11   -3 STO 12   1 STO 13                                                                ( control number = 11.013 )

4.009  ENTER^
11.013  XEQ "DIV"  >>>> ( in 6 seconds )    4.007  X<>Y  8.009        ( q and r have replaced a  ;  b is unchanged )

R04 = 1   R05 = 4   R06 = -5   R07 = 2   whence  q(x) = x3+4.x2-5.x+2
R08 = 14   R09 = 3                   whence  r(x) = 14.x + 3

Notes:

-When b(x) is constant, (bbb.eee)r = 1+eee.eee
-Roundoff errors may produce small leading coefficients instead of zero in the remainder.

"DIV" doesn't work if  deg(a) < deg(b)    but in this case,  q = 0 and r = a

j) Polynomial Multiplication

-Let  a(x) = a0.xn+a1.xn-1+ ... + an-1.x+an  and   b(x) = b0.xm+b1.xm-1+ ... + bm-1.x+bm   2 polynomials,

"PRO" computes and stores the coefficients of  p(x) = a(x).b(x) into contiguous registers. ( you have to specify the first of these registers )

Data Registers:  R00 to R03: temp                       ( Registers Rbb thru Ree & Rbb' thru Ree' are to be initialized before executing "PRO" )

•  Rbb = a0   •  Rbb+1 = a1  , ....... , •  Ree = an      bb > 03
•  Rbb' = b0  •  Rbb'+1 = b1 , ....... , •  Ree' = bm       bb' > 03
Flags: /
Subroutines: /

-If you don't have an HP-41CX replace line 23  by     0  LBL 00  STO IND Y  ISG Y  GTO 00

 01  LBL "PRO"   02  ENTER^ 03  R^ 04  STO 01 05  FRC 06  R^ 07  STO 02 08  FRC 09  + 10  X<>Y 11  RCL 01        12  INT 13  - 14  RCL 02 15  INT 16  - 17   E3 18  / 19  + 20  + 21  STO 00        22  STO 03 23  CLRGX   24  LBL 01 25  RCL 01 26  RCL 00 27  LBL 02 28  RCL IND Y 29  RCL IND 02 30  * 31  ST+ IND Y 32  ISG Y 33  RDN 34  ISG Y 35  GTO 02 36  ISG 00 37  CLX 38  ISG 02 39  GTO 01        40  RCL 03 41  END

( 60 bytes )

 STACK INPUTS OUTPUTS Z (bbb.eee)a Y (bbb.eee)b X bbbp (bbb.eee)p

where         (bbb.eee)a = the control number of  a(x)
(bbb.eee)b = the control number of  b(x)                           (bbb.eee)p = the control number of the product        ( all  bbb > 003 )

Example:     a(x) =  2.x5+5.x4-21.x3+23.x2+3.x+5      b(x) = 2.x2-3.x+1         Find  p(x) = a(x).b(x)

For instance:        2 STO 04   5 STO 05   -21 STO 06   23 STO 07   3 STO 08   5 STO 09     ( control number = 4.009 )
2  STO 11   -3 STO 12   1 STO 13                                                                ( control number = 11.013 )

and if we want to have p(x) in registers R21  R22  ... etc ...

4.009  ENTER^
11.013  ENTER^
21     XEQ "PRO"  yields  21.028   ( in 7 seconds )

-We have   R21 = 4   R22 = 4   R23 = -55   R24 = 114   R25 = -84   R26 = 24   R27 = -12   R28 = 5

whence  p(x) = 4.x7+4.x6-55.x5+114.x4-84.x3+24.x2-12.x+5

Note:    Neither (bbb.eee)a  &   (bbb.eee)p   nor   (bbb.eee)b  &   (bbb.eee)p   can overlap.

-Let  a(x) = a0.xn+a1.xn-1+ ... + an-1.x+an  and   b(x) = b0.xm+b1.xm-1+ ... + bm-1.x+bm   2 polynomials,

"ADD" computes and stores the coefficients

of  p(x) = a(x) + b(x)  if F01 is clear
or  p(x) = a(x) - b(x)  if F01 is set          into contiguous registers. ( you have to specify the first of these registers )

Data Registers:  R00 to R03: temp                       ( Registers Rbb thru Ree & Rbb' thru Ree' are to be initialized before executing "ADD" )

•  Rbb = a0   •  Rbb+1 = a1  , ....... , •  Ree = an      bb > 03
•  Rbb' = b0  •  Rbb'+1 = b1 , ....... , •  Ree' = bm       bb' > 03

Flags:    F01      -If flag F01 is clear, "ADD" calculates the sum of 2 polynomials
-If flag F01 is set,  "ADD" calculates the difference of 2 polynomials
Subroutines:     "DEL"

-If you don't have an HP-41CX replace line 28  by     0  LBL 00  STO IND Y  ISG Y  GTO 00   RCL 00

 01  LBL "ADD"   02  STO 00 03  RDN 04  STO 02 05  FRC 06   E3 07  * 08  RCL 02 09  INT 10  - 11  X<>Y 12  STO 01 13  FRC 14   E3 15  * 16  RCL 01         17  INT 18  - 19  XY 21  RCL 00 22  + 23   E3 24  / 25  RCL 00 26  + 27  STO 00 28  CLRGX 29  XEQ 01 30  LASTX 31  - 32  STO 03 33  RCL 01         34  XEQ 01 35  STO 01 36  RCL 02 37  XEQ 01 38  STO 02 39  GTO 02 40  LBL 01 41  INT 42  DSE X 43  LASTX 44  FRC 45   E3 46  ST/ Z 47  * 48  + 49  1 50  + 51  RTN 52  LBL 02 53  CLX 54  DSE 01 55  RCL IND 01 56  ST+ IND 03 57  CLX 58  DSE 02 59  RCL IND 02 60  FS? 01 61  CHS 62  ST+ IND 03 63  DSE 03 64  GTO 02 65  RCL 00 66  XEQ "DEL" 67  END

( 102 bytes )

 STACK INPUTS OUTPUTS Z (bbb.eee)a Y (bbb.eee)b X bbbp (bbb.eee)p

where         (bbb.eee)a = the control number of  a(x)
(bbb.eee)b = the control number of  b(x)                   (bbb.eee)p = the control number of the sum ( if CF 01 ) or the difference ( if SF 01 )

all  bbb's > 003

Example:

a(x) = 2.x3+4.x2+5.x+6     b(x) = 2.x3-3.x2+7.x+1    Find  a+b  and  a-b

2 STO 06   4 STO 07   5 STO 08   6 STO 09   ( control number = 6.009 )
2 STO 11   -3 STO 12   7 STO 13   1 STO 14  ( control number = 11.014 )   and if we want to have p(x) in registers R21  R22  ... etc ...

CF 01
6.009  ENTER^
11.014  ENTER^
21     XEQ "ADD"  gives  21.024   ( in 5 seconds )   R21 = 4  R22 = 1  R23 = 12  R24 = 7  whence  a(x)+b(x) = 4.x3+x2+12.x+7

SF 01
6.009  ENTER^
11.014  ENTER^
21     XEQ "ADD"  gives  22.024   ( in 5 seconds )   R22 = 7  R23 = -2  R24 = 5  whence  a(x)-b(x) = 7.x2-2.x+5

( without line 66 , we would get  21.024  with the superfluous  R21 = 0 )

-This short routine is useful to delete tiny dominant coefficients resulting from roundoff errors or occuring in additions or subtractions.

Data Registers:                           ( Registers Rbb thru Ree are to be initialized before executing "DEL" )

•  Rbb = a0   •  Rbb+1 = a1  , ....... , •  Ree = an
Flags: /
Subroutines: /

 01  LBL "DEL" 02  LBL 01 03  RCL IND X 04  ABS 05   E-7      06  X Z 09  ISG X 10  GTO 01 11  1 12  ST- Y 13  0 14  STO IND Z 15  LBL 02 16  X<> Z 17  END

( 35 bytes )

 STACK INPUTS OUTPUTS X (bbb.eee)p (bbb'.eee)p

Example:      If   p(x) = 10-9x2+3.x-7    and   R01 = 10-9   R02 = 3  R03 = -7     ( control number = 1.003 )

1.003   XEQ "DEL"  yields   2.003   meaning  p(x) = 3.x-7

Note:

-If p(x) = 10-9 is stored in R01     1.001  XEQ "DEL" gives  1.001  but 0 is stored into register R01

m) Composition of 2 Polynomials

-Let  f(x) = a0.xn+a1.xn-1+ ... + an-1.x+an  and   g(x) = b0.xm+b1.xm-1+ ... + bm-1.x+bm   2 polynomials,

"PCOMP" computes and stores the coefficients of  p(x) = f[g(x)]  into contiguous registers. ( you have to specify the first of these registers )

Data Registers:  R00 to R03: temp                       ( Registers Rbb thru Ree & Rbb' thru Ree' are to be initialized before executing "COMP" )

•  Rbb = a0   •  Rbb+1 = a1  , ....... , •  Ree = an      bb > 08
•  Rbb' = b0  •  Rbb'+1 = b1 , ....... , •  Ree' = bm       bb' > 08        REE+1 .... REE+n.m+1 are also used for temporary data storage.
Flags: /
Subroutines:   "PRO"  ( cf § 1°) j) )  &  "LCO" ( cf "Miscellaneous Short Routines for the HP-41" )

 01  LBL "PCOMP" 02  STO 07 03  .1 04  % 05  + 06  STO 06 07  RDN 08  STO 05 09  INT 10  LASTX 11  FRC 12   E3 13  * 14  - 15  X<>Y 16  STO 04         17  STO 08 18  INT 19  LASTX 20  FRC 21   E3 22  * 23  - 24  * 25  ST+ 07 26  RCL IND 04 27  STO IND 06 28  ISG 07 29  LBL 01 30  RCL 06 31  ISG 08 32  X=0? 33  RTN 34  RCL 05 35  RCL 07 36  XEQ "PRO"  37  RCL 06 38  INT 39  XEQ "LCO" 40  STO 06 41  FRC 42   E3 43  * 44  RCL IND 08 45  ST+ IND Y 46  GTO 01 47  END

( 77 bytes )

 STACK INPUTS OUTPUTS Z bbb.eee(f) / Y bbb.eee(g) / X BBB BBB.EEE(fog)

Example:        f(x) = 2 x2 + 3 x + 7    g(x) = 4 x3 + 5 x2 + 6 x + 1

-Store for instance   [ 2  3  7 ]  into  R11  R12  R13          control number = 11.013
and   [ 4  5  6  1 ] into  R14  R15  R16  R17  control number = 14.017

-If you choose  BBB = 18

11.013  ENTER^
14.017  ENTER^
18     XEQ "PRO"  >>>>   18.024         --- Execution time = 16 seconds ---

R18 = 32   R19 = 80   R20 = 146   R21 = 148   R22 = 107   R23 = 42   R24 = 12

-So  fog(x) = 32 x6 + 80 x5 + 146 x4 + 148 x3 + 107 x2 + 42 x + 12

-Likewise, you'll find in 21 seconds,  gof(x) = 32 x6 + 144 x5 + 572 x4 + 1176 x3 + 2129 x2 + 1992 x + 1660

2°) Complex Polynomials

-Most of the above methods are now applied to complex polynomials:

-"P2Z"  solves   z2+(a+ib).z+(c+id) = 0

Data Registers:  /
Flags: /
Subroutines:  /

 01  LBL "P2Z"  02  X<> Z 03  CHS 04  STO N 05  R^ 06  CHS 07  STO M 08  R-P 09  X^2 10  X<>Y 11  ST+ X 12  X<>Y 13  P-R 14  R^ 15  ST+ X 16  ST+ X     17  ST- Z 18  X<> T 19  4 20  * 21  - 22  R-P 23  4 24  ST/ Y 25  SQRT     26  ST/ M 27  ST/ N 28  ST/ Z 29  1/X 30  Y^X 31  P-R 32  RCL N     33  X<> Z 34  ST- Z 35  ST+ N 36  X<> N 37  RCL M     38  X<> Z 39  ST- Z 40  ST+ M 41  X<> M 42  CLA 43  END

( 73 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS T a y' Z b x' Y c y X d x

-The solutions are:   x+i.y ; x'+i.y'

Example:   Solve  z2-(6+10.i).z + (-13+26.i) = 0

-6   ENTER^
-10  ENTER^
-13  ENTER^
26  XEQ "P2Z"  >>>>   4   RDN   7   RDN   2   RDN   3    whence   z1 = 4 + 7.i   and   z2 = 2 + 3.i      ( in 4 seconds )

b) Complex cubic equation

-This program solves the equation:    (a+i.b).z3 + (c+i.d).z2 + (e+i.f).z + (g+i.h) = 0     assuming (a;b) # (0;0)

Data Registers:            R00 & R09: temp                        ( Registers R01 thru R08 are to be initialized before executing "P3Z" )

•  R01 = a    •  R03 = c   •  R05 = e   •  R07 = g
•  R02 = b    •  R04 = d   •  R06 = f   •  R08 = h

-When the program stops, the 3 solutions are in registers R01 thru R06

Flags: /
Subroutines:  "Z/Z"  "Z*Z"  "Z-Z"  "Z^2"  "SHZ"  "ASHZ"  "SQRZ"  "Z^X"  ( cf "Complex Functions for the HP-41" )

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

( 282 bytes / SIZE 010 )

 STACK INPUTS OUTPUTS T / D Z / C Y / B X / A

A+i.B & C+i.D  are 2 complex roots of the cubic.

Example:    Find the 3 complex roots of      (2+3.i).z3 + (4+5.i).z2 + (6+7.i).z + (8+9.i) = 0

-Store   2 , 3 , 4 , 5 , 6 , 7 , 8 , 9  into registers  R01 to R08

and  XEQ "P3Z"  >>>>   -0.179640712    RDN    -1.436921967   whence   z1 = -0.179640712 -1.436921967.i  =  R01+i.R02    ( in 35 seconds )
RDN   -0.061935981    RDN      1.506099080   whence   z2 = -0.061935981 + 1.506099080.i = R03+i.R04
[ RCL 05 ] -1.527654077 [ RCL 06 ] 0.084669040  whence   z3 =  -1.527654077 + 0.084669040.i = R05+i.R06

c) Complex quartic equation

-This program solves the equation:    (a+i.b).z4 + (c+i.d).z3 + (e+i.f).z2 + (g+i.h).z + (k+i.l) = 0     assuming (a;b) # (0;0)

Data Registers:            R00 & R11 to R15: temp                        ( Registers R01 thru R10 are to be initialized before executing "P4Z" )

•  R01 = a    •  R03 = c   •  R05 = e   •  R07 = g   •  R09 = k
•  R02 = b    •  R04 = d   •  R06 = f   •  R08 = h    •  R10 = l

-When the program stops, the 4 solutions are in registers R01 thru R08

Flags: /
Subroutines:  "P2Z"  "P3Z" &  "Z/Z"  "Z*Z"  "Z-Z"  "Z^2"   "SQRZ"   ( cf "Complex Functions for the HP-41" )

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

( 324 bytes / SIZE 016 )

 STACK INPUTS OUTPUTS T / D Z / C Y / B X / A

where  A+i.B & C+i.D  are 2 complex roots of the quartic.

Example:    Find the 4 complex roots of      (2+3.i).z4 + (4+5.i).z3 + (6+7.i).z2 + (8+9.i).z + (10+11.i) = 0

-Store   2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 , 10 , 11   into registers  R01 to R10

and  XEQ "P4Z"  >>>>   0.268462211      RDN      -1.342085496     whence  z1 = 0.268462211 -1.342085496.i  =  R01+i.R02    ( in 72 seconds )
RDN  -1.226930040      RDN      -0.762095736     whence  z2 = -1.226930040 -0.762095736.i =  R03+i.R04
[ RCL 05 ]   0.361168050   [ RCL 06 ]   1.374351077     whence  z3 =  0.361168050 + 1.374351077.i = R05+i.R06
[ RCL 07 ]  -1.171930991  [ RCL 08 ]   0.883676309     whence  z4 =  -1.171930991 + 0.883676309.i = R07+i.R08

Note:   Low accuracy is to be expected for multiple roots ( use MSRZ first )

d) Evaluating a  complex Polynomial

-"PLZ" computes   p(z) = (a0+i.b0).zn + (a1+i.b1).zn-1 + ......... + (an-1+i.bn-1).z + (an+i.bn)       with  z = x + i.y

Data Registers:                                                                  ( Registers Rbb thru Ree are to be initialized before executing "PLZ" )

•  Rbb = a0 •  Rbb+1 = b0 , •  Rbb+2 = a1 •  Rbb+3 = b1 , ....... , •  Ree-1 = an •  Ree = bn
Flags: /
Subroutines:  /

 01  LBL "PLZ" 02  R-P 03  STO M 04  RDN 05  STO N 06  CLX 07  ENTER^ 08  LBL 01 09  R-P 10  RCL M 11  * 12  X<>Y 13  RCL N 14  + 15  X<>Y 16  P-R 17  RCL IND Z 18  + 19  X<>Y 20  ISG Z 21  RCL IND Z 22  + 23  X<>Y 24  ISG Z 25  GTO 01 26  CLA 27  END

( 44 bytes )

 STACK INPUTS OUTPUTS Z bbb.eee / Y y B X x A

where  p(z) = p(x+i.y) = A+i.B

Example:       p(z) = (1+2i).z3+(4-7i).z2+(2-3i).z+(1-4i)   Evaluate  p(-3+5i)

-If we choose bbb = 001

1 STO 01       4 STO 03      2 STO 05      1 STO 07
2 STO 02     -7 STO 04    -3 STO 06     -4 STO 08    ( control number = 1.008 )

1.008  ENTER^
5      ENTER^
-3     XEQ "PLZ"  >>>>   -86  X<>Y  413   whence   p(-3+5i) =  -86 + 413.i     ( in 6.5 seconds )

Note:

-If you have a "Z*Z" M-code routine in your HP-41 ( cf for instance "A few M-Code Routines for the HP-41" ), "PLZ" may be improved as follows:

 01  LBL "PLZ" 02  STO M 03  RDN 04  STO N 05  X<>Y 06  STO O 07  CLST 08  LBL 01 09  RCL N 10  RCL M 11  Z*Z 12  RCL IND O 13  + 14  ISG O 15  X<>Y 16  RCL IND O 17  + 18  X<>Y 19  ISG O 20  GTO 01 21  CLA 22  END

( 41 bytes )

-With the same example, it yields:

1.008  ENTER^
5      ENTER^
-3     XEQ "PLZ"  >>>>   -86  X<>Y  413   whence   p(-3+5i) =  -86 + 413.i     ( in 2.2 seconds instead of 6.5 seconds )

-Moreover, there is no roundoff error.
-If the polynomial itself has only real coefficients, the routine may be simplified further:

 01  LBL "PLXZ" 02  STO M 03  RDN 04  STO N 05  X<>Y 06  STO O 07  CLST 08  LBL 01 09  RCL N 10  RCL M 11  Z*Z 12  RCL IND O 13  + 14  ISG O 15  GTO 01 16  CLA 17  END

( 35 bytes )

-For instance,   p(z) = -6 z3 + 5 z2 + 2 z - 4   Evaluate  p(-3+5i)

-6  STO 01   5  STO 02   2  STO 03   -4  STO 04     ( control number = 1.004 )

1.004  ENTER^
5      ENTER^
-3     XEQ "PLXZ"  >>>>   -1278  X<>Y  -200   whence   p(-3+5i) =  -1278 - 200 i

e) First derivative of a complex Polynomial

-"dPLZ" calculates  p'(z)  where   p(z) = (a0+i.b0).zn + (a1+i.b1).zn-1 + ......... + (an-1+i.bn-1).z + (an+i.bn)        ( z = x + i.y )

Data Registers:                           ( Registers Rbb thru Ree are to be initialized before executing "dPLZ" )

•  Rbb = a0 •  Rbb+1 = b0 , •  Rbb+2 = a1 •  Rbb+3 = b1 , ....... , •  Ree-1 = an •  Ree = bn
Flags: /
Subroutines:  /

 01  LBL "dPLZ" 02  R-P 03  STO M 04  RDN 05  STO N 06  X<>Y 07  STO O 08  FRC 09   E3 10  * 11  RCL O        12  INT 13  - 14  1 15  - 16  2 17  / 18  STO P 19  CLST 20  LBL 01 21  R-P 22  RCL M        23  * 24  X<>Y 25  RCL N 26  + 27  X<>Y 28  P-R 29  RCL IND O 30  RCL P 31  * 32  ST+ Y 33  X<> L 34  ISG O 35  RCL IND O 36  * 37  ST+ Z 38  ISG O 39  RDN 40  DSE P         41  GTO 01 42  CLA 43  END

( 70 bytes )

 STACK INPUTS OUTPUTS Z bbb.eee / Y y B X x A

where  p'(z) = p'(x+i.y) = A+i.B

Example:   With the same polynomial,  evaluate  p'(-3+5i)

1.008  ENTER^
5     ENTER^
-3    XEQ "dPLZ"  >>>>  180  X<>Y  -107   whence  p'(-3+5i) = 180-107.i

Note:

-If you have "Z*Z" & "DEG F" M-code routines in your HP-41 ( cf for instance "A few M-Code Routines for the HP-41" ), "dPLZ" may be improved as follows:
-Line 07 may be replaced by    FRC   E3   *   RCL O   INT   -

 01  LBL "dPLZ"  02  STO M 03  RDN 04  STO N 05  X<>Y 06  STO O 07  DEG F 08  1 09  - 10  2 11  / 12  STO P         13  CLST 14  LBL 01 15  RCL N 16  RCL M 17  Z*Z 18  RCL IND O 19  RCL P 20  * 21  + 22  ISG O 23  X<>Y 24  RCL IND O 25  RCL P 26  * 27  + 28  ISG O 29  X<>Y 30  DSE P 31  GTO 01       32  CLA 33  END

( 58 bytes )

f) Second derivative of a complex Polynomial

-"d2PLZ" calculates  p''(z)  where   p(z) = (a0+i.b0).zn + (a1+i.b1).zn-1 + ......... + (an-1+i.bn-1).z + (an+i.bn)        ( z = x + i.y )

Data Registers:                                                      ( Registers Rbb thru Ree are to be initialized before executing "d2PLZ" )

•  Rbb = a0 •  Rbb+1 = b0 , •  Rbb+2 = a1 •  Rbb+3 = b1 , ....... , •  Ree-1 = an •  Ree = bn
Flags: /
Subroutines:  /

Note:  Synthetic register P is used, therefore don't interrupt d2PLZ

 01  LBL "d2PLZ" 02  R-P 03  STO M 04  RDN 05  STO N 06  X<>Y 07  STO O 08  FRC 09   E3 10  * 11  RCL O         12  INT 13  - 14  1 15  - 16  2 17  ST- Y 18  / 19  STO P  20  CLST 21  LBL 01 22  R-P 23  RCL M         24  * 25  X<>Y 26  RCL N 27  + 28  X<>Y 29  P-R 30  RCL P 31  ENTER^ 32  ISG X 33  CLX 34  * 35  RCL IND O 36  X<>Y 37  * 38  ST+ Y 39  X<> L 40  ISG O 41  RCL IND O 42  * 43  ST+ Z 44  ISG O 45  RDN 46  DSE P 47  GTO 01 48  CLA 49  END

( 79 bytes )

 STACK INPUTS OUTPUTS Z bbb.eee / Y y B X x A

where  p''(z) = p''(x+i.y) = A+i.B

Example:   With the same polynomial,  evaluate  p''(-3+5i)

1.008  ENTER^
5     ENTER^
-3    XEQ "d2PLZ"  >>>>  -70  X<>Y  -20   whence  p''(-3+5i) = -70-20.i

g) Complex Polynomial Roots

-"PLRZ" solves the equation  p(z) = 0  where   p(z) = (a0+i.b0).zn + (a1+i.b1).zn-1 + ......... + (an-1+i.bn-1).z + (an+i.bn)        ( z = x + i.y )
-The Newton's formula is applied , which is probably not quite rigorous for complex equations.

-A better formula is used in the Laguerre's method.

Data Registers:   R00 thru R06: temp                           ( Registers Rbb thru Ree are to be initialized before executing "PLRZ" )

•  Rbb = a0 •  Rbb+1 = b0 , •  Rbb+2 = a1 •  Rbb+3 = b1 , ....... , •  Ree-1 = an •  Ree = bn     bb > 06
Flags: /
Subroutines:   "PLZ" "dPLZ"

 01  LBL "PLRZ" 02  STO 01 03  RDN 04  STO 02 05  X<>Y 06  STO 00 07  STO 05 08  2 09  + 10  STO 06 11  LBL 01 12  VIEW 01 13  RCL 05 14  RCL 02 15  RCL 01 16  XEQ "dPLZ" 17  R-P 18  STO 03 19  X<>Y 20  STO 04 21  RCL 05 22  RCL 02 23  RCL 01 24  XEQ "PLZ"  25  R-P 26  RCL 03 27  / 28  X<>Y 29  RCL 04 30  - 31  X<>Y 32  P-R 33  ST- 01 34  X<>Y 35  ST- 02 36  LASTX 37  3 E-9     38  X T 52  + 53  X<>Y 54  P-R 55  RCL IND 03 56  + 57  STO IND 03 58  X<>Y 59  ISG 03 60  RCL IND 03 61  + 62  STO IND 03 63  X<>Y 64  ISG 03 65  GTO 02 66  RCL 01 67  STO IND 03 68  ISG 03 69  CLX 70  RCL 02 71  STO IND 03 72  ISG 06 73  ISG 06 74  GTO 01 75  RCL 00 76  2 77  + 78  CLD 79  END

( 126 bytes )

 STACK INPUTS OUTPUTS Z bbb.eee / Y y0 / X x0 2+bbb.eee

where  bbb.eee  is the control number of the polynomial ( bbb > 006 )  and   z0 = x0+i.y0 is an initial guess.

Example:   Find the 3 roots of   p(z) = (1+2i).z3+(4-7i).z2+(2-3i).z+(1-4i)

-If we choose  bbb = 007    1 STO 07       4 STO 09      2 STO 11      1 STO 13
2 STO 08      -7 STO 10    -3 STO 12     -4 STO 14    ( control number = 7.014 )   and  with   z0 = 1 + i

7.014  ENTER^
1      ENTER^     XEQ "PLRZ"     the successive x-approximations are displayed and 3mn  later,
we get  9.014 = the control number of the solutions ( in R09 R10 R11 R12 R13 R14 )

R09 =   2.473329599    R10 =  2.966260806     whence   z1 =  2.473329599 + 2.966260806.i
R11 = -0.291909945    R12 = -0.629741372     whence   z2 =  -0.291909945 - 0.629741372.i
R13 = -0.181419655    R14 =  0.663480568     whence    z3 = -0.181419655 + 0.663480568.i

N.B:

-If all the coefficients are real, choose  y0 # 0   ( otherwise, no complex root could be found )
-If you get "DATA ERROR" at line 30 ( meaning  p'(z) = 0 )  change registers R01 & R02 and press  XEQ 01

-If you have "Z*Z" & "Z/Z" M-code functions, "PLRZ" may be improved:
-Line 31 RND   the accuracy is controlled by the display settings

 01  LBL "PLRZ"  02  STO 01 03  RDN 04  STO 02 05  X<>Y 06  STO 00 07  STO 05 08  2 09  + 10  STO 06 11  LBL 01 12  VIEW 01 13  RCL 05 14  RCL 02 15  RCL 01 16  XEQ "dPLZ" 17  STO 03 18  X<>Y 19  STO 04 20  RCL 05 21  RCL 02 22  RCL 01 23  XEQ "PLZ" 24  RCL 04 25  RCL 03 26  Z/Z 27  ST- 01 28  X<>Y 29  ST- 02 30  R-P 31  RND   32  X#0? 33  GTO 01        34  2 E-3 35  ST- 05 36  RCL 05 37  STO 03 38  CLST 39  LBL 02 40  RCL 02 41  RCL 01 42  Z*Z 43  RCL IND 03 44  + 45  STO IND 03 46  ISG 03 47  X<>Y 48  RCL IND 03 49  + 50  STO IND 03 51  X<>Y 52  ISG 03 53  GTO 02 54  RCL 01 55  STO IND 03 56  ISG 03 57  CLX 58  RCL 02 59  STO IND 03 60  ISG 06 61  ISG 06 62  GTO 01 63  2 64  RCL 00 65  + 66  CLD 67  END

( 108 bytes )

h) Laguerre's method

-We solve again   p(z) = (a0+i.b0).zn + (a1+i.b1).zn-1 + ......... + (an-1+i.bn-1).z + (an+i.bn)  = 0      ( z = x + i.y )
-In the Laguerre's method, the successive approximations are given by the formula  zk+1 = -n.p(zk)/[p'(zk)+E.((n-1)2(p'(zk))2-n(n-1).p(zk).p''(zk))1/2]
where  E = +1 or -1  is choosen to give the larger denominator.

-Convergence order = 3  for single roots.

Data Registers:   R00 thru R10: temp                           ( Registers Rbb thru Ree are to be initialized before executing "LAGZ" )

•  Rbb = a0 •  Rbb+1 = b0 , •  Rbb+2 = a1 •  Rbb+3 = b1 , ....... , •  Ree-1 = an •  Ree = bn        bb > 10

Flags: /
Subroutines:   "PLZ" "dPLZ" "d2PLZ"

-Lines 104-138 are three-byte  GTO 01

 01  LBL "LAGZ"   02  STO 01   03  RDN   04  STO 02   05  X<>Y   06  STO 00   07  STO 09   08  FRC   09   E3   10  *   11  RCL 00   12  INT   13  -   14  DSE X   15  2   16  /   17  STO 10   18  LBL 01   19  VIEW 01   20  RCL 09   21  RCL 02   22  RCL 01   23  XEQ "PLZ"   24  R-P   25  STO 03   26  X<>Y   27  STO 04   28  RCL 09   29  RCL 02 30  RCL 01   31  XEQ "dPLZ"   32  STO 05   33  X<>Y   34  STO 06   35  X<>Y   36  R-P   37  STO 07   38  X<>Y   39  STO 08   40  RCL 09   41  RCL 02   42  RCL 01   43  XEQ "d2PLZ"   44  R-P   45  1   46  RCL 10   47  ST* 03   48  -   49  *   50  RCL 03   51  *   52  X<>Y   53  RCL 04   54  +   55  X<>Y   56  P-R   57  RCL 10   58  1 59  -   60  RCL 07   61  *   62  X^2   63  RCL 08          64  ST+ X   65  X<>Y   66  P-R   67  ST+ Z   68  X<> T   69  +   70  X<>Y   71  R-P   72  .5   73  ST* Z   74  Y^X   75  P-R   76  RCL 06   77  RCL Z   78  ST- 06   79  +   80  RCL 05   81  RCL Z   82  ST- 05   83  +   84  R-P   85  X<>Y   86  RCL 06   87  RCL 05 88  R-P   89  R^   90  XY 100  ST- 02 101  LASTX 102  3 E-9 103  X T 117  + 118  X<>Y 119  P-R 120  RCL IND 03 121  + 122  STO IND 03 123  X<>Y 124  ISG 03 125  RCL IND 03 126  + 127  STO IND 03 128  X<>Y 129  ISG 03 130  GTO 02 131  RCL 01 132  STO IND 03 133  ISG 03 134  CLX 135  RCL 02 136  STO IND 03 137  DSE 10 138  GTO 01  139  RCL 00 140  2 141  + 142  CLD 143  END

( 209 bytes )

 STACK INPUTS OUTPUTS Z bbb.eee / Y y0 / X x0 2+bbb.eee

where  bbb.eee  is the control number of the polynomial ( bbb > 010 )  and   z0 = x0+i.y0 is an initial guess.

Example:   Find the 5 roots of   p(z) = (1+2i).z5+(4-7i).z4+(2-3i).z3+(1-4i).z2+(3+i).z+(7+2i)

-If we choose  bbb = 011

1 STO 11       4 STO 13      2 STO 15      1 STO 17      3 STO 19      7 STO 21
2 STO 12      -7 STO 14    -3 STO 16     -4 STO 18      1 STO 20      2 STO 22           ( control number = 11.022 )

and  with  z0 = 1 + i

11.022  ENTER^
1      ENTER^     XEQ "LAGZ"     the successive x-approximations are displayed and 9mn25s  later,
we get  13.022 = the control number of the solutions:

z1 =  2.502865969 + 2.951734364.i        ( R13 & R14 )
z2 =  0.698862128 - 0.516058660.i         ( R15 & R16 )
z3 = -0.549300937 - 0.866560849.i        ( R17 & R18 )
z4 = -0.778640429 + 0.313576664.i       ( R19 & R20 )
z5 =  0.126213268 + 1.117308483.i        ( R21 & R22 )

Notes:

-If you "DATA ERROR"  line 92 , change registers R01 & R02 and press  XEQ 01
-"PLRZ" solves the same equation in 12mn46s

i) Multiple Roots

-The same method used in "MSR" is now applied to complex polynomials:

Data Registers:  R00 thru R10: temp - when the program stops, R02 = R09 = the control number of GCD(p;p')

( Registers Rbb thru Ree are to be initialized before executing "MSRZ" )

•  Rbb = a0 •  Rbb+1 = b0 , •  Rbb+2 = a1 •  Rbb+3 = b1 , ....... , •  Ree-1 = an •  Ree = bn     bb > 10

Flags: /
Subroutines:  "DIVZ"

 01  LBL "MSRZ"   02  ENTER^   03  STO 08   04  FRC   05   E3   06  *   07  RCL 08   08  INT   09  -   10  STO 01   11  1   12  ST+ 01   13  -   14  2   15  /   16  X<> 01   17  .1   18  %   19  +   20  +   21  STO 09   22  LASTX   23  + 24  2 E-3   25  -   26  STO 10   27  RCL 08   28  RCL 09   29  LBL 00   30  RCL IND Y   31  STO IND Y   32  ISG Y   33  RDN   34  ISG Y   35  GTO 00   36  RCL 08   37  RCL 10   38  LBL 01   39  RCL IND Y   40  RCL 01   41  *   42  STO IND Y   43  ISG Z   44  ISG Y   45  CLX   46  RCL IND Z 47  LASTX   48  *   49  STO IND Y   50  ISG Z   51  ISG Y   52  RDN   53  DSE 01   54  GTO 01   55  LBL 02   56  RCL 09   57  RCL 10   58  XEQ "DIVZ"   59  CLX   60  2   61  -   62  STO 09   63  LBL 03   64  ISG 09   65  ISG 09   66  X<0?   67  GTO 04   68  RCL 09   69  ISG X 70  RCL IND X   71  RCL IND 09   72  R-P   73   E-4   74  X>Y?   75  GTO 03   76  LBL 04   77  RCL 10   78  X<> 09   79  STO 10   80  1   81  -   82  ISG X   83  GTO 02   84  RCL 09   85  STO 01   86  ISG X   87  RCL IND X   88  RCL IND 01   89  R-P   90  STO 02   91  X<>Y   92  STO 03 93  LBL 05   94  RCL 01   95  ISG 01   96  RCL IND 01   97  RCL IND Y   98  R-P   99  RCL 02 100  / 101  X<>Y 102  RCL 03 103  - 104  X<>Y 105  P-R 106  STO IND Z 107  X<>Y 108  STO IND 01 109  ISG 01 110  GTO 05 111  RCL 08 112  RCL 09 113  XEQ "DIVZ" 114  END

( 177 bytes / SIZE > 6n+14 )

 STACK INPUTS OUTPUTS X (bbb.eee)p (bbb.eee)s

( bbb > 010 )

Example:    Find all the roots of the polynomial  p(z) = z5+(1-12.i).z4+(-62-6.i).z3+(-10+170.i).z2+(245-10.i).z+(-31-142.i)

-If we choose  bbb = 011

1  STO 11     1   STO 13   -62  STO 15   -10  STO 17   245  STO 19    -31  STO 21
0  STO 12   -12  STO 14    -6   STO 16   170 STO 18   -10   STO 20  -142  STO 22     ( control number = 11.022 )

11.022  XEQ "MSRZ"  >>>>  11.016  ( in 83 seconds )

R11 = 1       R13 =  1    R15 =  -8
R12 = 0       R14 = -5    R16 =  -1

-Whence    s(z) =  z2+(1-5.i).z+(-8-i)         ( All the coefficients of p(z) are integers, so we can be sure all the coefficients of s(z) are integers too )

1  ENTER^
-5  ENTER^
-8  ENTER^
-1  XEQ "P2Z" >>>  1  RDN  2  RDN  -2  RDN  3    therefore, the 2 solutions are  1+2.i  and -2+3.i      ( Actually,  p(z) = ( z-1-2.i )3.( z + 2 - 3.i )2 )

-We have  R02 = 27.034   and the contents of registers R27 to R34 give  GCD(p;p') = z3-7.i.z2+(-19+2.i).z+( 6+17.i)

j) Euclidean Division

-Let    a(z) = (a0+i.b0).zn + (a1+i.b1).zn-1 + ......... + (an-1+i.bn-1).z + (an+i.bn)
and    b(z) = (a'0+i.b'0).zm + (a'1+i.b'1).zm-1 + ......... + (a'm-1+i.b'm-1).z + (a'm+i.b'm)        ( z = x + i.y )

-"DIVZ"  calculates  q(z) and r(z)  such that   a = b.q + r  with  deg(r) < deg(b)

Data Registers:     R00 thru R07: temp                ( Registers Rbb thru Ree & Rbb' thru Ree' are to be initialized before executing "DIVZ" )

•  Rbb = a0   •  Rbb+1 = b0  , •  Rbb+2 = a •  Rbb+3 = b1  , ....... ,  •  Ree-1 = an  •  Ree = bn        bb > 07
•  Rbb' = a'0 •  Rbb'+1 = b'0 , •  Rbb'+2 = a'1 •  Rbb'+3 = b'1 , ....... , •  Ree'-1 = a'n •  Ree' = b'n      bb' > 07
Flags: /
Subroutines:  /

 01  LBL "DIVZ"  02  STO 02 03  FRC 04  X<>Y 05  STO 01 06  FRC 07  X<>Y 08  - 09   E3 10  * 11  RCL 01 12  INT 13  STO 03 14  - 15  RCL 02 16  INT 17  + 18  2 19  / 20  STO 00 21  ISG 00 22  LBL 01 23  RCL 01 24  STO 04 25  RCL 02 26  STO 05 27  ISG Y 28  RCL IND Y 29  RCL IND 01 30  R-P 31  ISG Z 32  RCL IND Z 33  RCL IND 02 34  R-P 35  ST/ Z 36  X<> T 37  X<>Y 38  - 39  STO 07 40  X<>Y 41  STO 06 42  P-R 43  STO IND 04 44  ISG 04 45  X<>Y 46  STO IND 04 47  ISG 04 48  CLX 49  ISG 05 50  ISG 05 51  FS? 30 52  GTO 03 53  LBL 02 54  RCL IND 05 55  ISG 05 56  RCL IND 05 57  X<>Y 58  R-P 59  RCL 06 60  * 61  X<>Y 62  RCL 07 63  + 64  X<>Y 65  P-R 66  ST- IND 04 67  ISG 04 68  X<>Y 69  ST- IND 04 70  ISG 04 71  CLX 72  ISG 05 73  GTO 02 74  LBL 03 75  2 76  ST+ 01 77  DSE 00 78  GTO 01        79  RCL 01 80  RCL 03 81  RCL 01 82  INT 83  1 84  - 85   E3 86  / 87  + 88  END

( 128 bytes )

 STACK INPUTS OUTPUTS Y (bbb.eee)a (bbb.eee)r X (bbb.eee)b (bbb.eee)q

where         (bbb.eee)a = the control number of the dividend                        (bbb.eee)r  = the control number of the remainder    bbb > 007
(bbb.eee)b = the control number of the divisor                           (bbb.eee)q = the control number of the quotient

Example:     a(z) =  (-4+7.i).z5+(-15+12.i).z4+(-34+33.i).z3+(-48+11.i).z2+(-16+13.i).z+(-12+3.i)
b(z) = (1+2.i).z2+(2+3.i).z+(4+7.i)                           Find   q(z) and r(z)

-For instance:       -4 STO 08    -15 STO 10   -34 STO 12   -48 STO 14   -16 STO 16   -12 STO 18
7  STO 09     12  STO 11    33 STO 13    11 STO 15    13 STO 17      3  STO 19          ( control number = 8.019 )

and            1  STO 21    2  STO 23    4  STO 25
2  STO 22    3  STO 24    7  STO 26       ( control number = 21.026)

8.019  ENTER^
21.026  XEQ "DIVZ"  >>>> ( in 27 seconds )    8.015  X<>Y  16.019        ( q and r have replaced a  ;  b is unchanged )

R08 = 2   R09 = 3   R10 = -2   R11 = 4   R12 = 1   R13 = 3   R14 = -1  R15 = 2   whence  q(z) = (2+3.i)z3+(-2+4.i).z2+(1+3.i).z+(-1+2.i)
R16 = 9   R17 = -7  R18 = 6   R19 = 2              whence  r(z) = (9-7.i).z + (6+2.i)

Notes:

-When b(z) is constant, (bbb.eee)r = 1+eee.eee
-Roundoff errors may produce small leading coefficients instead of zero in the remainder.

"DIVZ" doesn't work if  deg(a) < deg(b)    but in this case,  q = 0 and r = a

-Like many other routines, "DIVZ" may be improved by the M-code routines "Z*Z"  "Z/Z"  "DEG F"
-Moreover, the following variant deletes the zero leading-coefficients.

-Do not key in lines 62 to 82 if  "DIVZ" is called by "MSRZ" as a subroutine ( or use a flag to skip these lines )
-Replace line 67 by    X<>Y   CLX   E-7   X<Y?    if you want to delete the leading coefficients when they are smaller than  E-7  in magnitude

 01  LBL "DIVZ"  02  STO 02 03  DEG F 04  X<>Y 05  STO 01 06  DEG F 07  X<>Y 08  - 09  2 10  / 11  STO 00 12  RCL 01 13  INT 14  STO 03 15  ISG 00 16  LBL 01 17  RCL 01 18  STO 04 19  RCL 02 20  STO 05 21  ISG Y 22  RCL IND Y 23  RCL IND 01 24  ISG Z 25  RCL IND Z 26  RCL IND 02 27  Z/Z 28  STO 06 29  STO IND 04 30  ISG 04 31  X<>Y 32  STO 07 33  STO IND 04 34  ISG 04 35  CLX 36  ISG 05 37  ISG 05 38  FS? 30 39  GTO 00 40  LBL 02 41  RCL IND 05 42  ISG 05 43  RCL IND 05 44  X<>Y 45  RCL 07 46  RCL 06 47  Z*Z 48  ST- IND 04 49  ISG 04 50  X<>Y 51  ST- IND 04 52  ISG 04 53  CLX 54  ISG 05 55  GTO 02 56  LBL 00 57  2 58  ST+ 01 59  DSE 00 60  GTO 01 61  RCL 01 62  LBL 03 63  RCL IND X 64  ISG Y 65  RCL IND Y 66  R-P     67  X#0?       68  GTO 00 69  X<> Z 70  ISG X 71  GTO 03 72  2 73  - 74  STO Y 75  0 76  STO IND Z 77  ISG Z 78  STO IND Z 79  LBL 00 80  X<> Z 81  1 82  - 83  RCL 03 84  RCL 01 85  INT 86  DSE X 87   E3 88  / 89  + 90  END

( 144 bytes )

k) Polynomial Multiplication

-Let    a(z) = (a0+i.b0).zn + (a1+i.b1).zn-1 + ......... + (an-1+i.bn-1).z + (an+i.bn)
and    b(z) = (a'0+i.b'0).zm + (a'1+i.b'1).zm-1 + ......... + (a'm-1+i.b'm-1).z + (a'm+i.b'm)

"PROZ" calculates and stores the coefficients of  p(z) = a(z).b(z)   into contiguous registers. ( you have to specify the first of these registers )

Data Registers:  R00 thru R07: temp                ( Registers Rbb thru Ree & Rbb' thru Ree' are to be initialized before executing "PROZ" )

•  Rbb = a0   •  Rbb+1 = b0  , •  Rbb+2 = a •  Rbb+3 = b1  , ....... ,  •  Ree-1 = an  •  Ree = bn        bb > 07
•  Rbb' = a'0 •  Rbb'+1 = b'0 , •  Rbb'+2 = a'1 •  Rbb'+3 = b'1 , ....... , •  Ree'-1 = a'n •  Ree' = b'n      bb' > 07
Flags: /
Subroutines:  /

-If you don't have an HP-41CX replace line 24  by     0  LBL 00  STO IND Y  ISG Y  GTO 00

 01  LBL "PROZ" 02  ENTER^ 03  R^ 04  STO 01 05  FRC 06  R^ 07  STO 02 08  FRC 09  + 10  X<>Y 11  DSE X 12  RCL 01 13  INT 14  - 15  RCL 02 16  INT 17  - 18   E3 19  / 20  + 21  + 22  STO 00        23  STO 03 24  CLRGX   25  LBL 01 26  RCL 00 27  STO 04 28  RCL 01 29  STO 05 30  RCL 02 31  ISG 02 32  RCL IND 02 33  RCL IND Y 34  R-P 35  STO 06 36  X<>Y 37  STO 07 38  LBL 02 39  RCL IND 05 40  ISG 05 41  RCL IND 05 42  X<>Y 43  R-P 44  RCL 06 45  * 46  X<>Y 47  RCL 07 48  + 49  X<>Y 50  P-R 51  ST+ IND 04 52  ISG 04 53  X<>Y 54  ST+ IND 04 55  ISG 04 56  CLX 57  ISG 05 58  GTO 02 59  2 60  ST+ 00 61  ISG 02 62  GTO 01 63  RCL 03 64  END

( 91 bytes )

 STACK INPUTS OUTPUTS Z (bbb.eee)a Y (bbb.eee)b X bbbp (bbb.eee)p

where         (bbb.eee)a = the control number of  a(z)
(bbb.eee)b = the control number of  b(z)                           (bbb.eee)p = the control number of the product        ( all  bbb > 007 )

Example:     a(z) = (2+3i).z2+(4+7i).z+(1-9i)      b(z) =  (4-6i).z3+(2-3i).z2+(5+7i).z+(1+2i)        Find  p(z) = a(z).b(z)

For instance:

2  STO 08   4  STO 10    1  STO 12                4  STO 14    2  STO 16   5  STO 18   1  STO 20
3  STO 09   7  STO 11   -9  STO 13              -6  STO 15   -3  STO 17   7  STO 19   2  STO 21

( control number = 8.013 )                                          ( control number = 14.021 )

and if we want to have p(z) in registers  R22  R23  ... etc ...

8.013  ENTER^
14.021 ENTER^
22     XEQ "PROZ"  yields  22.033  ( in 27 seconds )

-We find  R22 = 26   R23 = 0   R24 = 71   R25 = 4   R26 = -32   R27 = -11   R28 = -58   R29 = 49   R30 = 58   R31 = -23   R32 = 19   R33 = -7

whence  p(z) = 26.z5+(71+4i).z4+(-32-11i).z3+(-58+49i).z2+(58-23i).z+(19-7i)

Note:

-The following variant uses the M-code routines "Z*Z" and "DEG F"

 01  LBL "PROZ" 02  ENTER^ 03  R^ 04  STO 01 05  DEG F 06  R^ 07  STO 02 08  DEG F 09  + 10  + 11  DSE X 12   E3 13  / 14  + 15  STO 00        16  STO 03 17  CLRGX 18  LBL 01 19  RCL 00 20  STO 04 21  RCL 01 22  STO 05 23  RCL IND 02 24  ISG 02 25  RCL IND 02 26  X<>Y 27  LBL 02 28  RCL IND 05 29  ISG 05 30  RCL IND 05 31  X<>Y 32  Z*Z 33  ST+ IND 04 34  ISG 04 35  RDN 36  ST+ IND 04 37  ISG 04 38  RDN 39  ISG 05 40  GTO 02 41  2 42  ST+ 00 43  ISG 02 44  GTO 01        45  RCL 03 46  END

( 76 bytes )

-Let    a(z) = (a0+i.b0).zn + (a1+i.b1).zn-1 + ......... + (an-1+i.bn-1).z + (an+i.bn)
and    b(z) = (a'0+i.b'0).zm + (a'1+i.b'1).zm-1 + ......... + (a'm-1+i.b'm-1).z + (a'm+i.b'm)

"ADDZ" computes and stores the coefficients of  p(z) = a(z) + b(z)  if F01 is clear
or  p(z) = a(z) - b(z)  if F01 is set          into contiguous registers. ( you have to specify the first of these registers )

Data Registers:  R00 to R03: temp           ( Registers Rbb thru Ree & Rbb' thru Ree' are to be initialized before executing "ADDZ" )

•  Rbb = a0   •  Rbb+1 = b0  , •  Rbb+2 = a •  Rbb+3 = b1  , ....... ,  •  Ree-1 = an  •  Ree = bn        bb > 03
•  Rbb' = a'0 •  Rbb'+1 = b'0 , •  Rbb'+2 = a'1 •  Rbb'+3 = b'1 , ....... , •  Ree'-1 = a'n •  Ree' = b'n      bb' > 03

Flags:  F01  -If flag F01 is clear, "ADDZ" calculates the sum of 2 polynomials
-If flag F01 is set,     ---------------------- difference of 2 polynomials

Subroutines:  "DELZ"

..........................

Simply replace line 66 by XEQ "DELZ"

..........................
66  XEQ "DELZ"
67  END

( 104 bytes )

 STACK INPUTS OUTPUTS Z (bbb.eee)a Y (bbb.eee)b X bbbp (bbb.eee)p

where         (bbb.eee)a = the control number of  a(z)
(bbb.eee)b = the control number of  b(z)                           (bbb.eee)p = the control number of the sum ( if CF 01 ) or the difference ( if SF 01 )

( all  bbb > 003 )

Example:

a(z) = (1+4i).z2+(2-5i)z.+(1+6i)     b(z) = (2-7i).z+(3+4i)    Find  a+b  and  a-b

1 STO 04   4 STO 05   2 STO 06  -5 STO 07   1 STO 08  6 STO 09   ( control number = 4.009 )
2 STO 11  -7 STO 12  3 STO 13   4 STO 14  ( control number = 11.014 )   and if we want to have p(x) in registers R21  R22  ... etc ...

CF 01
4.009  ENTER^
11.014  ENTER^
21     XEQ "ADDZ"  gives  21.026   ( in 7 seconds )   R21 = 1  R22 = 4  R23 = 4  R24 = -12  R25 = 4  R26 = 10

whence  a(z)+b(z) = (1+4i).z2+(4-12i).z+(4+10i)

SF 01
4.009  ENTER^
11.014  ENTER^
21     XEQ "ADDZ"  gives  21.026   ( in 7 seconds )   R21 = 1  R22 = 4  R23 = 0  R24 = 2  R25 = -2  R26 = 2

whence  a(z)-b(z) = (1+4i).z2+2i.z+(-2+2i)

-This short routine deletes tiny dominant coefficients resulting from roundoff errors or occuring in additions or subtractions.

Data Registers:                                                           ( Registers Rbb thru Ree are to be initialized before executing "DELZ" )

•  Rbb = a0   •  Rbb+1 = b0  , •  Rbb+2 = a •  Rbb+3 = b1  , ....... ,  •  Ree-1 = an  •  Ree = bn
Flags: /
Subroutines: /

 01  LBL "DELZ" 02  LBL 01 03  RCL IND X 04  ISG Y 05  RCL IND Y 06  R-P 07   E-7     08  X

( 45 bytes )

 STACK INPUTS OUTPUTS X (bbb.eee)p (bbb'.eee)p

Example:      If   p(z) = (10-9+10-8i).z2+(3+2i).z+(3+7i)    and   R01 = 10-9   R02 = 10-8  R03 = 3  R04 = 2  R05 = 3  R06 =7    ( control number = 1.006 )

1.006   XEQ "DELZ"  yields   3.006   meaning  p(z) = (3+2i).z+(3+7i)

Note:

-If p(x) = (10-9+10-8i)  is stored in R01 & R02    1.002  XEQ "DELZ" gives  1.002  but 0 is stored into registers R01 & R02

3°) Three short routines

a) Extremum of the curve y = a.x2+b.x+c

Data Registers:  /
Flags:  /
Subroutines:  /

 01  LBL "EX2"  02  X<>Y  03  2  04  /  05  CHS  06  ENTER^  07  X^2  08  R^  09  ST/ Z  10  /  11  ST- Z  12  RDN  13  END

( 23 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Z a c Y b y X c x L / a

where A(x;y)  is the extremum

Example:     y = 3.x2-12.x +7

3   ENTER^
-12  ENTER^
7  XEQ "EX2"  yields   2   X<>Y   -5    whence  A(2;-5)

Exercise:   Write a 22-byte variant

b) Extrema of the curve y = a.x3+b.x2+c.x+d

Data Registers:  /
Flags:  F00
Subroutines:  "P2"

Note:  Synthetic register P is used, therefore don't interrupt  EX3

 01  LBL "EX3" 02  STO M 03  CLX 04  3 05  R^ 06  STO N 07  * 08  X<> Z 09  STO O 10  ST+ X 11  X<>Y 12  STO P   13  XEQ "P2" 14  RCL N 15  RCL Z 16  * 17  RCL O 18  + 19  R^ 20  * 21  RCL P 22  + 23  R^ 24  * 25  RCL M     26  + 27  X<> N 28  RCL Y 29  * 30  RCL O 31  + 32  RCL Y 33  * 34  RCL P 35  + 36  RCL Y 37  * 38  RCL M     39  + 40  RCL N 41  R^ 42  R^ 43  CLA 44  FC?C 00   45  X=Y? 46  SF 46 47  X<> Z 48  X<>Y 49  END

( 82 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS T a y2 Z b x2 Y c y1 X d x1

where  A(x1,y1)  and  B(x2,y2)  are the 2 extrema
If you get "NON EXISTENT" there is no extremum

Example1:     y = 2.x3 - 3.x2 - 12.x + 1

2   ENTER^
-3   ENTER^
-12  ENTER^
1   XEQ "EX3"  >>>>   -1  RDN  8   RDN  2   RDN  -19    whence  A(-1;8) &  B(2;-19)

Example2:     y = x3 - 3.x2 + 3.x + 4   gives "NON EXISTENT"  ( y'(x) = 0 has only one solution )

Example3:     y = x3 + x2 + x + 1  also produces "NON EXISTENT"  ( y'(x) = 0 has no real solution )

c) Center of symmetry of the curve:    y = a.x3+b.x2+c.x+d

Data Registers:  /
Flags:  /
Subroutines:  /

 01  LBL "INFL"  02  X<> Z  03  R^  04  /  05  3  06  ST/ Y  07  X<> L  08  ST+ X  09  X<>Y  10  CHS  11  ST* Y  12  ST* Y  13  RDN  14  -  15  R^  16  *  17  +  18  X<>Y  19  END

( 34 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS T a / Z b / Y c y X d x

where  A(x;y)  is the center of symmetry

Example:     y = 2.x3 + 3.x2 + 4.x + 5

2   ENTER^
3   ENTER^
4   ENTER^
5   XEQ "INFL"  yields  -0.5  X<>Y  3.5   whence  A(-1/2,7/2)