hp41programs

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.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<> 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.x² + 3.x - 4 = 0 and 2.x² + 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.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

-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.x³-5.x²-7.x+1 = 0 and 2.x³-5.x²+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 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"

-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 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 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 L-register.

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 ""
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 ( Registers Rbb thru Ree are to be initialized before executing "PLR" )

• Rbb = a₀ , • Rbb+1 = a₁ , ....... , • Ree = a_n 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	x₀	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 )

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

-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	v₀	/	y₁		y_k
X	u₀	bbb.eee	x₁		x_k

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 ) -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 = 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 ( Registers Rbb thru Ree & Rbb' thru Ree' are to be initialized before executing "DIV" )

• Rbb = a₀ • Rbb+1 = a₁ , ....... , • Ree = a_n bb > 03
• Rbb' = b₀ • Rbb'+1 = b₁ , ....... , • Ree' = b_m 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.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 ( Registers Rbb thru Ree & Rbb' thru Ree' are to be initialized before executing "PRO" )

• Rbb = a₀ • Rbb+1 = a₁ , ....... , • Ree = a_n bb > 03
• Rbb' = b₀ • Rbb'+1 = b₁ , ....... , • Ree' = b_m 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	bbb_p	(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.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 ( Registers Rbb thru Ree & Rbb' thru Ree' are to be initialized before executing "ADD" )

• Rbb = a₀ • Rbb+1 = a₁ , ....... , • Ree = a_n bb > 03
• Rbb' = b₀ • Rbb'+1 = b₁ , ....... , • Ree' = b_m 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	bbb_p	(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.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 ...

      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

      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: ( Registers Rbb thru Ree are to be initialized before executing "DEL" )

• Rbb = a₀ • Rbb+1 = a₁ , ....... , • Ree = a_n
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^-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

m) Composition of 2 Polynomials

-Let f(x) = a₀.xⁿ+a₁.x^n-1+ ... + a_n-1.x+a_n and g(x) = b₀.x^m+b₁.x^m-1+ ... + b_m-1.x+b_m 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 = a₀   • Rbb+1 = a₁ , ....... , • Ree = a_n      bb > 08
                          • Rbb' = b₀ • Rbb'+1 = b₁ , ....... , • Ree' = b_m       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 x² + 3 x + 7 g(x) = 4 x³ + 5 x² + 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 x⁶ + 80 x⁵ + 146 x⁴ + 148 x³ + 107 x² + 42 x + 12

-Likewise, you'll find in 21 seconds, gof(x) = 32 x⁶ + 144 x⁵ + 572 x⁴ + 1176 x³ + 2129 x² + 1992 x + 1660

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 ( 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).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   whence   z₂ = -0.061935981 + 1.506099080.i = R03+i.R04
                       [ RCL 05 ] -1.527654077 [ RCL 06 ] 0.084669040 whence   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 ( 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).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     whence z₂ = -1.226930040 -0.762095736.i = R03+i.R04
                    [ RCL 05 ]   0.361168050   [ RCL 06 ]   1.374351077     whence z₃ = 0.361168050 + 1.374351077.i = R05+i.R06
                    [ RCL 07 ] -1.171930991 [ RCL 08 ]   0.883676309     whence 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: ( Registers Rbb thru Ree are to be initialized before executing "PLZ" )

• Rbb = a₀ • Rbb+1 = b₀ , • Rbb+2 = a₁ • Rbb+3 = b₁ , ....... , • Ree-1 = a_n • Ree = b_n
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 )

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 z³ + 5 z² + 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) = (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: ( Registers Rbb thru Ree are to be initialized before executing "dPLZ" )

• Rbb = a₀ • Rbb+1 = b₀ , • Rbb+2 = a₁ • Rbb+3 = b₁ , ....... , • Ree-1 = a_n • Ree = b_n
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) = (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: ( Registers Rbb thru Ree are to be initialized before executing "d2PLZ" )

• Rbb = a₀ • Rbb+1 = b₀ , • Rbb+2 = a₁ • Rbb+3 = b₁ , ....... , • Ree-1 = a_n • Ree = b_n
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) = (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 ( Registers Rbb thru Ree are to be initialized before executing "PLRZ" )

• Rbb = a₀ • Rbb+1 = b₀ , • Rbb+2 = a₁ • Rbb+3 = b₁ , ....... , • Ree-1 = a_n • Ree = b_n 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	y₀	/
X	x₀	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     whence   z₂ = -0.291909945 - 0.629741372.i
R13 = -0.181419655    R14 = 0.663480568     whence    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

-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) = (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 ( Registers Rbb thru Ree are to be initialized before executing "LAGZ" )

• Rbb = a₀ • Rbb+1 = b₀ , • Rbb+2 = a₁ • Rbb+3 = b₁ , ....... , • Ree-1 = a_n • Ree = b_n 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	y₀	/
X	x₀	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')

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

• Rbb = a₀ • Rbb+1 = b₀ , • Rbb+2 = a₁ • Rbb+3 = b₁ , ....... , • Ree-1 = a_n • Ree = b_n 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) = 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 ( Registers Rbb thru Ree & Rbb' thru Ree' are to be initialized before executing "DIVZ" )

• Rbb = a₀ • Rbb+1 = b₀ , • Rbb+2 = a₁ • Rbb+3 = b₁ , ....... , • Ree-1 = a_n • Ree = b_n bb > 07
• Rbb' = a'₀ • Rbb'+1 = b'₀ , • Rbb'+2 = a'₁ • Rbb'+3 = b'₁ , ....... , • 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).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

-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

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

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

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

"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 = a₀ • Rbb+1 = b₀ , • Rbb+2 = a₁ • Rbb+3 = b₁ , ....... , • Ree-1 = a_n • Ree = b_n bb > 03
• Rbb' = a'₀ • Rbb'+1 = b'₀ , • Rbb'+2 = a'₁ • Rbb'+3 = b'₁ , ....... , • 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	bbb_p	(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).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 ...

       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)

       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 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 = a₀ • Rbb+1 = b₀ , • Rbb+2 = a₁ • Rbb+3 = b₁ , ....... , • Ree-1 = a_n • Ree = b_n
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).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
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	y₂
Z	b	x₂
Y	c	y₁
X	d	x₁

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