Polynomials2

Discriminant for the HP-41

Overview

1°) Polynomials of degree 3
2°) Polynomials of degree 4
3°) Polynomials of degree n

1°) Polynomials of Degree 3

-The discriminant of a polynomial of degree 2: a x² + b x + c is easy to compute: D = b² - 4 a c
-With a polynomial of degree 3: a x³ + b x² + c x + d , the formula is already more complicated:

D = b² c² - 4 a c³ - 4 b³ d - 27 a² d² + 18 a b c d

Data Registers: R00 unused

R01 = a R02 = b R03 = c R04 = d
Flags: /
Subroutines: /

01 LBL "DIS3"
02 STO 04
03 X<>Y
04 STO 03
05 *
06 X<>Y
07 STO 02
08 *

09 X<>Y
10 STO 01
11 *
12 ST+ X
13 RCL 01
14 RCL 04
15 *
16 X^2

17 3
18 *
19 -
20 9
21 *
22 RCL 02
23 3
24 Y^X

25 RCL 04
26 *
27 RCL 03
28 3
29 Y^X
30 RCL 01
31 *
32 +

33 4
34 *
35 -
36 RCL 02
37 RCL 03
38 *
39 X^2
40 +
41 END

( 51 bytes / SIZE 005 )

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

Examples:

  • 2 x³ - 5 x² - 7 x + 1

   2 ENTER^
-5 ENTER^
-7 ENTER^
1 XEQ "DIS3" >>>> D = 5621 ---Execution time = 2s---

D > 0 so there are 3 distinct roots

• 2 x³ - 5 x² + 7 x + 1

   2 ENTER^
-5 ENTER^
7 ENTER^
1 R/S >>>> D = -2387

D < 0 so there is one real root

• x³ + x² - 5 x + 3

   1 ENTER^ ENTER^
-5 ENTER^
3 R/S >>>> D = 0

D = 0 so there are ( at least ) 2 equal roots

2°) Polynomials of Degree 4

-Of course, with a polynomial of degree 4: a x⁴ + b x³ + c x² + d x + e , the formula is much more complicated:

D = 256 a³ e³ -192 a² b d e² - 128 a² c² e² + 144 a² c d² e - 27 a² d⁴ + 144 a b² c e² - 6 a b² d² e
- 80 a b c² d e + 18 a b c d³ + 16 a c⁴ e - 4 a c³ d² - 27 b⁴ e² + 18 b³ c d e - 4 b³ d³ - 4 b² c³ e + b² c² d²

Data Registers: R00: temp ( Registers R01 thru R05 are to be initialized before executing "DIS4" )

• R01 = a • R02 = b • R03 = c • R04 = d • R05 = e
Flags: /
Subroutines: /

01 LBL "DIS4"
02 RCL 01
03 RCL 05
04 *
05 STO 00
06 16
07 *
08 RCL 02
09 RCL 04
10 *
11 12
12 *
13 -
14 RCL 03
15 X^2
16 8
17 *
18 -
19 RCL 00
20 X^2
21 *
22 RCL 00
23 RCL 03

24 X^2
25 X^2
26 *
27 +
28 RCL 00
29 RCL 02
30 *
31 RCL 03
32 X^2
33 *
34 RCL 04
35 *
36 5
37 *
38 -
39 16
40 *
41 RCL 01
42 RCL 04
43 *
44 X^2
45 RCL 03
46 *

47 RCL 05
48 *
49 RCL 02
50 LASTX
51 *
52 X^2
53 RCL 01
54 *
55 RCL 03
56 *
57 +
58 144
59 *
60 +
61 RCL 01
62 RCL 04
63 X^2
64 *
65 X^2
66 RCL 02
67 X^2
68 RCL 05
69 *

70 X^2
71 +
72 27
73 *
74 -
75 RCL 02
76 RCL 04
77 *
78 X^2
79 RCL 00
80 *
81 6
82 *
83 -
84 RCL 01
85 RCL 02
86 *
87 RCL 03
88 *
89 RCL 04
90 3
91 Y^X
92 *

93 RCL 02
94 3
95 Y^X
96 RCL 03
97 *
98 RCL 04
99 *
100 RCL 05
101 *
102 +
103 18
104 *
105 +
106 RCL 02
107 RCL 04
108 *
109 3
110 Y^X
111 RCL 03
112 LASTX
113 Y^X
114 STO 00
115 RCL 02

116 X^2
117 *
118 RCL 05
119 *
120 +
121 RCL 00
122 RCL 01
123 *
124 RCL 04
125 X^2
126 *
127 +
128 4
129 *
130 -
131 RCL 02
132 RCL 03
133 *
134 RCL 04
135 *
136 X^2
137 +
138 END

( 154 bytes / SIZE 005 )

STACK	INPUT	OUTPUT
X	/	D

Examples:

  • x⁴ + 3 x³ + x² - 3 x - 2

1 STO 01 STO 03
3 STO 02 CHS STO 04
-2 STO 05

XEQ "DIS4" >>>> D = 0 ( at least 2 equal roots ) ---Execution time = 6s---

  • x⁴ - 2 x³ - 13 x² + 38 x - 24

1 STO 01
-2 STO 02
-13 STO 03
38 STO 04
-24 STO 05 R/S >>>> D = 176400 ( All roots are real or all roots are non-real )

  • x⁴ - 13 x³ + 60 x² - 112 x + 64

1    STO 01
-13 STO 02
60 STO 03
-112 STO 04
64 STO 05 R/S >>>> D = 10

-Unfortunately, this result is wrong because it involves numbers > 10¹⁰
-So, roundoff-errors yield 10 wheras the exact result is D = 0
-However, this program gives the correct discriminant if we use free42 !

-Thus, with an HP41, it's better to employ the following routine which calculates D with a determinant of a square matrix.

3°) Polynomials of Degree n

-If p(x) = a_n xⁿ + a_n-1 x^n-1 + ..... + a₁ x + a₀ is a polynomial of degree n > 1, the discriminant may be computed by the following formula:

D = (-1)^n(n-1)/2 (1/a_n) Det (M) where M is the square matrix of order 2n - 1

a_n 0 0 ..... b_n-1 0 0 ........
a_n-1 a_n 0 ..... b_n-2 b_n-1 0 ........ where the 1st derivative of p(x) is p'(x) = b_n-1 x^n-1 + ..... + b₁ x + b₀
a_n-2 a_n-1 a_n ..... b_n-3 b_n-2 b_n-1
...............................................................

0 ......... 0 ......a₀ 0 0 0 ... b₀

with (n-1) columns and n columns
with p(x) with p'(x)
coefficients coefficients
and zeros and zeros

-For instance, if p(x) = x⁴ - 13 x³ + 60 x² - 112 x + 64 we have to calculate the determinant of the 7x7 matrix:

1 0 0 4    0 0 0
-13 1 0 -39 4 0    0
60 -13 1 120 -39 4    0
-112 60 -13 -112 120 -39 4
64 -112 60 0 -112 120 -39
0 64 -112 0 0 -112 120
0 0 64 0 0 0 -112

Data Registers: R00: temp          ( Registers R01 thru Rn+1 are to be initialized before executing "DISC" )

• R01 = a_n • R02 = a_n-1 ...... • Rn+1 = a₀

Flags:    CF 00 = a pivot p is regarded as zero if | p | < 10^-7                   CF 01 = partial pivoting
               SF 00 =  a pivot p is regarded as zero if    p = 0                         SF 01 = no pivoting

Subroutines: "DET" ( cf "Determinant for the HP41" §4a) ) after adding RCL M just after LBL 08

01 LBL "DISC"
02 STO 00
03 STO N
04 ST+ X
05 DSE X
06 STO M
07 X^2
08 STO O
09 LBL 01
10 RCL 00
11 1
12 SIGN
13 LBL 02
14 CLX
15 RCL IND Y
16 LASTX
17 *
18 STO IND O

19 DSE O
20 ISG L
21 CLX
22 DSE Y
23 GTO 02
24 DSE N
25 FS? 30
26 GTO 00
27 RCL 00
28 SIGN
29 CLX
30 LBL 03
31 STO IND O
32 DSE O
33 DSE L
34 GTO 03
35 GTO 01
36 LBL 00

37 RCL 00
38 DSE X
39 STO N
40 LBL 04
41 RCL 00
42 1
43 +
44 SIGN
45 DSE N
46 X=0?
47 GTO 00
48 LBL 05
49 RCL IND L
50 STO IND O
51 DSE O
52 DSE L
53 GTO 05
54 RCL 00

55 DSE X
56 SIGN
57 CLX
58 LBL 06
59 STO IND O
60 DSE O
61 DSE L
62 GTO 06
63 GTO 04
64 LBL 00
65 RCL M
66 X^2
67 1
68 ST+ Y
69 CHS
70 RCL 00
71 X^2
72 LASTX
73 -

74 2
75 /
76 Y^X
77 RCL 01
78 *
79 STO IND Y
80 X<> M
81 XEQ "DET"
82 X<>Y
83 FRC
84 E3
85 *
86 INT
87 1
88 +
89 X<>Y
90 RCL IND Y
91 /
92 END

( 145 bytes / SIZE (2n-1)²+2 )

STACK	INPUT	OUTPUT
X	n	D

Examples:

  • 2 x³ - 5 x² - 7 x + 1

CF 00 CF 01

   2 STO 01
-5 STO 02
-7 STO 03
1 STO 04

3 XEQ "DISC" >>>> D = 5621 ---Execution time = 36s---

  • x⁴ - 13 x³ + 60 x² - 112 x + 64

CF 00 CF 01

1    STO 01
-13 STO 02
60 STO 03
-112 STO 04
64 STO 05

4 XEQ "DISC" >>>> D = 0 ---Execution time = 75s---

Notes:

-We can also use a program that solves linear systems.
-For instance, with "LS3" listed in "Linear & Non-Linear Systems for the HP41" §1c) , replace lines 81 to 92 with

81 ENTER^
82 XEQ "LS3"
83 R^
84 INT
85 X^2
86 1
87 +
88 X<>Y
89 RCL IND Y
90 /
91 END

-With the examples above, execution times become 48 & 92 seconds.
-We can also insert "DET" in the same program: it saves one data register:

Data Registers: R00: temp          ( Registers R01 thru Rn+1 are to be initialized before executing "DISC" )

• R01 = a_n • R02 = a_n-1 ...... • Rn+1 = a₀

Flags:    CF 00 = a pivot p is regarded as zero if | p | < 10^-7
               SF 00 =  a pivot p is regarded as zero if    p = 0

Subroutines: /

-Line 176 is a three-byte GTO 07

01 LBL "DISC"
02 STO 00
03 STO N
04 ST+ X
05 DSE X
06 STO M
07 X^2
08 STO O
09 LBL 01
10 RCL 00
11 1
12 SIGN
13 LBL 02
14 CLX
15 RCL IND Y
16 LASTX
17 *
18 STO IND O
19 DSE O
20 ISG L
21 CLX
22 DSE Y
23 GTO 02
24 DSE N
25 FS? 30
26 GTO 00

27 RCL 00
28 SIGN
29 CLX
30 LBL 03
31 STO IND O
32 DSE O
33 DSE L
34 GTO 03
35 GTO 01
36 LBL 00
37 RCL 00
38 DSE X
39 STO N
40 LBL 04
41 RCL 00
42 1
43 +
44 SIGN
45 DSE N
46 X=0?
47 GTO 00
48 LBL 05
49 RCL IND L
50 STO IND O
51 DSE O
52 DSE L

53 GTO 05
54 RCL 00
55 DSE X
56 SIGN
57 CLX
58 LBL 06
59 STO IND O
60 DSE O
61 DSE L
62 GTO 06
63 GTO 04
64 LBL 00
65 1
66 CHS
67 RCL 00
68 X^2
69 LASTX
70 -
71 2
72 /
73 Y^X
74 RCL 01
75 /
76 STO 00
77 X<> M
78 STO Y

79 .1
80 %
81 +
82 STO N
83 FRC
84 STO O
85 *
86 1
87 ST+ O
88 LASTX
89 %
90 +
91 +
92 LBL 07
93 STO M
94 INT
95 RCL O
96 FRC
97 +
98 ENTER
99 ENTER
100 CLX
101 LBL 08
102 RCL IND Z
103 ABS
104 X>Y?

105 STO Z
106 X>Y?
107 +
108 RDN
109 ISG Z
110 GTO 08
111 RCL M
112 ENTER
113 FRC
114 R^
115 INT
116 +
117 X=Y?
118 GTO 00
119 LBL 09
120 RCL IND X
121 X<> IND Z
122 STO IND Y
123 ISG Y
124 RDN
125 ISG Y
126 GTO 09
127 RCL 00
128 CHS
129 STO 00
130 LBL 00

131 CLX
132 FC? 00
133 E-7
134 RCL IND M
135 ST* 00
136 ABS
137 X<=Y?
138 GTO 12
139 ISG O
140 X<0?
141 GTO 13
142 RCL O
143 STO P
144 LBL 10
145 RCL M
146 ENTER
147 FRC
148 RCL P
149 INT
150 +
151 RCL IND X
152 RCL IND Z
153 /
154 SIGN
155 ISG Y
156 CLX

157 ISG Z
158 LBL 11
159 CLX
160 RCL IND Z
161 LASTX
162 *
163 ST- IND Y
164 ISG Y
165 CLX
166 ISG Z
167 GTO 11
168 ISG P
169 GTO 10
170 RCL N
171 ST+ O
172 SIGN
173 RCL M
174 +
175 ISG X
176 GTO 07
177 LBL 12
178 CLX
179 STO 00
180 LBL 13
181 RCL 00
182 CLA
183 END

( 275 bytes / SIZE (2n-1)²+1 )

STACK	INPUT	OUTPUT
X	n	D

Examples:

• 2 x³ - 5 x² - 7 x + 1

CF 00

   2 STO 01
-5 STO 02
-7 STO 03
1 STO 04

3 XEQ "DISC" >>>> D = 5621 ---Execution time = 35s---

  • x⁴ - 13 x³ + 60 x² - 112 x + 64

1    STO 01
-13 STO 02
60 STO 03
-112 STO 04
64 STO 05

4 R/S >>>> D = 0    ---Execution time = 74s---

• 2 x⁷+ 3 x⁶ + 4 x⁵ + 5 x⁴ + 6 x³ + 7 x² + 8 x + 9

2 STO 01 3 STO 02 4 STO 03 5 STO 04 6 STO 05 7 STO 06 8 STO 07 9 STO 08

7 R/S >>>> D = -6.302249843 10¹² ---Execution time = 5mn59s---

Notes:

-Since there are at most 319 data registers, the maximum degree is 9
-But we can modify this program if you want to use it with free42:

Data Registers: R00.... R09: temp R10 = D ( Registers R11 thru Rn+11 are to be initialized before executing "DISC" )

• R11 = a_n • R12 = a_n-1 ...... • Rn+11 = a₀

Flags: /
Subroutines: /

-Line 207 is a three-byte GTO 07

01 LBL "DISC"
02 STO 00
03 STO 01
04 ST+ X
05 DSE X
06 STO 04
07 X^2
08 10
09 +
10 STO 02
11 LBL 01
12 RCL 00
13 10.01
14 +
15 STO 03
16 1
17 SIGN
18 LBL 02
19 CLX
20 RCL IND Y
21 LASTX
22 *
23 STO IND 02
24 DSE 02
25 ISG L
26 CLX
27 DSE Y
28 GTO 02
29 DSE 01
30 FS? 30

31 GTO 00
32 RCL 03
33 SIGN
34 CLX
35 LBL 03
36 STO IND 02
37 DSE 02
38 DSE L
39 GTO 03
40 GTO 01
41 LBL 00
42 RCL 00
43 DSE X
44 STO 01
45 LBL 04
46 RCL 03
47 1
48 +
49 SIGN
50 DSE 01
51 X=0?
52 GTO 00
53 LBL 05
54 RCL IND L
55 STO IND 02
56 DSE 02
57 DSE L
58 GTO 05
59 RCL 00
60 DSE X

61 SIGN
62 CLX
63 LBL 06
64 STO IND 02
65 DSE 02
66 DSE L
67 GTO 06
68 GTO 04
69 LBL 00
70 1
71 CHS
72 RCL 00
73 X^2
74 LASTX
75 -
76 2
77 /
78 Y^X
79 RCL 11
80 /
81 STO 10
82 RCL 04
83 STO 02
84 STO 04
85 STO 08
86 X^2
87 11
88 STO 09
89 +
90 1

91 -
92 STO 03
93 RCL 04
94 E5
95 /
96 +
97 LBL 07
98 STO 01
99 INT
100 ENTER
101 ENTER
102 RCL 04
103 STO 05
104 CLX
105 LBL 08
106 RCL IND Z
107 ABS
108 X>Y?
109 STO Z
110 X>Y?
111 +
112 RDN
113 DSE Z
114 STO X
115 DSE 05
116 GTO 08
117 RCL 01
118 ENTER
119 FRC
120 R^

121 INT
122 +
123 X=Y?
124 GTO 00
125 RCL 08
126 STO 05
127 RDN
128 LBL 09
129 RCL IND X
130 X<> IND Z
131 STO IND Y
132 DSE Y
133 RDN
134 DSE Y
135 STO X
136 DSE 05
137 GTO 09
138 RCL 10
139 CHS
140 STO 10
141 LBL 00
142 E-7
143 RCL IND 01
144 ST* 10
145 ABS
146 X<=Y?
147 CLX
148 X=0?
149 STO 10
150 X=0?

151 GTO 00
152 RCL 08
153 STO 05
154 RCL 01
155 LASTX
156 LBL 10
157 ST/ IND Y
158 DSE Y
159 STO X
160 DSE 05
161 GTO 10
162 RCL 02
163 STO 05
164 RCL 03
165 STO 06
166 LBL 11
167 RCL 01
168 ENTER
169 FRC
170 RCL 06
171 +
172 X=Y?
173 GTO 13
174 RCL 08
175 STO 07
176 CLX
177 RCL IND Y
178 SIGN
179 RDN
180 LBL 12

181 RCL IND Y
182 LASTX
183 *
184 ST- IND Y
185 DSE Y
186 RDN
187 DSE Y
188 STO X
189 DSE 07
190 GTO 12
191 LBL 13
192 DSE 06
193 DSE 05
194 GTO 11
195 LBL 00
196 RCL 02
197 ST- 03
198 RCL 01
199 1
200 -
201 DSE X
202 STO X
203 DSE 08
204 FS? 30
205 GTO 14
206 DSE 04
207 GTO 07
208 LBL 14
209 RCL 10
210 END

( 302 bytes / SIZE (2n-1)²+11 )

STACK	INPUT	OUTPUT
X	n	D = R10

Examples: With free42, line 142 should be replaced by a smaller number ( perhaps E-20 )

• 2 x⁷+ 3 x⁶ + 4 x⁵ + 5 x⁴ + 6 x³ + 7 x² + 8 x + 9

2 STO 11 3 STO 12 4 STO 13 5 STO 14 6 STO 15 7 STO 16 8 STO 17 9 STO 18

7 XEQ "DISC" >>>> D = -6302249844736

• 2 x¹⁰+ 3 x⁹ + 4 x⁸ + 5 x⁷ + 6 x⁶ + 7 x⁵ + 8 x⁴ + 9 x³ + 10 x² + 11 x + 12

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

10 R/S >>>> D = -3727403657953361647360

Notes:

-I've rounded the results above: in the 2nd example, X-output = -3727403657953361647359,999999999999
-With free42, we can can compute the discriminant of polynomials of degree much larger than 9...

Reference:

[1] https://en.wikipedia.org/w/index.php?title=Discriminant&oldid=1057323016

hp41programs

Discriminant for the HP-41