The Museum of HP Calculators

Polynomials for the HP-41

This program is Copyright © 2004, 2005 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Overview
 

1°) Real Polynomials

     a) Quadratic equations  ( improved )
     b) Cubic equations  ( improved )
     c) Quartic equations  ( improved )
     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

2°) Complex Polynomials

     a) Simplified quadratic equations   ( leading coefficient = 1 )
     b) Cubic equations  ( new )
     c) Quartic equations  ( new )
     d) Evaluating a polynomial
     e) First derivative of a polynomial
     f) Second derivative of a polynomial ( a subroutine used in the Laguerre program )
     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.x²+b.x+c
    b) Extrema of the curve y =  a.x³+b.x²+c.x+d
    c) Center of symmetry of the curve y = a.x³+b.x²+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.x²+b.x+c = 0    with a # 0

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

01  LBL "P2"
02  CF 00
03  X#0?                      Theoretically, lines 03 to 09 are superfluous,
04  GTO 00                  Practically, they are useful if c = 0 to obtain the exact solution  x = 0
05  X<> Z
06  /
07  CHS
08  RTN
09  LBL 00
10  X<> Z
11  ST/ Z
12  ST+ X
13  /
14  CHS
15  ENTER^
16  ENTER^
17  X^2
18  R^
19  -
20  X<0?
21  SF 00
22  ABS
23  SQRT
24  FC? 00
25  ST+ Z
26  FC? 00
27  -
28  X<>Y
29  END

( 46 bytes / SIZE 000 )
 
 
 

STACK	INPUTS	OUTPUTS
Z	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
 

Example:     Solve  2.x² + 3.x - 4 = 0   and   2.x² + 3.x + 4 = 0

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

  b)       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.x³+b.x²+c.x+d  by the Cardan's ( or Tartaglia's ) formulae:  ( with a # 0 )
-First, the term in x² is removed by a change of argument, leading to  x³+p.x+q = 0
-Then, x = u+v  with u.v = -p/3 leads to a quadratic equation in u³

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

  01  LBL "P3"
  02  DEG
  03  CF 01
  04  X#0?                             Lines 04 to 22 are useful to produce exactly  x = 0 if  d = 0
  05  GTO 00                        This is important when "P3" is called as a subroutine by "P4" or some "cosmological" programs
  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	/
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.x³-5.x²-7.x+1 = 0  and   2.x³-5.x²+7.x+1 = 0

  a)   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.

  b)   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   x⁴+a.x³+b.x²+c.x+d = 0   ( if the leading coefficient is not 1 , divide all the equation by this coefficient ).
-First, the term in x³ is removed by a change of argument, leading to  x⁴+p.x²+q.x+r = 0  (E')
-Then, the resolvant  z³+p.z²/2+(p²-4r).z/16-q²/64 = 0 is solved by "P3" and if we call  z₁ , z₂ , z₃  the 3 roots of this equation, the zeros of (E') are

          x = z₁^1/2 sign(-q) +/- ( z₂^1/2+z₃^1/2 )    ;    x = -(z₁^1/2) sign(-q) +/- ( z₂^1/2-z₃^1/2 )

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

  01  LBL "P4"
  02  CF 02
  03  X#0?                                 Lines 03 to 12 are only useful if  d = 0  to compute exactly the solution  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  x⁴-2.x³-35.x²+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  x⁴-5.x³+11.x²-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  x⁴-8.x³+26.x²-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 new version is less simple but much more reliable than the previous program.
 

        d) Evaluating a real Polynomial

-The following program calculates  p(x) = a₀.xⁿ+a₁.x^n-1+ ... + a_n-1.x+a_n  for any given x-value.
-The coefficients   a₀ ; a₁ ;  ...... ; a_n-1 ; a_n  are to be stored into contiguous registers.

Data Registers:   Rbb = a₀ ; Rbb+1 = a₁ ; ....... ; Ree = a_n   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.x³+3.x²-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 LastX
 

        e) First derivative of a Polynomial

-"dPL" calculates  p'(x) = n.a₀.x^n-1+(n-1).a₁.x^n-2+ ... + a_n-1  for any given x-value.
-The coefficients   a₀ ; a₁ ;  ...... ; a_n-1 ; a_n  are to be stored into contiguous registers.
 

Data Registers:   Rbb = a₀ ; Rbb+1 = a₁ ; ....... ; Ree = a_n  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  ""                   (  TEXT 0  or another NOP instruction like  LBL 02 ...  )
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.x³+3.x²-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) = a₀.xⁿ+a₁.x^n-1+ ... + a_n-1.x+a_n = 0  provided all roots are real.
-The coefficients   a₀ ; a₁ ;  ...... ; a_n-1 ; a_n  are to be stored into contiguous registers ( Rbb thru Ree ).

-Starting with an initial guess  x₀ , the Newton's formula   x_k+1 = x_k-p(x_k)/p'(x_k) is applied until  | p(x_k)/p'(x_k) | 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
                      and  Rbb = a₀ ; Rbb+1 = a₁ ; ....... ; Ree = a_n   these n+1 registers are to be initialized before executing "PLR"    ( 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  x₀ is an initial guess.

Example:   Find all the roots of   p(x) = 2.x⁵+3.x⁴-35.x³-10.x²+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  x₀ = 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) = a₀.xⁿ+a₁.x^n-1+ ... + a_n-1.x+a_n  into quadratic factors and solves  p(x) = 0      ( n > 1 )
-If deg(p) is odd, we have   p(x) =    (a₀.x+b).(x²+u₁.x+v₁)...............(x²+u_m.x+v_m)      where m = (n-1)/2
-If deg(p) is even                p(x) = (a₀x²+u₁.x+v₁)(x²+u₂.x+v₂)..........(x²+u_m.x+v_m)      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 (x²+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 = a₀ ; Rbb+1 = b ; Rbb+2 = u₁ ; Rbb+3 = v₁ ; ........ ; Ree-1 = u_m ; Ree = v_m
  If deg(p) is even         Rbb = a₀ ; Rbb+1 = u₁ ; Rbb+2 = v₁ ; ............................ ; Ree-1 = u_m ; Ree = v_m

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

Data Registers:  R00 to R08: temp.   ( R06 = u ; R07 = v ; R08 = bbb.eee )
                      and  Rbb = a₀ ; Rbb+1 = a₁ ; ....... ; Ree = a_n   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 )
 

001  LBL "BRST"
002  STO 06
003  RDN
004  STO 07
005  X<>Y
006  STO 08
007  STO 04
008  LBL 01
009  VIEW 06
010  CLX
011  STO 00
012  STO 01
013  STO 02
014  STO 03
015  RCL 04
016  STO 05
017  LBL 02
018  RCL IND 05
019  RCL 00
020  RCL 07
021  *
022  -
023  RCL 01
024  RCL 06
025  *
026  -
027  X<> 01
028  STO 00
029  RCL 02
030  RCL 07
031  *
032  -
033  RCL 03
034  RCL 06
035  *
036  -
037  X<> 03
038  X<> 02
039  ISG 05
040  GTO 02
041  STO 05
042  RCL 01
043  *
044  RCL 00
045  RCL 02
046  ST* 01
047  *
048  -
049  RCL 03
050  RCL 05
051  *
052  RCL 02
053  X^2
054  -
055  STO 05
056  /
057  ST+ 06
058  RCL 00
059  RCL 03
060  *
061  RCL 01
062  -
063  RCL 05
064  /
065  ST+ 07
066  R-P
067   E-8           ( or another "small" number )
068  X<Y?
069  GTO 01
070  SIGN
071  STO 05
072  RCL 04
073  5.007
074  XEQ "DIV"
075  STO 04
076  RCL 06
077  STO IND Z
078  ISG Z
079  RCL 07
080  STO IND T
081  RCL 04
082  2
083  +
084  ISG X
085  GTO 01
086  CLD
087  RCL 08
088  BEEP
089  STOP
090  LBL 10
091  RCL 08
092  STO 04
093  FRC
094   E3
095  *
096  RCL 04
097  INT
098  -
099  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  u₀ and  v₀  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:     x_k  ;  y_k   if flag F00 is clear
                                                                            or        x_k+i.y_k  ;   x_k-i.y_k  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.x⁵+7.x⁴+20.x³+81.x²+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  u₀ = v₀ = 0

   11.016  ENTER^
     0         ENTER^
                XEQ "BRST"  the successive u_k 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).(x²-2x+10).(x²+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      x⁶-6.x⁵+8.x⁴+64.x³-345.x²+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  u₀ = v₀ = 0

   11.017  ENTER^
     0         ENTER^
                XEQ "BRST"  the successive u_k 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) = (x²+x-12).(x²-4x+13).(x²-3x+2)

-Then     R/S  >>>  ( TONE 9 )   3      X<>Y   -4   with  F00 clear  the first pair of roots are   3 and  -4
              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   3 ; -4 ; 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')
                      and  Rbb = a₀ ; Rbb+1 = a₁ ; ....... ; Ree = a_n   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) = x⁴+2.x³+5.x²+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) = x²+x+2      ( in fact,  p(x) = ( x²+x+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) = x¹⁰-17.x⁹+127.x⁸-549.x⁷+1521.x⁶-2823.x⁵+3557.x⁴-3007.x³+1634.x²-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) = x³-6.x²+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)⁵(x-2)³(x-3)²  )

-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') =  x⁷-11.x⁶+50.x⁵-122.x⁴+173.x³-143.x²+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) = x³+0.0001.x²
   but if you replace line 57 by  E-9 the good result  s(x) = x²+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) = a₀.xⁿ+a₁.x^n-1+ ... + a_n-1.x+a_n  and   b(x) = b₀.x^m+b₁.x^m-1+ ... + b_m-1.x+b_m  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
                             Rbb = a₀  ; Rbb+1 = a₁  ; ....... ; Ree = a_n    these n+1 registers are to be initialized before executing "DIV"     ( bb > 03 )
                      and  Rbb' = b₀ ; Rbb'+1 = b₁ ; ....... ; Ree' = b_m   these m+1 registers are to be initialized before executing "DIV"    ( 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 = ------------------------- divisor                           (bbb.eee)_q = ------------------------- quotient
 
 

Example:     a(x) =  2.x⁵+5.x⁴-21.x³+23.x²+3.x+5      b(x) = 2.x²-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) = x³+4.x²-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) = a₀.xⁿ+a₁.x^n-1+ ... + a_n-1.x+a_n  and   b(x) = b₀.x^m+b₁.x^m-1+ ... + b_m-1.x+b_m   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
                             Rbb = a₀  ; Rbb+1 = a₁  ; ....... ; Ree = a_n    these n+1 registers are to be initialized before executing "PRO"     ( bb > 03 )
                      and  Rbb' = b₀ ; Rbb'+1 = b₁ ; ....... ; Ree' = b_m   these m+1 registers are to be initialized before executing "PRO"    ( bb' > 03 )
Flags: /
Subroutines: /
 

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                If you don't have an HP-41CX replace line 23  by     0  LBL 00  STO IND Y  ISG Y  GTO 00
24  LBL 01
25  RCL 00
26  RCL 01
27  LBL 02
28  RCL IND X
29  RCL IND 02
30  *
31  ST+ IND Z
32  ISG Z
33  RDN
34  ISG X
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 = ---------------------   b(x)                           (bbb.eee)_p = the control number of the product        ( all  bbb > 003 )
 

Example:     a(x) =  2.x⁵+5.x⁴-21.x³+23.x²+3.x+5      b(x) = 2.x²-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.x⁷+4.x⁶-55.x⁵+114.x⁴-84.x³+24.x²-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) = a₀.xⁿ+a₁.x^n-1+ ... + a_n-1.x+a_n  and   b(x) = b₀.x^m+b₁.x^m-1+ ... + b_m-1.x+b_m   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
                             Rbb = a₀  ; Rbb+1 = a₁  ; ....... ; Ree = a_n    these n+1 registers are to be initialized before executing "ADD"     ( bb > 03 )
                      and  Rbb' = b₀ ; Rbb'+1 = b₁ ; ....... ; Ree' = b_m   these m+1 registers are to be initialized before executing "ADD"    ( bb' > 03 )
Flags:  F01  -If flag F01 is clear, "ADD" calculates the sum of 2 polynomials
                     -If flag F01 is set,     -------------------- difference of 2 polynomials
Subroutines:  "DEL"
 

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                If you don't have an HP-41CX replace line 28  by     0  LBL 00  STO IND Y  ISG Y  GTO 00   RCL 00
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 = ---------------------   b(x)                           (bbb.eee)_p = the control number of the sum ( if CF 01 ) or the difference ( if SF 01 )
                                                                                ( all  bbb > 003 )

Example:

    a(x) = 2.x³+4.x²+5.x+6     b(x) = 2.x³-3.x²+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 ...

1°)            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.x³+x²+12.x+7

2°)             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.x²-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:     Rbb = a₀  ; Rbb+1 = a₁  ; ....... ; Ree = a_n  these n+1 registers are to be initialized before executing "DEL"
Flags: /
Subroutines: /
 

01  LBL "DEL"
02  LBL 01
03  RCL IND X
04  ABS
05   E-7                      ( or another "small" number )
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^-9x²+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
 

2°) Complex Polynomials

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

        a) Complex quadratic equation

-"P2Z"  solves   z²+(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  z²-(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   z₁ = 4 + 7.i   and   z₂ = 2 + 3.i      ( in 4 seconds )
 

        b) Complex cubic equation
 

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

Data Registers:            R00 & R09: temp

                                   •  R01 = a    •  R03 = c   •  R05 = e   •  R07 = g
                                   •  R02 = b    •  R04 = d   •  R06 = f   •  R08 = h           ( these 8  registers are to be initialized before executing "P3Z" )

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

Flags: /
Subroutines:  "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/"
  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/"
  19  STO 05
  20  X<>Y
  21  STO 06
  22  RCL 08
  23  RCL 07
  24  RCL 02
  25  RCL 01
  26  XEQ "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-"
  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*"
  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*"
  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*"
  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*"
113  STO 03
114  X<>Y
115  X<> 04
116  RCL 07
117  XEQ "SHZ"
118  RCL 02
119  RCL 01
120  XEQ "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

( 273 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).z³ + (4+5.i).z² + (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   z₁ = -0.179640712 -1.436921967.i  =  R01+i.R02    ( in 35 seconds )
                             RDN   -0.061935981    RDN      1.506099080   -------    z₂ = -0.061935981 + 1.506099080.i = R03+i.R04
                       [ RCL 05 ] -1.527654077 [ RCL 06 ] 0.084669040  -------    z₃ =  -1.527654077 + 0.084669040.i = R05+i.R06
 

        c) Complex quartic equation
 

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

Data Registers:            R00 & R11 to R15: temp

                                   •  R01 = a    •  R03 = c   •  R05 = e   •  R07 = g   •  R09 = k
                                   •  R02 = b    •  R04 = d   •  R06 = f   •  R08 = h    •  R10 = l       ( these 10  registers are to be initialized before executing "P4Z" )

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

Flags: /
Subroutines:  "P2Z"  "P3Z"
                        "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/"
  07  STO 00
  08  X<>Y
  09  STO 09
  10  RCL 04
  11  RCL 03
  12  RCL 02
  13  RCL 01
  14  XEQ "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/"
  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/"
  37  STO 14
  38  X<>Y
  39  STO 15
  40  X<>Y
  41  RCL 11
  42  RCL 10
  43  XEQ "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-"
  64  RCL 09
  65  RCL 00
  66  XEQ "Z*"
  67  RCL 15
  68  RCL 14
  69  XEQ "Z^2"
  70  XEQ "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*"
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/"
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

( 314 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).z⁴ + (4+5.i).z³ + (6+7.i).z² + (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   z₁ = 0.268462211 -1.342085496.i  =  R01+i.R02    ( in 72 seconds )
                             RDN  -1.226930040      RDN      -0.762095736        -------   z₂ = -1.226930040 -0.762095736.i =  R03+i.R04
                    [ RCL 05 ]   0.361168050   [ RCL 06 ]   1.374351077       -------     z₃ =  0.361168050 + 1.374351077.i = R05+i.R06
                    [ RCL 07 ]  -1.171930991  [ RCL 08 ]   0.883676309        -------    z₄ =  -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) = (a₀+i.b₀).zⁿ + (a₁+i.b₁).z^n-1 + ......... + (a_n-1+i.b_n-1).z + (a_n+i.b_n)       with  z = x + i.y
 

Data Registers:   Rbb = a₀ ; Rbb+1 = b₀ ; Rbb+2 = a₁ ; Rbb+3 = b₁ ; ....... ; Ree-1 = a_n ; Ree = b_n
                          these 2n+2 registers are to be initialized before executing "PLZ"
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).z³+(4-7i).z²+(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 )
 

        e) First derivative of a complex Polynomial
 

-"dPLZ" calculates  p'(z)  where   p(z) = (a₀+i.b₀).zⁿ + (a₁+i.b₁).z^n-1 + ......... + (a_n-1+i.b_n-1).z + (a_n+i.b_n)        ( z = x + i.y )

Data Registers:   Rbb = a₀ ; Rbb+1 = b₀ ; Rbb+2 = a₁ ; Rbb+3 = b₁ ; ....... ; Ree-1 = a_n ; Ree = b_n
                          these 2n+2 registers are to be initialized before executing "dPLZ"
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
 

        f) Second derivative of a complex Polynomial

-"d2PLZ" calculates  p''(z)  where   p(z) = (a₀+i.b₀).zⁿ + (a₁+i.b₁).z^n-1 + ......... + (a_n-1+i.b_n-1).z + (a_n+i.b_n)        ( z = x + i.y )

Data Registers:   Rbb = a₀ ; Rbb+1 = b₀ ; Rbb+2 = a₁ ; Rbb+3 = b₁ ; ....... ; Ree-1 = a_n ; Ree = b_n
                          these 2n+2 registers are to be initialized before executing "dPLZ"
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        ( not  STOP  but synthetic  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) = (a₀+i.b₀).zⁿ + (a₁+i.b₁).z^n-1 + ......... + (a_n-1+i.b_n-1).z + (a_n+i.b_n)        ( 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
                              Rbb = a₀ ; Rbb+1 = b₀ ; Rbb+2 = a₁ ; Rbb+3 = b₁ ; ....... ; Ree-1 = a_n ; Ree = b_n     ( bb > 06 )
                          these 2n+2 registers are to be initialized before executing "PLRZ"
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           ( or another "small" number )
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   z₀ = x₀+i.y₀ is an initial guess.

Example:   Find the 3 roots of   p(z) = (1+2i).z³+(4-7i).z²+(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   z₀ = 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   z₁ =  2.473329599 + 2.966260806.i
  R11 = -0.291909945  R12 = -0.629741372     -------    z₂ =  -0.291909945 - 0.629741372.i
  R13 = -0.181419655  R14 =  0.663480568      -------    z₃ = -0.181419655 + 0.663480568.i

N.B:   -If all the coefficients are real, choose  y₀ # 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
 

        h) Laguerre's method

-We solve again   p(z) = (a₀+i.b₀).zⁿ + (a₁+i.b₁).z^n-1 + ......... + (a_n-1+i.b_n-1).z + (a_n+i.b_n)  = 0      ( z = x + i.y )
-In the Laguerre's method, the successive approximations are given by the formula  z_k+1 = -n.p(z_k)/[p'(z_k)+E.((n-1)²(p'(z_k))²-n(n-1).p(z_k).p''(z_k))^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
                              Rbb = a₀ ; Rbb+1 = b₀ ; Rbb+2 = a₁ ; Rbb+3 = b₁ ; ....... ; Ree-1 = a_n ; Ree = b_n     ( bb > 10 )
                          these 2n+2 registers are to be initialized before executing "LAGZ"
Flags: /
Subroutines:   "PLZ" "dPLZ" "d2PLZ"
 

001  LBL "LAGZ"
002  STO 01
003  RDN
004  STO 02
005  X<>Y
006  STO 00
007  STO 09
008  FRC
009   E3
010  *
011  RCL 00
012  INT
013  -
014  DSE X
015  2
016  /
017  STO 10
018  LBL 01
019  VIEW 01
020  RCL 09
021  RCL 02
022  RCL 01
023  XEQ "PLZ"
024  R-P
025  STO 03
026  X<>Y
027  STO 04
028  RCL 09
029  RCL 02
030  RCL 01
031  XEQ "dPLZ"
032  STO 05
033  X<>Y
034  STO 06
035  X<>Y
036  R-P
037  STO 07
038  X<>Y
039  STO 08
040  RCL 09
041  RCL 02
042  RCL 01
043  XEQ "d2PLZ"
044  R-P
045  1
046  RCL 10
047  ST* 03
048  -
049  *
050  RCL 03
051  *
052  X<>Y
053  RCL 04
054  +
055  X<>Y
056  P-R
057  RCL 10
058  1
059  -
060  RCL 07
061  *
062  X^2
063  RCL 08
064  ST+ X
065  X<>Y
066  P-R
067  ST+ Z
068  X<> T
069  +
070  X<>Y
071  R-P
072  .5
073  ST* Z
074  Y^X
075  P-R
076  RCL 06
077  RCL Z
078  ST- 06
079  +
080  RCL 05
081  RCL Z
082  ST- 05
083  +
084  R-P
085  X<>Y
086  RCL 06
087  RCL 05
088  R-P
089  R^
090  X<Y?
091  RCL Y
092  ST/ 03
093  R^
094  ST- 04
095  RCL 04
096  RCL 03
097  P-R
098  ST- 01
099  X<>Y
100  ST- 02
101  LASTX
102  3 E-9                ( or another "small" number )
103  X<Y?
104  GTO 01            ( a three-byte GTO )
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              ( a three-byte GTO )
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   z₀ = x₀+i.y₀ is an initial guess.

Example:   Find the 5 roots of   p(z) = (1+2i).z⁵+(4-7i).z⁴+(2-3i).z³+(1-4i).z²+(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   z₀ = 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:

   z₁ =  2.502865969 + 2.951734364.i        ( R13 & R14 )
   z₂ =  0.698862128 - 0.516058660.i         ( R15 & R16 )
   z₃ = -0.549300937 - 0.866560849.i        ( R17 & R18 )
   z₄ = -0.778640429 + 0.313576664.i       ( R19 & R20 )
   z₅ =  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') )
                             Rbb = a₀ ; Rbb+1 = b₀ ; Rbb+2 = a₁ ; Rbb+3 = b₁ ; ....... ; Ree-1 = a_n ; Ree = b_n     ( bb > 10 )
                         these 2n+2 registers are to be initialized before executing "MSRZ"
Flags: /
Subroutines:  "DIVZ"
 

001  LBL "MSRZ"
002  ENTER^
003  STO 08
004  FRC
005   E3
006  *
007  RCL 08
008  INT
009  -
010  STO 01
011  1
012  ST+ 01
013  -
014  2
015  /
016  X<> 01
017  .1
018  %
019  +
020  +
021  STO 09
022  LASTX
023  +
024  2 E-3
025  -
026  STO 10
027  RCL 08
028  RCL 09
029  LBL 00
030  RCL IND Y
031  STO IND Y
032  ISG Y
033  RDN
034  ISG Y
035  GTO 00
036  RCL 08
037  RCL 10
038  LBL 01
039  RCL IND Y
040  RCL 01
041  *
042  STO IND Y
043  ISG Z
044  ISG Y
045  CLX
046  RCL IND Z
047  LASTX
048  *
049  STO IND Y
050  ISG Z
051  ISG Y
052  RDN
053  DSE 01
054  GTO 01
055  LBL 02
056  RCL 09
057  RCL 10
058  XEQ "DIVZ"
059  CLX
060  2
061  -
062  STO 09
063  LBL 03
064  ISG 09
065  ISG 09
066  X<0?
067  GTO 04
068  RCL 09
069  ISG X
070  RCL IND X
071  RCL IND 09
072  R-P
073   E-4
074  X>Y?
075  GTO 03
076  LBL 04
077  RCL 10
078  X<> 09
079  STO 10
080  1
081  -
082  ISG X
083  GTO 02
084  RCL 09
085  STO 01
086  ISG X
087  RCL IND X
088  RCL IND 01
089  R-P
090  STO 02
091  X<>Y
092  STO 03
093  LBL 05
094  RCL 01
095  ISG 01
096  RCL IND 01
097  RCL IND Y
098  R-P
099  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) = z⁵+(1-12.i).z⁴+(-62-6.i).z³+(-10+170.i).z²+(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) =  z²+(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 )³.( z + 2 - 3.i )² )

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

        j) Euclidean Division

-Let    a(z) = (a₀+i.b₀).zⁿ + (a₁+i.b₁).z^n-1 + ......... + (a_n-1+i.b_n-1).z + (a_n+i.b_n)
 and    b(z) = (a'₀+i.b'₀).z^m + (a'₁+i.b'₁).z^m-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
       Rbb = a₀ ; Rbb+1 = b₀ ; Rbb+2 = a₁ ; Rbb+3 = b₁ ; ....... ; Ree-1 = a_n ; Ree = b_n     these 2n+2 registers are to be initialized before executing "DIVZ"
     Rbb' = a'₀ ; Rbb'+1 = b'₀ ; Rbb'+2 = a'₁ ; Rbb'+3 = b'₁ ; ....... ; Ree'-1 = a'_m ; Ree' = b'_m these 2m+2 registers are to be initialized before executing "DIVZ"
Flags: /                                             ( bb > 07 & bb' > 07 )
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        ( all  bbb > 007 )
                  (bbb.eee)_b = ------------------------- divisor                           (bbb.eee)_q = ------------------------- quotient
 
 

Example:     a(z) =  (-4+7.i).z⁵+(-15+12.i).z⁴+(-34+33.i).z³+(-48+11.i).z²+(-16+13.i).z+(-12+3.i)
                    b(z) = (1+2.i).z²+(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)z³+(-2+4.i).z²+(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
 

        k) Polynomial Multiplication
 

-Let    a(z) = (a₀+i.b₀).zⁿ + (a₁+i.b₁).z^n-1 + ......... + (a_n-1+i.b_n-1).z + (a_n+i.b_n)
 and    b(z) = (a'₀+i.b'₀).z^m + (a'₁+i.b'₁).z^m-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
       Rbb = a₀ ; Rbb+1 = b₀ ; Rbb+2 = a₁ ; Rbb+3 = b₁; ....... ; Ree-1 = a_n ; Ree = b_n     these 2n+2 registers are to be initialized before executing "PROZ"
     Rbb' = a'₀ ; Rbb'+1 = b'₀ ; Rbb'+2 = a'₁ ; Rbb'+3 = b'₁ ; ....... ; Ree'-1 = a'_m ; Ree' = b'_m these 2m+2 registers are to be initialized before executing "PROZ"
Flags: /                                             ( bb > 07 & bb' > 07 )
Subroutines:  /
 

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                If you don't have an HP-41CX replace line 24  by     0  LBL 00  STO IND Y  ISG Y  GTO 00
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 = ---------------------   b(z)                           (bbb.eee)_p = the control number of the product        ( all  bbb > 007 )
 

Example:     a(z) = (2+3i).z²+(4+7i).z+(1-9i)      b(z) =  (4-6i).z³+(2-3i).z²+(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.z⁵+(71+4i).z⁴+(-32-11i).z³+(-58+49i).z²+(58-23i).z+(19-7i)
 
 

        l) Addition/Subtraction of 2 Polynomials
 

-Let    a(z) = (a₀+i.b₀).zⁿ + (a₁+i.b₁).z^n-1 + ......... + (a_n-1+i.b_n-1).z + (a_n+i.b_n)
 and    b(z) = (a'₀+i.b'₀).z^m + (a'₁+i.b'₁).z^m-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      and, with  all bb > 03

       Rbb = a₀ ; Rbb+1 = b₀ ; Rbb+2 = a₁ ; Rbb+3 = b₁ ; ....... ; Ree-1 = a_n ; Ree = b_n     these 2n+2 registers are to be initialized before executing "ADDZ"
     Rbb' = a'₀ ; Rbb'+1 = b'₀ ; Rbb'+2 = a'₁ ; Rbb'+3 = b'₁ ; ....... ; Ree'-1 = a'_m ; Ree' = b'_m these 2m+2 registers are to be initialized before executing "ADDZ"

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 = ---------------------   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).z²+(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 ...

1°)            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).z²+(4-12i).z+(4+10i)

2°)             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).z²+2i.z+(-2+2i)
 

        m) 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:
    Rbb = a₀ ; Rbb+1 = b₀ ; Rbb+2 = a₁ ; Rbb+3 = b₁ ; ....... ; Ree-1 = a_n ; Ree = b_n        these 2n+2 registers are to be initialized before executing "DELZ"
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                       ( or another "smal" number )
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).z²+(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.x²+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.x²-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.x³+b.x²+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        ( not  STOP  but synthetic  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(x₁;y₁)  and  B(x₂;y₂)  are the 2 extrema
   If you get "NON EXISTENT" there is no extremum
 

Example1:     y = 2.x³ - 3.x² - 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 = x³ - 3.x² + 3.x + 4   gives "NON EXISTENT"  ( y'(x) = 0 has only one solution )

Example3:     y = x³ + x² + x + 1  also produce "NON EXISTENT"  ( y'(x) = 0 has no real solution )
 

        c) Center of symmetry of the curve y = a.x³+b.x²+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.x³ + 3.x² + 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)
 
 
 
 

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall