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 [1]. 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 [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. 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}
$$
Quadratic polynomial (2nd order)
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[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}\) gives one of the roots
$$
\begin{aligned}
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{aligned}
$$
Polynomial division removes the root, and leaves a quadratic polynomial
5th order
The 5th root is numerically factored, although there are other methods. [3]
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.
| 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.
| 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.
| 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
| [1] | Title Unknown Michiel Niemeijer, approx 1985 HP User News, TH Boekhandel Prins, Binnenwatersloot 30, Delft, the Netherlands |
|
| [2] | 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 |
|
| [3] | Quintic Equation (5th degree) Wolfram |