Finding one root The most widely used method for computing a root of any differentiable function f is
Newton's method, in which an initial guess x_0 is iteratively refined. At each iteration the
tangent line to f at x_n is used as a
linear approximation to f, and its root is used as the succeeding guess x_{n+1}: :x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}, In general, the value of x_n will converge to a root of f. In particular, the method can be applied to compute a root of polynomial functions. In this case, the computations in Newton's method can be accelerated using
Horner's method or
evaluation with preprocessing for computing the polynomial and its derivative in each iteration. Though the rate of convergence of Newton's method is generally
quadratic, it might converge much slowly or even not converge at all. In particular, if the polynomial has no real root, and x_0 is chosen to be a real number, then Newton's method cannot converge. However, if the polynomial has a real root, which is larger than the larger real root of its derivative, then Newton's method converges quadratically to this largest root if x_0 is larger than this larger root (there are easy ways for computing an upper bound of the roots, see
Properties of polynomial roots). This is the starting point of
Horner's method for computing the roots. Closely related to Newton's method are
Halley's method and
Laguerre's method. Both use the polynomial and its two first derivations for an iterative process that has a
cubic convergence. 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. This gives a slight advantage to these methods (less clear for Laguerre's method, as a square root has to be computed at each step). When applying these methods to polynomials with real coefficients and real starting points, Newton's and Halley's method stay inside the real number line. One has to choose complex starting points to find complex roots. In contrast, the Laguerre method with a square root in its evaluation will leave the real axis of its own accord.
Finding all complex roots Methods using complex-number arithmetic Both the
Aberth method and the similar yet simpler
Durand–Kerner method simultaneously find all of the roots using only simple
complex number arithmetic. The Aberth method is presently the most efficient method. 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. A
free implementation of Aberth's method is available under the name of
MPSolve. 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. Another method with this style is the
Dandelin–Gräffe method (sometimes also ascribed to
Lobachevsky), which uses
polynomial transformations to repeatedly and implicitly square the roots. This greatly magnifies variances in the roots. Applying
Viète's formulas, one obtains easy approximations for the modulus of the roots, and with some more effort, for the roots themselves.
Methods using linear algebra Arguably, the most reliable method to find all roots of a polynomial is to find the eigenvalues of the
companion matrix of monic polynomial, which coincides with the roots of the polynomial. There are plenty of algorithms for computing the eigenvalue of matrices. The standard method for finding all roots of a polynomial in MATLAB uses the Francis
QR algorithm to compute the
eigenvalues of the corresponding companion matrix of the polynomial. By the use of the sparse structure of the companion matrix, some special QR iteration methods give all complex roots in O(n^3) arithmetics and O(n) storage. In principle, one can use any
eigenvalue algorithm to find the roots of the polynomial. However, for efficiency reasons one prefers methods that employ the structure of the matrix, that is, can be implemented in matrix-free form. Among these methods are the
power method, whose application to the transpose of the companion matrix is the classical
Bernoulli's method to find the root of greatest modulus. 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. Additionally, it has fast convergence with order 1+\varphi\approx 2.6 (where \varphi is the
golden ratio) even in the presence of clustered roots. This fast convergence comes with a cost of three polynomial evaluations per step, resulting in a residual of , that is a slower convergence than with three steps of Newton's method.
Limitations of iterative methods for finding all roots The oldest method of finding all roots is to start by finding a single root. When a root has been found, it can be removed from the polynomial by dividing out the binomial . The resulting polynomial contains the remaining roots, which can be found by iterating on this process. This idea, despite being common in theoretical deriviations, does not work well in numerical computations because of the phenomenon of
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). Also, even with a good approximation, when one evaluates a polynomial at an approximate root, one may get a result that is far too close to zero. For example, if a polynomial of degree 20 (the degree of Wilkinson's polynomial) has a root close to 10, the derivative of the polynomial at the root may be of the order of 10^{20}; this implies that an error of 10^{-10} on the value of the root may produce a value of the polynomial at the approximate root that is of the order of 10^{10}.
Finding all real roots Finding the real roots of a polynomial with real coefficients is a problem that has received much attention since the beginning of 19th century, and is still an active domain of research. Methods for finding all complex roots can provide the real roots. However, because of the numerical instability of polynomials, it may need
arbitrary-precision arithmetic to decide whether a root with a small imaginary part is real or not. Moreover, as the number of the real roots is, on the average, proportional to the logarithm of the degree, it is a waste of computer resources to compute the non-real roots when one is interested in real roots. The standard way of computing real roots is to compute first disjoint intervals, called
isolating intervals, such that each one contains exactly one real root, and together they contain all the roots. This computation is called
real-root isolation. Having an isolating interval, one may use fast numerical methods, such as
Newton's method for improving the precision of the result. 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 its extensions—
Budan's and
Vincent's theorems. These methods divide into two main classes, one using
continued fractions and the other using bisection. Both method have been dramatically improved since the beginning of 21st century. With these improvements they reach a
computational complexity that is similar to that of the best algorithms for computing all the roots (even when all roots are real). These algorithms have been implemented and are available in
Mathematica (continued fraction method) and
Maple (bisection method), as well as in other main
computer algebra systems (
SageMath,
PARI/GP) . Both implementations can routinely find the real roots of polynomials of degree higher than 1,000.
Finding roots in a restricted domain Several fast tests exist that tell if a segment of the real line or a region of the complex plane contains no roots. By bounding the modulus of the roots and recursively subdividing the initial region indicated by these bounds, one can isolate small regions that may contain roots and then apply other methods to locate them exactly. All these methods involve finding the coefficients of shifted and scaled versions of the polynomial. For large degrees,
FFT-based accelerated methods become viable. The
Lehmer–Schur algorithm uses the
Schur–Cohn test for circles; a variant,
Wilf's global bisection algorithm uses a winding number computation for rectangular regions in the complex plane. The
splitting circle method uses FFT-based polynomial transformations to find large-degree factors corresponding to clusters of roots. The precision of the factorization is maximized using a Newton-type iteration. This method is useful for finding the roots of polynomials of high degree to arbitrary precision; it has almost optimal complexity in this setting.
Finding complex roots in pairs If the given polynomial only has real coefficients, one may wish to avoid computations with complex numbers. To that effect, one has to find quadratic factors for pairs of conjugate complex roots. The application of the
multidimensional Newton's method to this task results in
Bairstow's method. The real variant of
Jenkins–Traub algorithm is an improvement of this method.
Polynomials with rational coefficients For polynomials whose coefficients are exactly given as
integers or
rational numbers, there is an efficient method to factorize them into factors that have only simple roots and whose coefficients are also given in precise terms. This method, called
square-free factorization, is based on the multiple roots of a polynomial being the roots of the
greatest common divisor of the polynomial and its derivative. The square-free factorization of a polynomial
p is a factorization p=p_1p_2^2\cdots p_k^k where each p_i is either 1 or a polynomial without multiple roots, and two different p_i do not have any common root. An efficient method to compute this factorization is
Yun's algorithm. == See also==