In 1996, V. Y. Pan wrote a paper that described "algorithms for approximating all the complex zeros of an $n^{th}$ degree polynomial $p(x)$". His algorithms assume that all zeros are within the unit disc. He notes:
[O]ur algorithms (which promise to be practically effective) approximate all the zeros of $p(x)$ within the absolute error bound $2^{-b}$, by using order of $n$ arithmetic operations and order of $(b + n)n^2$ Boolean (bitwise) operations (in both cases up to within polylogarithmic factors).
He also compares this algorithm to a couple others:
In particular, for practical approximation of complex polynomial zeros, the most promising alternative to the approach of this paper probably comes from the Durand-Kerner algorithm and its various modifications (such as Aberth's...), which rely on Newton's iteration for Viète's system of polynomial equations for the zeros of $p(x)$. The absence of global convergence proofs and of any reasonably good computational complexity estimates for these iterative algorithms is partly compensated by their very good record in numerical experiments. On the other hand, these algorithms require us to use either the order of $n^2$ arithmetic operations per iteration (which is roughly $n$ times as many as we use in our entire algorithms) or a much higher precision of computing (to support application of fast multipoint polynomial evaluation, which is a basic step of these algorithms).
However, more than 20 years later, I'm having trouble finding an implementation of this algorithm. Mathematica's NRoots allows one to specify one of three possible Methods: "Aberth", "CompanionMatrix", or "JenkinsTraub". In Python, mpmath's polyroots uses the Durand-Kerner method.
I have a couple questions:
- Why aren't these algorithms more commonly used? Why aren't they the standard in Mathematica or mpmath? Is there some practical reason why these algorithms perform poorly? (That is, do they perform worse in practice than in theory?) Is there some other theoretical issue I'm missing for why these algorithms are not the standard?
- Does anyone know of any implementation of these algorithms? Are there other algorithms that are provably as good? If I want to get the $\widetilde{\mathcal O}(n^3 + n^2b)$ running time, do I have to implement these algorithms myself?
Answers to any parts of either of these questions would be appreciated.