.

Polynomial Factorization (to 5th degree) for HP-41cv/cx

This program calculates the roots of a polynomial up to the 5th degree (quintic).\(\require{cancel}\)

 

History

Around 1985, Michiel Niemeijer published a similar program for the HP-41 in the Dutch publication HP User News [1].  I am not sure if that was the source.  My changes were aimed at making the user experience more intuitive. This article lists my version 1.1 from October 12, 1986.

Now, twenty-five years later, I restored my HP-41CV and use the program once again. The HP-41 version is also verified to run under Windows using Warren Furlow’s V41 simulator when compiled with Leo Duran’s hp41uc.

Theory

Consider the monic fourth-order polynomial [3]

$$p(x)=x^4+x^3+x^2+11x+10$$

A root of polynomial \(p(x)\) is a number \(x_i\) such that \(p(x_i)=0\) .  Each nth-order polynomial has exactly n roots.  These roots may be real or complex.  Then we have \(p(x_i)=0\), and we can write

$$p(x)=(x-x_1)(x-x_2)(x-x_3)(x-x_4)$$

This is the factored form of the monic polynomial.  Multiplying out the symbolic factored form gives a nonlinear system of four equations and four unknowns.  Solving it gives

$$\begin{cases}
x_1=-1\\
x_2=-2\\
x_3=1+2j\\
x_4=1-2j
\end{cases}$$

Quadratic polynomial (2nd order)

The generic form is

$$\begin{array}{cr}
f(x)=ax^2+bx+c,&a\neq0
\end{array}$$

Normalize and find the roots where \(f(x)=0\)

$$\begin{array}{ll}
&ax^2+bx+c=0\\
\Rightarrow &x^2+\frac{b}{a}x+\frac{c}{a}=0\\
\end{array}$$

By completing the square, we can use the algebra identity \((x+q)^2=x^2+2qx+q^2\) to find the roots [mathisfun]

$$\begin{array}{ll}
&\left(x+\frac{b}{2a}\right)^2-\left(\frac{b}{2a}\right)^2+\frac{c}{a}=0\\
\Rightarrow &x+\frac{b}{2a}=\pm\sqrt{\left(\frac{b}{2a}\right)^2-\frac{c}{a}}\\
\Rightarrow &x=\frac{-b\pm\sqrt{b^2-4ac}}{2a}\\
\end{array}\label{eq:quadratic}\tag{1}$$

Cubic polynomial (3rd order)

In 1545, Gerolamo Cardano published the method for finding the roots of [wiki]

$$\begin{array}{cr}
f(x)=ax^3+bx^2+cx+d,&a\neq0
\end{array}$$

In 1557, Tartaglia showed that the \(x^2\) term can be eliminated by substituting \(t-\frac{b}{3a}\) for \(x\) and applying the algebra identifies \((a-b)^3=a^3-3a^2b+3ab^2-b^3\) and \((a-b)^2=a^2-2ab+b^2\)

$$\begin{array}{ll}
&a\left(t-\frac{b}{3a}\right)^3+b\left(t-\frac{b}{3a}\right)^2+c\left(t-\frac{b}{3a}\right)+d=0\\
\Rightarrow& a\left[
t^3
– 3t^2\left(\frac{b}{3a}\right)
+ 3t\left(\frac{b}{3a}\right)^2
– \left(\frac{b}{3a}\right)^3
\right]
+ b\left[
t^2-2t\left(\frac{b}{3a}\right)+\left(\frac{b}{3a}\right)^2
\right]
+ c\left[t-\frac{b}{3a}\right]
+ d=0\\
\Rightarrow &
at^3
– bt^2
+ \frac{b^2}{3a}t
– \frac{b^3}{27a^3}
+ bt^2
– \frac{2b^2}{3a}t + \frac{b^3}{9a^2}
+ ct
-\frac{bc}{3a}
+ d=0\\
\Rightarrow&
at^3
+\left(b-b\right)t^2
+\left(
\frac{b^2}{3a}
– \frac{2b^2}{3a}
+c
\right)t
+\left(
d
+ \frac{2b^3}{27a^2}
-\frac{bc}{3a}
\right)=0\\
\Rightarrow&
at^3
+\left( c -\frac{b^2}{3a} \right) t
+\left( d + \frac{2b^3}{27a^2} -\frac{bc}{3a} \right)=0\\
\end{array}$$

That leaves us with a depressed cubic in the form \(t^3+pt+q=0\)

$$\begin{array}{lll}
t^3+pt+q=0,&
p=\frac{c-\frac{b^2}{3a}}{a}, &
q=\frac{d + \frac{2b^3}{27a^2} -\frac{bc}{3a}}{a}\\
t^3+pt+q=0,&
p=\frac{3ac-b^2}{3a^2}, &
q=\frac{27a^2d + 2b^3 -9abc}{27a^3}\\
\end{array}$$

Solve the depressed cubic

Starting with this depressed cubic

$$t^3+pt+q=0\label{eq:cubic_depressed}\tag{2}$$

Use the Vieta substitution

$$t=w-\frac{p}{3w}\label{eq:cubic_vieta}\tag{3}$$

to solve the depressed cubic

$$\begin{array}{lll}
&\begin{array}{lll}
t^3+pt+q=0,&
t=w-\frac{p}{3w}\\
\end{array}\\
\Rightarrow &\left(w-\frac{p}{3w}\right)^3
+ p\left(w-\frac{p}{3w}\right)
+ q = 0\\
\Rightarrow&w^3
-3w^2\frac{p}{3w}
+3w\frac{p^2}{3^2w^2}
-\frac{p^3}{3^3w^3}
+pw
-\frac{p^2}{3w}
+q = 0 \\
\Rightarrow&w^3
\cancel{-pw}
\bcancel{+\frac{p^2}{3w}}
-\frac{p^3}{3^3w^3}
\cancel{+pw}
\bcancel{-\frac{p^2}{3w}}
+q = 0\\
\Rightarrow&w^3
-\frac{p^3}{27w^3}
+q = 0\\
\Rightarrow&w^6
+qw^3
-\frac{p^3}{27}=0
\end{array}$$

Solve this \(w^3\) using sagemath (the quadratic formula \(\eqref{eq:quadratic}\))

w = var(‘w’)
p = var(‘p’)
q = var(‘q’)
for s in solve(w^3 -q/2 + sqrt(q^2/4 – p^3/27), w):
show(s)
#print latex(s)
for s in solve(w^3 -q/2 – sqrt(q^2/4 – p^3/27), w):
show(s)
#print latex(s)

yields

$$\begin{array}{lll}
&w^3=-\frac{q}{2}\pm\sqrt{\frac{q^2}{4}-\frac{p^3}{27}}\\
\Rightarrow &\begin{cases}
w_1 = {\left(\frac{1}{2} \, q + \sqrt{-\frac{1}{27} \, p^{3} + \frac{1}{4} \, q^{2}}\right)}^{\frac{1}{3}}\\
w_2 = {\left(\frac{1}{2} \, q – \sqrt{-\frac{1}{27} \, p^{3} + \frac{1}{4} \, q^{2}}\right)}^{\frac{1}{3}}\\
w_3 = \frac{1}{2} i \, \sqrt{3} {\left(\frac{1}{2} \, q – \sqrt{-\frac{1}{27} \, p^{3} + \frac{1}{4} \, q^{2}}\right)}^{\frac{1}{3}} – \frac{1}{2} \, {\left(\frac{1}{2} \, q – \sqrt{-\frac{1}{27} \, p^{3} + \frac{1}{4} \, q^{2}}\right)}^{\frac{1}{3}}\\
w_4 = \frac{1}{2} i \, \sqrt{3} {\left(\frac{1}{2} \, q + \sqrt{-\frac{1}{27} \, p^{3} + \frac{1}{4} \, q^{2}}\right)}^{\frac{1}{3}} – \frac{1}{2} \, {\left(\frac{1}{2} \, q + \sqrt{-\frac{1}{27} \, p^{3} + \frac{1}{4} \, q^{2}}\right)}^{\frac{1}{3}}\\
w_5 = -\frac{1}{2} i \, \sqrt{3} {\left(\frac{1}{2} \, q – \sqrt{-\frac{1}{27} \, p^{3} + \frac{1}{4} \, q^{2}}\right)}^{\frac{1}{3}} – \frac{1}{2} \, {\left(\frac{1}{2} \, q – \sqrt{-\frac{1}{27} \, p^{3} + \frac{1}{4} \, q^{2}}\right)}^{\frac{1}{3}}\\
w_6 = -\frac{1}{2} i \, \sqrt{3} {\left(\frac{1}{2} \, q + \sqrt{-\frac{1}{27} \, p^{3} + \frac{1}{4} \, q^{2}}\right)}^{\frac{1}{3}} – \frac{1}{2} \, {\left(\frac{1}{2} \, q + \sqrt{-\frac{1}{27} \, p^{3} + \frac{1}{4} \, q^{2}}\right)}^{\frac{1}{3}}\\
\end{cases}
\end{array}$$

Take a real root

$$\begin{array}{lll}
w_1 = \sqrt[3]{\frac{1}{2} \, q + \sqrt{-\frac{1}{27} \, p^{3} + \frac{1}{4} \, q^{2}}}\\
\end{array}$$

Put this back in the substitution \(\eqref{eq:cubic_vieta}\)

$$t_1=w_1-\frac{p}{3w_1}$$

Put this back in the substitution \(x=t-\frac{b}{3a}\) yields one of the roots

$$\begin{array}{lll}
x_1=w_1-\frac{p}{3w_1}-\frac{b}{3a},&
w_1 = \sqrt[3]{\frac{1}{2} \, q + \sqrt{-\frac{1}{27} \, p^{3} + \frac{1}{4} \, q^{2}}},&
p=\frac{3ac-b^2}{3a^2}, &
q=\frac{27a^2d + 2b^3 -9abc}{27a^3}
\end{array}$$

Polynomial division removes the root, and leaves a quadratic polynomial

For example

Solve

$$2x^3-30x^2+162x-350=0$$

 

 

((All at once))

$$\begin{array}{llll}
x=\sqrt[3]{q+\sqrt{q^2+(r-p^2)^3}},&
p=-\frac{b}{3a},&
q=p^3+\frac{bc-3ad}{6a^2},&
r=\frac{c}{3a}
\end{array}$$

Quartic polynomial (4th order)

Quintic polynomial (5th order)

Can only be solved numerical

Examples

Assumes FIX 3 and SIZE 15.

  1. Calculate the roots for: x5 – x4 – 101.x3 + 101.x2 + 100.x – 100 = 0
    entry display
    XEQ “POL5” ORDER?
    5 R/S X^5= ?
    1 R/S X^4= ?
    1 CHS
    R/S X^3= ?
    101 CHS
    R/S X^2= ?
    101 R/S X^1=?
    100 R/S CONST?
    100 CHS
    R/S 1.000
    R/S -1.000
    R/S 10.000
    R/S 1.000
    R/S -10.000
    R/S END
  2. Calculate the roots for: x4 + x3 + x2 + 11.x + 10 = 0
    entry display
    R/S ORDER?
    4 R/S X^4= ?
    1 R/S X^3= ?
    1 R/S X^2= ?
    1 R/S X^1= ?
    11 R/S CONST= ?
    10 R/S -1.000
    R/S -2.000
    R/S 1.000+2.000J
    R/S 1.000-2.000J
    R/S END
  3. Calculate the roots for: x3 – x2 + x – 1 = 0
    entry display
    XEQ “POL5” ORDER?
    3 R/S X^3= ?
    1 R/S X^2= ?
    1 CHS
    R/S X^1= ?
    1 R/S CONST?
    1 CHS
    R/S 1.000E-10+1.000J
    R/S 1.000E-10-1.000J
    R/S 1
    R/S END

Method

The 5th root is numerically factored, although there are other methods. [3]

Program listing

  • Requires
    • X-Functions module on the HP-41cv
  • Available as
  • Size
    • The program takes 527 bytes including the END, and fits on 5 magnetic tracks.
  • Registers used
    • R00, constant
    • R01, 1st order coefficient
    • :
    • R05, 4th order coefficient
    • R06, highest order
    • R07, loop variable for coefficient input
    • R08, subroutine label (highest coefficient + 10)
    • R09, display loop variable (dd)
  • Flags used
    • flag 00, when set R00/R01 holds complex roots (R00+-R01.j), otherwise R00/R01 holds real roots
    • flag 02, when set, R02/R03 holds complex roots (R02+-R03.j), otherwise R02/R03 holds real roots
01	LBL "POL5"
02	SF 00		; no complex root in R00/R01
03	FIX 0
04	CF 29
05	"ORDER ?"
06	PROMPT
07	STO 06		; highest order coefficient
08	STO 07		; loop variable for coefficient input
09	 E1
10	+
11	STO 08

12	LBL 21		; coefficients input
13	"X^"
14	ARCL 07
15	>"= ?"
16	PROMPT
17	FS?C 00		; complex root in R00/R01?
18	GTO 20
19	RCL IND 06	; convert to depressed form
20	/
21	LBL 20
22	STO IND 07
23	DSE 07		; decrement loop variable
24	GTO 21
25	"CONST ?"
26	PROMPT
27	RCL IND 06
28	/
29	STO 00

30	SF 29		; "digit grouping"
31	FIX 3
32	 E		; x=1
33	STO IND 06	; highest order coefficient is 1 (because depressed form)
34	-1
35	STO 09		; init dd to -1
36	GTO IND 08	; jump to label matching highest order of polynomial

; SOLVE 2nd ORDER

37	LBL 12
38	XEQ 08		; find two roots
39	XEQ 09		; display two roots
40	GTO 11

; SOLVE 3rd ORDER

41	LBL 13
42	XEQ 10		; Find 3 roots
43	XEQ 09		; display first two roots
44	XEQ 16		; display next root
45	GTO 11

; SOLVE 4th ORDER

46	LBL 14
47	XEQ 18		; find four roots
48	XEQ 09		; display first two roots
49	XEQ 09		; display next two roots
50	GTO 11

; SOLVE 5th ORDER

51	LBL 15
52	XEQ 19		; find five roots
53	XEQ 09		; display first two roots
54	XEQ 09		; display next two roots
55	XEQ 16		; display next root

; FINISH UP, and get ready to do it again

56	LBL 11
57	"END"
58	PROMPT
59	XEQ "POL5"

; DISPLAY TWO ROOTS

60	LBL 09
61	CLA
62	ISG 09		; dd++
63	""              ; NOP
64	FC?C IND 09	; if they are not a complex roots,
65	GTO 17		;   then display two real roots

66	RCL IND 09	; display complex root in rectangular form
67	ISG 09
68	""
69	RCL IND 09
70	ABS
71	ARCL Y
72	>"+"
73	ARCL X
74	>"J"
75	AVIEW
76	STOP
77	CLA		; display conjugate of that root in rectangular form
78	ARCL Y
79	CHS
80	ARCL X
81	>"J"
82	PROMPT
83	RTN

; DISPLAY TWO REAL ROOTS

84	LBL 17
85	DSE 09		; dd
86	""              ; NOP
87	XEQ 16		; display next root, then fall through to show the next one

; DISPLAY ONE REAL ROOT

88	LBL 16
89	CLA
90	ISG 09		; dd points to next root
91	""              ; NOP
92	ARCL IND 09	; display register pointed to by dd
93	PROMPT
94	RTN

; FIND FIVE ROOTS
;
; on entry:
;   depressed form coefficients a0, a1, a2, a3, a4 in R00 .. R04
;
; on exit:
;   if FS?00 then root1 = (R00 + R01.j) and root2 = (R01 - R01.j)
;            else root1 = (R00 +   0.j) and root2 = (R01 -   0.j)
;   if FS?02 then root3 = (R02 + R03.j) and root4 = (R02 - R03.j)
;            else root3 = (R02 +   0.j) and root4 = (R03 -   0.j)
;   root5 = (R04 + 0.j)

95	LBL 19
96	RCL 00
97	STO N		; N = a0
98	 E
99	%
100	STO M		; M = a0 / 100
101	CLST

102	LBL 00		; on entry Z=1
103	RCL Z		; X=Z
104	STO O		; O = starts with 1
105	4   		; fn = a0 + n*(a1 + n*(a2 + n*(a3 + n*(a4 + n))))
106	RCL N
107	LBL 01
108	RCL IND Y
109	+
110	RCL N
111	*
112	DSE Y
113	GTO 01
114	RCL 00
115	+
116	ST* M		; M = M * fn
117	ST- O		; O = O - fn
118	RCL M
119	RCL O
120	X#0?		; if (O <> 0 )
121	/		;   then M = M / O
122	STO M		;   else M = 0

123	X<> N           ; N = N + M
124	ST+ N
125	RCL N
126	X#Y?		; if ( N <> M )
127	GTO 00		;   then next iteration

128	RCL O		; Eliminate 5th root
129	R^  		;   a3' = a4 + root5->re
130	*		;   a2' = a3 + root5->re * a3
131	+		;   a1' = a2 + root5->re * a2
132	X<> 04          ;   a0' = a1 + root5->re * a1
133	3
134	X<>Y            ; root5 = N + O * fn
135	 E
136	LBL 02
137	RCL 04
138	*
139	+
140	ENTER^
141	X<> IND Z
142	X<>Y
143	DSE Z
144	GTO 02
145	RCL 04
146	*
147	+
148	STO 00		; rolls over to "FIND FOUR ROOTS"

; FIND FOUR ROOTS
;
; on entry:
;   depressed form coefficients a0, a1, a2, a3 in R00 .. R03
;
; on exit:
;   if FS?00 then root1 = (R00 + R01.j) and root2 = (R01 - R01.j)
;            else root1 = (R00 +   0.j) and root2 = (R01 -   0.j)
;   if FS?02 then root3 = (R02 + R03.j) and root4 = (R02 - R03.j)
;            else root3 = (R02 +   0.j) and root4 = (R03 -   0.j)

149	LBL 18
150	RCL 01
151	STO N		; N = a1
152	X^2
153	RCL 00
154	STO M		; M = a0
155	RCL 02
156	CHS
157	STO O		; O = -a2
158	STO 02
159	4
160	*
161	RCL 03
162	ST* 01
163	X^2
164	+
165	*
166	+
167	CHS
168	X<> 00          ; Now:
169	4		;   a2' = -a2                  (R02)
170	*		;   a1' = a1*a3-4*a0	       (R01)
171	ST- 01		;   a0' = a0*(4*a2-a3^2)-a1^2  (R00)

172	XEQ 07		; Find 3 roots

173	RCL 02		; Determine largest root
174	FS? 00		;   if complex root in R00/R01, then
175	GTO 06          ;     X = root3->re
176	RCL 01		;   else
177	X<Y?            ;     X = maximum of root1->re, root2->re, root3->re
178	X<>Y
179	RCL 00
180	X<Y?
181	X<>Y            ; X = largestRoot

182	LBL 06 		; Eliminate the largest root
183	STO 00		;   A = a3/2
184	STO 02     	;   B = largestRoot/2
185	RCL O		;   C = sqrt(A*A-a2+largestRoot);
186	+   		;   D = sqrt(B*B-a0);
187	RCL 03		;   S = sign(A*B-a1/2))    (sign returns 1 if >=0, -1 if  02
225	STO 00
226	RCL 01
227	X<> 03
228	STO 01
229	CF 02
230	FS? 00
231	SF 02		; Find the next two roots
232	CLA		;   a1' = A - S.C;
233	GTO 03 		;   a0' = B - D;

; FIND THREE ROOTS
;
; on entry:
;   depressed form coefficients a0, a1, a2 in R00 .. R02
;
; on exit:
;   if FS?00 then root1 = (R00 + R01.j) and root2 = (R01 - R01.j)
;            else root1 = (R00 +   0.j) and root2 = (R01 -   0.j)
;   root3 = (R02 + 0.j)

234	LBL 10		; Find the 3rd root
235	LBL 07
236	XEQ 04
237	RCL 02
238	3
239	/
240	-

241	STO 00		; Eliminate the 3rd root
242	ST+ 02
243	RCL 02
244	*
245	ST+ 01
246	RCL 02
247	X<> 01          ; Find the next two roots (rolls over to LBL 08)
248	X<> 00          ;     a1' = a2 + root3->re                           (R01)
249	STO 02		;     a0' = a1 + root3->re * ( a2 + root3->re)       (R00)

; FIND TWO ROOTS
;
; on entry:
;   depressed for coefficients a0, a1 in R00 .. R01
;
; on exit:
;   if FS?00 then root1 = (R00 + R01.j) and root2 = (R01 - R01.j)
;            else root1 = (R00 +   0.j) and root2 = (R01 -   0.j)

250	LBL 08
251	LBL 03
252	CF 00
253	RCL 00
254	4
255	*
256	RCL 01
257	STO 00
258	X^2
259	-		; discr = (a1^2)/4 - a0;
260	X>0?            ; discr > 0?
261	SF 00
262	ABS
263	SQRT
264	2
265	CHS
266	ST/ 00
267	/
268	CHS
269	STO 01
270	FS? 00		; complex root in R00/R01?
271	RTN 		; X = sqrt(abs(a0-a1^2/4))

272	RCL 00
273	STO 01
274	RDN
275	ST+ 00
276	ST- 01
277	RTN 		; X = sqrt(abs(a0-a1^2/4))

; FIND 3RD ROOT HELPER, returns the real part + a2/3 in X
;
; on entry:
;   depressed form coefficients a0, a1, a2 in R00 .. R02
;
; on exit:
;   Note the ' sign, because to get root3, 2/3 needs to be subtracted
;   X register holds (root3 + a2/3)

278	LBL 04
279	RCL 01	        ; Method:
280	RCL 02	        ;   p = a1/3 - a2^2/9
281	X^2 		;   q = a1*a2/6 - a2^3/27 - a0/2;
282	3		;   tmp = p^3 + q^2
283	/
284	-
285	3
286	/
287	RCL 01
288	RCL 02
289	*
290	6
291	/
292	RCL 02
293	3
294	/
295	LASTX
296	Y^X
297	-
298	RCL 00
299	2
300	/
301	-
302	X=0?		; if (q==0) then
303	RTN		;     X = (0 + 0.j)

304	STO Z  	        ; else
305	X^2
306	X<>Y
307	3
308	Y^X
309	+		;     if (tmp < 0)
310	Xre' = 2*(q^2+sqrt(-tmp)^2)^1/3 *
311	GTO 05		;             cos(phi(q,sqrt(-tmp))/3)

312	SQRT		;     else
313	ST- Z		;         X = sign(q-sqrt(tmp))*abs(q-sqrt(tmp))^1/3 +
314	+   		;             sign(q+sqrt(tmp))*abs(q+sqrt(tmp))^1/3
315	SIGN
316	LASTX
317	ABS
318	3
319	1/X
320	Y^X
321	*
322	X<>Y
323	SIGN
324	LASTX
325	ABS
326	3
327	1/X
328	Y^X
329	*
330	+
331	RTN

332	LBL 05
333	CHS
334	SQRT
335	X<>Y
336	R-P
337	3
338	1/X
339	Y^X
340	X<>Y
341	3
342	/
343	COS
344	*
345	ST+ X
346	RTN

347	END

References

[1] Title Unknown
Michiel Niemeijer, approx 1985
HP User News, TH Boekhandel Prins, Binnenwatersloot 30, Delft, the Netherlands
[2] Mathematics of the Discrete Fourier Fransform (DFT), Factoring a Polynomial
Julius O. Smith III, Stanford University, April 2007 (2nd edition)
W3K Publishing, ISDN 978-0974560748
[3] Quintic Equation (5th degree)
Wolfram
Coert Vonk

Coert Vonk

Independent Firmware Engineer at Los Altos, CA
Welcome to the things that I couldn’t find.This blog shares some of the notes that I took while deep diving into various fields.Many such endeavors were triggered by curious inquiries from students. Even though the notes often cover a broader area, the key goal is to help the them adopt, flourish and inspire them to invent new technology.
Coert Vonk

Latest posts by Coert Vonk (see all)

Leave a Reply

You can use these HTML tags

<a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

  

  

  

Protected with IP Blacklist CloudIP Blacklist Cloud