Solving secular and polynomial equations: A multiprecision algorithm

We present an algorithm for the solution of polynomial equations and secular equations of the form S(x)=0S(x)=0 for S(x)=∑i=1naix−bi−1=0, which provides guaranteed approximation of the roots with any desired number of digits. It relies on the combination of two different strategies for dealing with the precision of the floating point computation: the strategy used in the package MPSolve of D. Bini and G. Fiorentino [D.A. Bini, G. Fiorentino, Design, analysis and implementation of a multi-precision polynomial rootfinder, Numer. Algorithms 23 (2000) 127–173] and the strategy used in the package Eigensolve of S. Fortune [S. Fortune, An iterated eigenvalue algorithm for approximating the roots of univariate polynomials, J. Symbolic Comput. 33 (5) (2002) 627–646]. The algorithm is based on the Ehrlich–Aberth (EA) iteration, and on several results introduced in the paper. In particular, we extend the concept and the properties of root-neighborhoods from polynomials to secular functions, provide perturbation results of the roots, obtain an effective stop condition for the EA iteration and guaranteed a posteriori error bounds. We provide an implementation, released in the package MPSolve 3.0, based on the GMP library. From the many numerical experiments it turns out that our code is generally much faster than MPSolve 2.0 and of the package Eigensolve. For certain polynomials, like the Mandelbrot or the partition polynomials the acceleration is dramatic. The algorithm exploits the parallel architecture of the computing platform.