# Polynomial root finders with consistent root ordering

When it comes to solving simple polynomials with complex coefficients and complex roots, we can use explicit formulae. For example, the roots for polynomial $z^2-c$ are the list $(\sqrt c,-\sqrt c)$. A small perturbation in the coefficients will result in a small perturbation of the roots, but the list will not be permuted.

For polynomials of higher degree we will use some kind of numerical procedure that gives us a list of roots. But running the numerical procedure on the same polynomial a second time may give us the same roots in a different order (because it may have used randomness in selecting the starting point of an iteration). And even if that's not the case, there is no guarantee that perturbing the polynomial slightly will yield the perturbed roots in the same order as the roots of the original.

For example, a numerical solver might give the roots of $z^2-c$ as $(\operatorname N(\sqrt c),-\operatorname N(\sqrt c))$ and the roots of $z^2-(c+\epsilon)$ as $(-\operatorname N(\sqrt {c+\epsilon},\operatorname N(\sqrt {c+\epsilon}))$, where $\operatorname N()$ is a numerical approximation.

Are there numerical root finders that preserve root order as well as formula based solvers? (There are some complexities that are discussed in the comments ... it impossible for roots to vary continuously with coefficients everywhere - even formula based solvers will have branch cuts)

The polynomials I'm interested in are of small degree - less than 10. And I would be happiest if the answer also pointed to a reliable open source implementation.

• You are assuming that the roots are real and will stay real under small perturbation. You should read about Wilkinson's polynomial, though it is of degree $20$.. All its roots are real and well separated, but perturbing one coefficient by a relative error of $2^{-29}$ sends the roots off into the complex plane, four with imaginary parts greater than $2$. – Ross Millikan Apr 14 '18 at 3:13
• The perturbation-invariant part of this is impossible for topological reasons. Consider the polynomial family $z^2-e^{it}$. If the algorithm produces the result $(1,-1)$ when $t=0$ and is invariant under small perturbations of $t$, it must produce the result $(e^{it/2},-e^{it/2})$. But that means that when you get to $t=2\pi$ you get the same roots in the opposite order, for the same polynomial. – Micah Apr 14 '18 at 3:16
• @Ross, I'm not assuming the roots are real. If some part of the question implies that, let me know and I will try to fix it. – brainjam Apr 14 '18 at 3:37
• If the roots aren't real the concept of order becomes more difficult because the complex numbers do not have an obvious order. In the Wilkinson example it is not easy to say which root in the original polynomial corresponds to which root in the perturbed polynomial, so it is not clear what order you want the roots of the perturbed polynomial produced in. – Ross Millikan Apr 14 '18 at 3:46
• @Micah, good point. I'll admit I'm struggling to formulate the right question, and maybe the perturbation idea is the wrong way to go. I've got a vague idea that if the polynomial is an $n$-fold covering with $n$ sheets that the roots should be in 'sheet order' if that makes any sense. So it's ok for there to be the kind of discontinuity that your example would imply - at cuts between the sheets. – brainjam Apr 14 '18 at 3:47