+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: HP-41C Software Library (/forum-11.html)
+--- Thread: Polynomial Roots Reloaded (yes, seriously...) (/thread-7763.html)
Polynomial Roots Reloaded (yes, seriously...) - Ángel Martin - 02-12-2017 05:07 PM
Yes I now, this has been rehashed a few dozen times and yet some more - but I didn't want to let this pass go unnoticed since it's also a nice application of the function set on the 41Z module - which packs so many functions that even I lose track of them.
The routine below is not a new conceptual method or a cleverly design algorithm, not at all. If anything it's a straight application of the oldest "brutish-force"approach there is, i.e. the Newton method - with a deflation method after each root is found using approximations. Why bother then?
Well, for once as a nice example of the 41Z function set, which includes MCODE functions to evaluate polynomials, their primitives and their first & second derivatives (see program steps 12 & 15 below). If anything this should make the implementation much faster than other previous attempts.
Second this is a nice tribute to JM Baillard's continuous dedication to the 41 legacy. The routine is a straight adaptation of his program dealing with real roots of polynomials, as can be seen at paragraph #1-f of this page.. If you compare both versions you'll see a direct translation from "standard" functions into 41Z-based, for the complex domain. The original routine only works to obtain real roots, but the complex-version doesn't have that limitation and works equally for real and complex cases - even if the coefficients are all real.
Here's the core routine, which assumes the polynomial coefficients are stored in Complex Data registers CR(bbb) to CR(eee) - the initial guess is the {Z,Y} stack registers, and the polynomial control word bbb.eee in the X-register (using Complex Data register indexes).
Code:
01 LBL "ZPRT bbb,000(eee)
02 STO 02
03 STO 03 reg range
04 STO 04
05 RDN
06 ZSTO (00)
07 ISG 04
08 LBL 01
09 ZRCL (00)
10 ZAVIEW
11 RCL 03 reg range
12 ZPL
13 LASTZ
14 RCL 03 reg range
15 ZPLD1 (ZF# 47)
16 Z/
17 ZST- (00)
18 ZRCL (00)
19 Z=0?
20 SIGN
21 Z/
22 ZMOD
23 E-8 tolerance
24 X<Y?
25 GTO 01
26 E-3 0,001
27 ST- 03 deflate pol
28 RCL 03 reduced degree
29 STO 05 used as index
30 CLZ
31 LBL 02
32 ZRC* (00)
33 ZST+ IND 05
34 ZRCL IND 05
35 ISG 05
36 GTO 02
37 ZRCL (00)
38 ZSTO IND 05
39 ISG 04
40 GTO 01
41 RCL 02
42 E
43 + (bbb+1),00(eee)
44 ENDAs an example let's use the same 5-degree polynomial from the 41Z manual, page.
P(x) = (1+2i)*z^4 + (-1-2i)*z^3 + (3-3i)*z^2 + z – 1
First store the five coefficients in complex registers CR03 to CR07, (can use the sub-function ZINPT). Then using (1+i) as initial guess and the control word 3,007 in the X-register, the four solutions are obtained in a few seconds on the 41-CL.
1, ENTER^, 1, 3.007, XEQ "ZPLRT"
1,698+J0,802
-0,400-J0,859
0,358+J0,130
-0,656-J0,073
Not bad for a (roughly) 44-step routine, won't you agree?
Happy root finding in the Complex realm...
ÁM.