hp41programs

Polynomials for the HP-41

Polynomials for the HP-41


Overview
 

1°) Real Polynomials

     a) Quadratic equations
     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
     k) Addition & Subtraction
     l)  Deleting tiny leading coefficients
    m) Composition of 2 Polynomials

2°) Complex Polynomials

     a) Simplified quadratic equations   ( leading coefficient = 1 )
     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
     l)  Addition & Subtraction
     m) Deleting tiny leading coefficients

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
 

        a) Quadratic Equations
 

-"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  X<Y?
  14  X<>Y
  15  RDN
  16  X<Y?
  17  X<>Y
  18  R^
  19  X<Y?
  20  X<>Y
  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<Y?
27  GTO 01       
28   E-3
29  ST- 03
30  RCL 03
31  STO 02
32  CLX
33  LBL 02
34  RCL 01
35  *
36  ST+ IND 02
37  RCL IND 02
38  ISG 02
39  GTO 02
40  RCL 01
41  STO IND 02
42  ISG 04
43  GTO 01
44  RCL 00
45  1
46  +
47  CLD
48  END

 
  ( 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
 

        g) Quadratic Factors:  Bairstow's method
 

-"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  X<Y? 
  69  GTO 01
  70  SIGN
  71  STO 05
  72  RCL 04
  73  5.007
  74  XEQ "DIV"
  75  STO 04
  76  RCL 06
  77  STO IND Z 
  78  ISG Z
  79  RCL 07
  80  STO IND T
  81  RCL 04
  82  2
  83  +
  84  ISG X
  85  GTO 01
  86  CLD
  87  RCL 08
  88  BEEP
  89  RTN
  90  LBL 10
  91  RCL 08
  92  STO 04
  93  FRC
  94   E3
  95  *
  96  RCL 04
  97  INT
  98  -
  99  2
100  MOD
101  X#0?
102  GTO 03
103  RCL IND 04
104  ISG 04
105  GTO 05
106  LBL 03
107  0
108  RCL IND 04
109  ISG 04
110  RCL IND 04
111  X<>Y
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.
 

        k) Addition/Subtraction of 2 Polynomials
 

-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  X<Y?
20  X<>Y
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 )
 

        l) Deleting tiny leading coefficients
 

-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<Y?
07  GTO 02
08  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:
 

        a) Complex quadratic equation
 

-"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<Y?
39  GTO 01       
40  2 E-3
41  ST- 05
42  RCL 05
43  STO 03
44  CLST
45  LBL 02
46  R-P
47  RCL 02
48  RCL 01
49  R-P
50  ST* Z
51  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  X<Y?
  91  RCL Y
  92  ST/ 03
  93  R^
  94  ST- 04
  95  RCL 04       
  96  RCL 03
  97  P-R
  98  ST- 01
  99  X<>Y
100  ST- 02
101  LASTX
102  3 E-9
103  X<Y?
104  GTO 01  
105  2 E-3
106  ST- 09
107  RCL 09
108  STO 03
109  CLST
110  LBL 02
111  R-P
112  RCL 02
113  RCL 01
114  R-P
115  ST* Z
116  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 )
 

        l) Addition/Subtraction of 2 Polynomials
 

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

01  LBL "ADDZ"
..........................

the same listing as "ADD"
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)
 

        m) Deleting tiny leading coefficients
 

-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<Y?
09  GTO 02      
10  R^
11  ISG X
12  GTO 01      
13  2
14  -
15  STO Z
16  0
17  STO IND T
18  ISG T
19  STO IND T
20  LBL 02
21  R^
22  1
23  -
24  END         

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