(71B) Aberth method

[robve]

ABERTH POLYNOMIAL ROOT FINDER
Aberth method: Durand-Kerner with Newton's method
Estimate the complex roots of a polynomial with real coefficients
Note: Durand-Kerner is not guaranteed to converge

Optimized implementation:
- Finds complex roots z of polynomial p such that |p(z)|<T given tolerance T
- Prevents overflow/nan/inf by bounding the root shift adjustment magnitude during polishing
- Tested with a large sample of polynomials

For details see: https://www.hpmuseum.org/forum/thread-18051.html
See also: https://en.wikipedia.org/wiki/Durand–Kerner_method and https://en.wikipedia.org/wiki/Aberth_method

For HP-71B + Math Pac

Degree: max polynomial degree limited by time and memory only
Tolerance: T=1E-6 (adjustable)
Levels: up to L=20 iterations (adjustable)

Example:

(Solve x^2+2x+3=0)
[RUN]
N=2 A(0)=1 A(1)=2 A(2)=3
T=1E-6 L=20
(-1,-1.414214)
(-1,1.414214)

(Solve x^3-2x^2-x+2=0)
[RUN]
N=3 A(0)=1 A(1)=-2 A(2)=-1 A(3)=2
T=1E-6 L=20
(1,0)
(2,0)
(-1,0)

(Solve x^3+3x^2+x+3=0)
[RUN]
N=3 A(0)=1 A(1)=3 A(2)=1 A(3)=3
T=1E-6 L=20
(0,1)
(0,-1)
(-3,0)

(Solve x^4+2999x^3-10003e3x^2-2399e7x+24e9=0)
[RUN]
N=4 A(0)=1 A(1)=2999 A(2)=-10003E3 A(3)=-2399E7 A(4)=24E9
T=1E-6 L=20
(1,0)
(-4000,0)
(-2000,0)
(3000,0)

(Solve 5x^6-45x^5+225x^4-425x^3+170x^2+370x-500=0)
[RUN]
N=6 A(0)=5 A(1)=-45 A(2)=225 A(3)=-425 A(4)=170 A(5)=370 A(6)=-500
T=1E-6 L=20
(1,1)
(3,4)
(-1,0)
(3,-4)
(1,-1)
(2,0)

Code:

' ABERTH POLYNOMIAL ROOT FINDER
' Estimate the complex roots of a polynomial with real coefficients

' For HP-71B + Math Pac by Robert van Engelen

' Optimized implementation:
' - Finds complex roots z of polynomial p such that |p(z)|<T given tolerance T
' - Prevents overflow/nan/inf by bounding the root stride modulus during polishing
' - Tested with a large sample of polynomials
'
' For details see: https://www.hpmuseum.org/forum/thread-18051.html
' See also: https://en.wikipedia.org/wiki/Durand–Kerner_method
'           https://en.wikipedia.org/wiki/Aberth_method
'
' VARIABLES
'  N         degree
'  T         tolerance (1E-6)
'  L         max iterations (20), updated to the iterations executed
'  A(0..N)   N+1 real coefficients indexed from highest to lowest
'  Z(0..N-1) N complex roots
'  D         denominator
'  E         estimated error
'  H         root stride modulus
'  P         p(z) and root stride p(z)/d
'  Q         p'(z)
'  R         root
'  U         max root stride modulus
'  I,J,K     loop iterators
'  M         number of coefficients - 1 (M<=N-1)
'  F         scratch
'
' DELAY 1 or smaller recommended

100 OPTION BASE 0
110 DESTROY ALL
120 INTEGER I,J,K,L,M,N @ REAL E,H,T,U
130 INPUT "N=";N @ IF N<2 THEN 130
140 REAL A(N) @ M=N-1 @ COMPLEX D,F,P,Q,R,Z(M)
150 MAT INPUT A
160 INPUT "T=","1E-6";T @ INPUT "L=","20";L
' remove leading zero coefficients (trailing A)
170 IF M>0 AND A(M+1)=0 THEN Z(M)=0 @ M=M-1 @ GOTO 170
180 IF M<0 OR A(0)=0 THEN DISP "ERROR" @ END
' generate initial roots
190 R=1 @ FOR I=0 TO M @ R=(.4,.9)*R @ Z(I)=R @ NEXT I
' max root stride modulus bounded to prevent overflow
200 U=MAXREAL^(1/(M+1)) @ IF U<1 THEN U=1
' iterative root polishing until convergence or max iterations reached
210 FOR K=1 TO L @ E=0
220 FOR I=0 TO M
' compute p(z) and p'(z) using Horner's rule
230 P=A(0) @ Q=0 @ R=Z(I) @ D=0
240 FOR J=0 TO M @ Q=Q*R+P @ P=P*R+A(J+1) @ NEXT J
' optimized: if |p(z)|<tolerance then do not update root z
250 IF ABS(P)<T THEN 310
' compute denominator
260 FOR J=0 TO M @ IF I=J THEN 280
270 F=R-Z(J) @ IF ABS(F)>=T THEN D=D+P/F
280 NEXT J
' compute p(z)/d if d is not too small to prevent overflow
290 D=Q-D @ IF ABS(D)>T THEN P=P/D
' polish root if its root stride modulus is within max bounds to prevent overflow
300 H=ABS(P) @ E=E+H @ IF H<U THEN Z(I)=R-P
310 NEXT I
320 DISP "K=";K;"E=";E
330 IF E<=T THEN J=L @ L=K @ K=J
340 NEXT K
' round results (optional)
350 U=1/T
360 FOR I=0 TO M @ Z(I)=(T*INT(U*REPT(Z(I))+.5),T*INT(U*IMPT(Z(I))+.5)) @ NEXT I
' estimate actual error
370 E=E+(E=0)*T
' display results
380 MAT DISP Z