# Finding all roots of a polynomial

Is it possible, for an arbitrary polynomial in one variable with integer coefficients, to determine the roots of the polynomial in the Complex Field to arbitrary accuracy? When I was looking into this, I found some papers on homotopy continuation that seem to solve this problem (for the Real solutions at least), is that correct? Or are there restrictions on whether homotopy continuation will work? Does the solution region need to be bounded?

-
The answer is "yes" and modern computer algebra systems have already done this for you. I confess I don't know how---but you don't make it clear whether you want to know how or you just want to know the answer. If you have a particular polynomial in mind, fire up the free maths package pari, set the precision to 1000 with \p 1000, and then use the polroots command. – Kevin Buzzard Mar 31 '10 at 20:53
The fact that it's implemented doesn't meant it is solved! CASes have all sorts of routines which 'solve' undecidable problems... because all undecidable problems have (often large) sub-classes which are semi-decidable. It turns out that, for this problem, there is a complete algorithm which is guaranteed to terminate and find all roots. As far as I know, none of the CASes actually implement that (it's much too slow), instead they all implement algorithms which might fail (but with extremely low probability). – Jacques Carette Mar 31 '10 at 21:29
Dror, Newton-Raphson is not guaranteed to converge. Even if it does, it finds only ONE solution. The question asked for ALL solutions. – user1855 Apr 1 '10 at 2:17
Well, from a theoretical perspective, this follows from the decidability of the theory of the real numbers as an ordered field, as proved by Tarski. I agree of course that if you want an efficient algorithm, that's a separate question. – Pete L. Clark Apr 1 '10 at 7:16
@Dror: dividing out by the root you found is known as 'deflation' as is amazingly badly behaved numerically. After you've deflated out about 10 roots, what you're left with is usually a total mess (from an error analysis point of view) and in practice the 'roots' you get after deflation are useless. – Jacques Carette Apr 29 '10 at 1:11

This argument is problematic; see Andrej Bauer's comment below.

Sure. I have no idea what an efficient algorithm looks like, but since you only asked whether it's possible I'll offer a terrible one.

Lemma: Let $f(z) = z^n + a_{n-1} z^{n-1} + ... + a_0$ be a complex polynomial and let $R = \text{max}(1, |a_{n-1}| + ... + |a_0|)$. Then all the roots of $f$ lie in the circle of radius $R$ centered at the origin.

Proof. If $|z| > R$, then $|z|^n > R |z|^{n-1} \ge |a_{n-1} z^{n-1}| + ... + |a_0|$, so by the triangle inequality no such $z$ is a root.

Now subdivide the disk of radius $R$ into, say, a mesh of squares of side length $\epsilon > 0$ and evaluate the polynomial at all the lattice points of the mesh. As the mesh size tends to zero you'll find points that approximate the zeroes to arbitrary accuracy.

There are also lots of specialized algorithms for finding roots of polynomials at the Wikipedia article.

-
To speed up your "Now subdivide..." argument you'd compute some Lipschitz bounds for the polynomial's derivative on your subdivisions and apply Kantorovich's theorem (using Newton's method to find the roots). In practice this is very fast. – Ryan Budney Mar 31 '10 at 23:04
:) It really is terrible... – Dror Speiser Mar 31 '10 at 23:27
A somewhat similar idea, due to Mike Meylan, follows: bound the circle in a square. Then, recursively, subdivide the square into 4 squares. Now approximately compute the argument principle integral around each square, and zoom in on any square that had a value larger than 0. This reduces the complexity dependence on $R$ from quadratic to linear. – Dror Speiser Apr 1 '10 at 14:00
The proposed answer does not work. If the taks is to list all the zeroes then the algorithm must decide which squares of side length $\epsilon$ contain zeroes and which do not. How is it supposed to do that? Just because an approximate value at a point in the mesh is close to zero does not mean there is an actual zero there. What you are proposing is to compute a sequence of nested compact sets (finite unions of squares) whose intersection is the set of zeroes. But the trouble is that some of the squares may "disappear" after a while, so it's hard to tell where the zeroes actuall are. – Andrej Bauer Apr 29 '10 at 11:59
Thanks, Andrej. Do you know if Dror's improvement avoids this problem? – Qiaochu Yuan Apr 29 '10 at 15:31

Homotopy continuation method is good for finding all COMPLEX solutions to arbitrary accuracy, and it is implemented in the Numerical Algebraic Geometry package in Macaulay 2, for example. The method is more general. It can solve a system of polynomial equations in many variables. In fact, it is a more difficult problem to find all REAL solutions WITHOUT finding all complex solutions.

From what I understand, the solution region does not need to be bounded for homotopy continuation to work. You can also "projectify" your problem if necessary, so that you don't have to worry about homotopy paths going off to infinity. Some methods assume that the solutions are all simple, but there're ways to work around it. One is the method of "deflation".

-

For univariate polynomials you should look at "An Efficient Algorithm for the Complex Roots Problem" by Andy Neff and John Reif http://citeseerx.ist.psu.edu/viewdoc/download;jsessionid=5E9156BAF80D8D6AEDCA2F42C11AB4B2?doi=10.1.1.33.3353&rep=rep1&type=pdf

-

The wikipedia article http://en.wikipedia.org/wiki/Root-finding_algorithm gives links to many different methods for finding roots of polynomials. (Start at the section entitled "Finding roots of polynomials".) Many of the methods are incomparable, in the sense that they work faster or slower than others depending on the specific polynomial.

-

One of the semi-recommended ones for finding roots in the complex plane is Laguerre's method, which for some reason is not included in the Wikipedia article on root-finding.

The reason I know of this is a colloquium lecture long ago by Steven Smale on the complexity of Newton's method, during which William Kahan stood up and held forth on why Newton's method was worthless and Laguerre's was much better.

I cannot tell whether you insist on finding all roots to high accuracy. One could perhaps divide out by $(x - r_k)^{n_k}$ each time a root $r_k$ with multiplicity $n_k$ is found, and search for roots for the new polynomial, using those results as seed values for finding accurate roots using the original polynomial.

-
Never bet against Kahan. He has been shown right an astonishing number of times. [Same with David Parnas for software engineering issues]. – Jacques Carette Apr 29 '10 at 1:14

This can be done. Check this article by Hubbard, Schleicher, and Sutherland, entitled "How to find all roots of complex polynomials by Newton's method".

-

Although it's not specific to polynomials with integer coefficients, have a lot at "Computing the Zeros of Analytic Functions".

-

At least for real roots it can be completely solved by bracketing zeroes with Sturm sequences.

-

You have already seen McNamee's excellent bibliography on polynomial root-finding methods?

Personally I have a preference for the "simultaneous iteration" methods (of which Durand-Kerner and Ehrich-Aberth are two of the simplest and most well-known); all you need to start from is a set of points equispaced around a circle in the complex plane (as to the radius of this circle, there are a number of suggestions in the literature; alternatively, formulas in Marden's "Geometry of Polynomials" might be of use here).

-

A completely ineffective theoretical method goes as follows: Write $f(z)$ as $f_r(x,y)+if_i(x,y)$ where $f_r,f_i\in \mathbb R[x,y]$ are real polynomials in the real and complex part $x,y$ of $z=x+i y$. Compute a Groebner basis of the ideal $(f_r,f_i)$ with respect to an order which eliminates one of the variables in the first element of the basis and use real techniques (based on Sturm sequences) to compute, say, the real parts of all solutions. Use another element of the Groebner basis (or again real techniques) to compute the corresponding imaginary part and test for multiplicities (which can be avoided by computing first gcd$(f,f')$).

Completely useless (and equivalent) variation: Study the intersection giving the zeroes of $f$ of the two real curves determined by $f_r$ and by $f_i$.

-