Numerical computation of polynomial zeros by means of Aberth's method (Q676928)

The method of \textit{O. Aberth} [Math. Comput. 27, 339-344 (1973; Zbl 0282.65037)] is a modification of the Laguerre method that avoids the computation of the second derivative and allows the simultaneous determination of all roots of a polynomial. The author introduces several improvements to the method. First, instead of looking for good unique starting values, he uses known inequalities for the coefficients to determine circular rings in which a number of roots is guaranteed to be found and, therefore, starts his method at multiple places. Second, the essential quantity to be computed in each iteration is the Newton quotient \(p(x)/p'(x)\) that can be cheaply obtained by Horner's method but that method is numerically unstable for \(|x|>1\). Let the degree of \(p\) be \(n\) and \(q(n)=x^np(1/x)\) the reverse polynomial. Then \(p(x)/p'(x)=[ny-y^2q(y)/q'(y)]^{-1}\) where \(y=1/x\). Hence, the Horner scheme can be used for numbers \(|x|\leq 1\) and is always numerically stable. The complete algorithm is outlined and it is claimed that in numerical implementation it improves even on Jenkins-Traub, in particular for very large \(n\) \((>850)\). A stopping criterion is developed that guarantees the accuracy of the compound roots. \{Now what the criterion really says is that the computed value of the root is the exact root of a slightly perturbed polynomial. But the famous Wilkinson polynomial shows that a minute change of coefficients can have catastrophic consequences for the roots so that the problem of determination of roots of polynomials of high degree is intrinsically extremely badly conditioned.\}.