# Polynomial factorization on HP-41

Scientific calculator software for polynomial factorization on HP-41. A polynomial of degree n has always n roots. Up to 5th degree.

## History

Around 1985, Michiel Niemeyer published a similar program for the HP-41 in the Dutch publication HP User News . I changed it to make the UX more intuitive. This article lists my version 1.1 as of 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 

$$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. The polynomial can be factorized as in

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

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}$$

The generic form for a 2nd order polynomial 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 third order polynomials [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

$$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

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{\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}$$ gives one of the roots

\begin{aligned} x_1&=w_1-\frac{p}{3w_1}-\frac{b}{3a}\\ w_1&= \sqrt{\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{aligned}

Polynomial division removes the root, and leaves a quadratic polynomial

### 5th order

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

## Examples

Assumes FIX 3 and SIZE 15.

### Example 1

Calculate the roots for

$$x^5-x^4-101x^3+101x^2+100x-100 = 0$$

The solution is shown below.

Keystrokes
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

### Example 2

Calculate the roots for

$$x^4+x^3+x^2+11x+10 = 0$$

The solution is shown below.

Keystrokes
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

### Example 3

Calculate the roots for

$$x^3-x^2+x-1 = 0$$

The solution is shown below.

Keystrokes
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

## Source code

• Requires X-Functions module on the HP-41cv
• Available as source code, raw for the V41 emulator and bar code for the HP Wand (HP82153A)
• .
• 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

### Listing

Available through the repository

 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

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

This site uses Akismet to reduce spam. Learn how your comment data is processed.