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)=0$ for $S(x)=sum_{i=1}^nrac{a_i}{x-b_i}-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 [Numer.~Algo.~23, 2000] and the strategy used in the package Eigensolve of S.~Fortune [J. Symb. Comput.~33, 2002]. 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 ounds. 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.