Design, analysis, and implementation of a multiprecision polynomial rootfinder

We present the design, analysis, and implementation of an algorithm for the computation of any number of digits of the roots of a polynomial with complex coefficients. The real and the imaginary parts of the coefficients may be integer, rational, or floating point numbers represented with an arbitrary number of digits. The algorithm has been designed to deal also with numerically hard polynomials like those arising from the symbolic preprocessing of systems of polynomial equations, where the degree and the size of the coefficients are typically huge. The algorithm is based on an adaptive strategy which automatically exploits any specific feature of the input polynomial, like its sparsity or the conditioning of its roots, in order to speed up the computation. We introduce different concepts and tools suitably designed to arrive at an adaptive implementation, such as the concepts of epsilon-root neighborhood, epsilon-inclusion discs and some inclusion and conditioning theorems for their determination. The main engine for shrinking the inclusion discs is the simultaneous iteration method of Ehrlich-Aberth, complemented with a suitable technique for cluster analysis that is used for getting rid of the slow convergence in case of clustered or multiple roots. The algorithm, implemented in C, relies on the GNU multiprecision package GMP and allows many options. Counting, isolating and approximating all roots in a given set \mathcalS are the main goals that the algorithm provides. Automatic determination of multiplicities and the detection of real or imaginary roots can be selected as well. Polynomials having coefficients with a bounded precision may be processed too. Comparisons with the polynomial rootfinders of the packages Mathematica(TM), Maple(TM) and Pari, performed on a wide class of test polynomials show that our algorithm is generally much faster: in most cases the speed-up factor is greater than 10 and, for certain polynomials, it is greater than 1000. The MPSolve package can be downloaded from the numeralgo library of netlib.