Numerical computation of polynomial zeros by means of Aberth's method

An algorithm for computing polynomial zeros, based on Aberth's method, is presented. The starting approximations are chosen by means of a suitable application of Rouche's theorem. More precisely, an integer q greater than or equal to 1 and a set of annuli A(i), i = 1,...,q, in the complex plane, are determined together with the number ici of zeros of the polynomial contained in each annulus A(i). As starting approximations we choose Ici complex numbers lying on a suitable circle contained in the annulus A(i), for i = 1,...,q. The computation of Newton's correction is performed in such a way that overflow situations are removed. A suitable stop condition, based on a rigorous backward rounding error analysis, guarantees that the computed approximations are the exact zeros of a ''nearby'' polynomial. This implies the backward stability of our algorithm. We provide a Fortran 77 implementation of the algorithm which is robust against overflow and allows us to deal with polynomials of any degree, not necessarily monic, whose zeros and coefficients are representable as floating point numbers. In all the tests performed with more than 1000 polynomials having degrees from 10 up to 25,600 and randomly generated coefficients, the Fortran 77 implementation of our algorithm computed approximations to all the zeros within the relative precision allowed by the classical conditioning theorems with 11.1 average iterations. In the worst case the number of iterations needed has beer, at most 17. Comparisons with available public domain software and with the algorithm PA16AD of Harwell are performed and show the effectiveness of our approach. A multiprecision implementation in MATHEMATICA(TM) is presented together with the results of the numerical tests performed.