Polynomial Roots Reloaded (yes, seriously...)

[Ángel Martin]

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  END

As 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.