Finding polynomial roots is a long-standing problem that has been the object of much research throughout history.
Finding the root of a linear polynomial (degree one) is easy and needs only one division: the general equation
For polynomials of degree five or higher Abel–Ruffini theorem asserts that there is, in general, no radical expression of the roots.
By the fundamental theorem of algebra, a polynomial of degree n has exactly n real or complex roots counting multiplicities.
However, except for low degrees, this does not work well because of the numerical instability: Wilkinson's polynomial shows that a very small modification of one coefficient may change dramatically not only the value of the roots, but also their nature (real or complex).
For avoiding these problems, methods have been elaborated, which compute all roots simultaneously, to any desired accuracy.
This is a reference implementation, which can find routinely the roots of polynomials of degree larger than 1,000, with more than 1,000 significant decimal digits.
The oldest method for computing the number of real roots, and the number of roots in an interval results from Sturm's theorem, but the methods based on Descartes' rule of signs and its extensions—Budan's and Vincent's theorems—are generally more efficient.
The main computer algebra systems (Maple, Mathematica, SageMath, PARI/GP) have each a variant of this method as the default algorithm for the real roots of a polynomial.
However, for efficiency reasons one prefers methods that employ the structure of the matrix, that is, can be implemented in matrix-free form.
The inverse power method with shifts, which finds some smallest root first, is what drives the complex (cpoly) variant of the Jenkins–Traub algorithm and gives it its numerical stability.
In fact, the problem of finding the roots of a polynomial from its coefficients can be highly ill-conditioned.
Combining two consecutive steps of these methods into a single test, one gets a rate of convergence of 9, at the cost of 6 polynomial evaluations (with Horner's rule).
On the other hand, combining three steps of Newtons method gives a rate of convergence of 8 at the cost of the same number of polynomial evaluation.
In contrast, the Laguerre method with a square root in its evaluation will leave the real axis of its own accord.
To that effect, one has to find quadratic factors for pairs of conjugate complex roots.
Accelerated algorithms for multi-point evaluation and interpolation similar to the fast Fourier transform can help speed them up for large degrees of the polynomial.
It is advisable to choose an asymmetric, but evenly distributed set of initial points.
The implementation of this method in the free software MPSolve is a reference for its efficiency and its accuracy.
Several fast tests exist that tell if a segment of the real line or a region of the complex plane contains no roots.
All these methods involve finding the coefficients of shifted and scaled versions of the polynomial.
The splitting circle method uses FFT-based polynomial transformations to find large-degree factors corresponding to clusters of roots.
This method is useful for finding the roots of polynomials of high degree to arbitrary precision; it has almost optimal complexity in this setting.
The oldest complete algorithm for real-root isolation results from Sturm's theorem.
However, it appears to be much less efficient than the methods based on Descartes' rule of signs and Vincent's theorem.
Both implementations can routinely find the real roots of polynomials of degree higher than 1,000.
Multiple roots are highly sensitive, known to be ill-conditioned and inaccurate in numerical computation in general.
The first step determines the multiplicity structure by applying square-free factorization with a numerical greatest common divisor algorithm.
The second step applies the Gauss-Newton algorithm to solve the overdetermined system for the distinct roots.
"Algorithm 835: MultRoot – A Matlab package for computing polynomial roots and multiplicities".