Solving for the roots of a cubic polynomial

[Namir]

This is a Free42 printout of an HP-41C program that uses the instructions found on page 75 of the HP-4 Application Book. The program calculates the real and complex roots of a cubic polynomial:

x^3 + a*x^2 + b^x + c = 0

The program prompts you for the coefficients, b, a, and c in that order. The HP-41C adaptation uses labels and flags 00 and 01 to direct the flow of the calculations. The output includes the determinant d and the roots x1, x2, and x3. When the roots are complex the output appears as <real,imaginary>.

The listing (revised from HP-42s to HP-41C) is:

Code:

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

[Albert Chan]

Hi, Namir

x^3 + x^2 = 0, but did not get back roots -1, 0, 0

Instead, I get d=0, X1 = X2 = X3 = -1/3

Try its depressed cubic, with x = y - 1/3

y^3 - y/3 + 2/27 = 0

Sure enough, I get d=0, X1 = X2 = X3 = 0

[Namir]

(06-30-2024 11:49 PM)Albert Chan Wrote:  Hi, Namir

x^3 + x^2 = 0, but did not get back roots -1, 0, 0

Instead, I get d=0, X1 = X2 = X3 = -1/3

Try its depressed cubic, with x = y - 1/3

y^3 - y/3 + 2/27 = 0

Sure enough, I get d=0, X1 = X2 = X3 = 0

I think the HP author of the program assumed that parameter c should not be 0, because in this case the cubic equation s easily reduced to a quadratic or linear equation (as is your example).

[Albert Chan]

(07-01-2024 04:13 AM)Namir Wrote:  

(06-30-2024 11:49 PM)Albert Chan Wrote:  Try its depressed cubic, with x = y - 1/3

y^3 - y/3 + 2/27 = 0

Sure enough, I get d=0, X1 = X2 = X3 = 0

I think the HP author of the program assumed that parameter c should not be 0, because in this case the cubic equation s easily reduced to a quadratic or linear equation (as is your example).

Depressed cubic does not look easily reduced to quadratic.
Fix should not be hard. Here is the trig solution.

[Namir]

As I mentioned before the original HP-45 key sequence that I translated into HP-41C code seem to have been meant to solve a cubic equation x^3+a*x^2+b*x+c =0, for typical cases where a, b, and c are not zero. You are welcome to enhance the HP-41C code for various combinations of zeros in these three parameters. You can download the HP-45 Application Book from Eric Rechlin's documents web pages and look at pages 74 and 75 to see the equations used and the keystroke sequence. An alternative to all this is to use Lin-Bairstow's algorithm that solves for real and complex roots of polynomials with real coefficients.

[Albert Chan]

Hi, Namir

It is not zeros that caused the problem. You could even try x^3 = 0, and it still work.

The problem is cube root of negative numbers. (note: all coefficients are real)
Code seems to assume first cube root is always positive.
2nd cube root use flag 01 to get sign(x) * ³√|x|

This is why y^3 - y/3 + 2/27 = 0 roots are 0

y1 = ³√(-1/27+0) + ³√(-1/27−0), mistakenly calculated as 1/3 - 1/3 = 0

We could reuse flag 01 for both cube roots, or do without. ³√(x) = sign(x) * ³√|x|

Note: I only touched cube roots. (LINE 65 .. 71, copy/paste to LINE 76 .. 82)

Code:

00 { 281-Byte Prgm }
01▸LBL "CUBE"
02▸LBL A
03 CF 00
04 "b?"
05 PROMPT
06 STO 07
07 "a?"
08 PROMPT
09 STO 08
10 ENTER
11 ×
12 3
13 ÷
14 -
15 STO 01
16 RCL 08
17 3
18 ÷
19 ENTER
20 ENTER
21 ENTER
22 ×
23 ×
24 2
25 ×
26 X<>Y
27 RCL 07
28 ×
29 -
30 "c?"
31 PROMPT
32 +
33 STO 02
34 ENTER
35 ×
36 4
37 ÷
38 RCL 01
39 ENTER
40 ENTER
41 ×
42 ×
43 27
44 ÷
45 +
46 RND
47 STO 03
48 "d="
49 ARCL ST X
50 PROMPT
51▸LBL 11
52 DEG
53 X=0?
54 SF 00
55 X<0?
56 GTO 31
57 SQRT
58 STO 03
59 RCL 02
60 +/-
61 2
62 ÷
63 STO 02
64 +
65 SIGN
66 LASTX
67 ABS
68 3
69 1/X
70 Y↑X
71 ×
72 STO 04
73 RCL 02
74 RCL 03
75 -
76 SIGN
77 LASTX
78 ABS
79 3
80 1/X
81 Y↑X
82 ×
83 STO 05
84 RCL 04
85 +
86 STO 04
87 RCL 08
88 3
89 ÷
90 STO 08
91 -
92 "X1="
93 ARCL ST X
94 PROMPT
95 RCL 04
96 2
97 ÷
98 +/-
99 STO 06
100 RCL 08
101 -
102 STO 11
103 FS?C 00
104 GTO 00
105 RCL 04
106 RCL 05
107 2
108 ×
109 -
110 2
111 ÷
112 3
113 SQRT
114 ×
115 STO 05
116 STO 12
117 CLA
118 ARCL 11
119 ├","
120 ARCL 12
121 PROMPT
122 CLA
123 RCL 12
124 +/-
125 ARCL 11
126 ├","
127 ARCL ST X
128 PROMPT
129 RTN
130▸LBL 31
131 RCL 02
132 RCL 01
133 ENTER
134 ENTER
135 ×
136 ×
137 +/-
138 27
139 ÷
140 SQRT
141 2
142 ×
143 ÷
144 +/-
145 ACOS
146 3
147 ÷
148 STO 02
149 COS
150 RCL 01
151 3
152 ÷
153 X↑2
154 SQRT
155 SQRT
156 2
157 ×
158 STO 01
159 ×
160 RCL 08
161 3
162 ÷
163 STO 08
164 -
165 "X1="
166 ARCL ST X
167 PROMPT
168 RCL 02
169 120
170 +
171 COS
172 RCL 01
173 ×
174 RCL 08
175 -
176 "X2="
177 ARCL ST X
178 PROMPT
179 RCL 02
180 240
181 +
182 COS
183 RCL 01
184 ×
185 RCL 08
186 -
187 "X3="
188 ARCL ST X
189 PROMPT
190 RTN
191▸LBL 00
192 "X2="
193 ARCL 11
194 PROMPT
195 "X3="
196 ARCL 11
197 PROMPT
198 RTN
199 END

[Namir]

Thanks Albert for the updated code.

Namir

[floppy]

What is the mathematical background? like Page 8 of https://arxiv.org/pdf/2206.03855 ?