# State of art numerical algorithm for all complex roots of high degree polynomial in $\mathbb Z[x]$

I have many polynomial with high (100+) degree. The coefficients are integer. I would like to find all the roots in $\mathbb C$ of each polynomial as fast as possible. Best complexity and best implementation are needed because I need solve very large set of them.

Note that algorithm needs to return all roots, not just one. (Well actually I would also be interest if there is some algorithm that finds not all but most root, say 80-90%, if enough speed increase to be worth.)

I have looked into several options from wikipedia's list but there are many with differing complexity and for condition other than mine.

What is the state of the art fastest method for solving all root at once of super high degree, integer polynomial?

Second part and very related question.

What is the fastest package or programming language to run above algorithm?

• See the comment at the top of the latest version of the answer below. I'm starting to think that the silly algorithm I made up is actually pretty good... – David C. Ullrich Aug 12 '18 at 17:48

Edit: I posted a much improved version of the algorithm below on MO, together with a Python implementation and empirical evidence that the running time is $O(\log(1/\epsilon))$ to approximate the zeroes within $\epsilon$. State of the art or not, if it's actually $\log(1/\epsilon)$ it seems like it should be good enough; hard to imagine doing much better than that.

Edit: I just made up the algorithm below. A comment leads to web pages showing that sure enough it's something people actually do. No surprise there, hard to see what else you're going to do. But note that looking at a few web sites I noticed that people worry about an important issue that I hadn't thought of: If $p$ has multiple roots (points where $p(z)=p'(z)=0$) the efficiency goes totally to hell. See Multiple roots at the bottom for why this is a problem and what to do about it.

What's below is surely not state-of-the-art (edit: maybe closer than I thought the other day). But it's very simple, and it works; also I suspect it's not too bad, efficiency-wise, so since you have no other answers I thought I might mention it.

At first it seemed like finding all the zeroes was going to involve some moderately sophisticated analysis, like calculating integrals that count zeroes or something. But, no, I realized last night that a very simple "bisection" method works.

Note it seems likely that what's below could be improved, at the cost of adding complexity to the code: Use the bisection below until you're fairly close to all the zeroes and then switch to Newton's method. This would require a little math to devise a test for when to switch - you don't want to switch until you get to a point when Newton is guaranteed to converge. More on that later if I figure it out. Meanwhile, the bisection works like so:

Say $$p(z)=\sum_{j=0}^n a_j z^j$$with $a_n\ne 0$. Let $$R=\max\left(1, |a_n|^{-1}\sum_{j=0}^{n-1}|a_j|\right).$$Then $|z|>R$ implies $p(z)\ne0$, so if you set $$Q_0=[-R,R]\times[-R,R]$$then you know all the zeroes lie in $Q_0$.

Now set $$\delta=\sum_{j=0}^n j|a_j|(\sqrt 2R)^{j-1},$$so that $$|p'(z)|\le\delta\quad(z\in Q_0).$$

Now you maintain a finite collection of squares. At each iteration you refine this collection: Divide each square into four subsquares. "Test" each subsquare, rejecting it or accepting it; the new collection consists of all the accepted subsquares.

To "test" a square $Q$: Say $a$ is the center of $Q$ and $r$ is the distance from $a$ to a corner of $Q$. If $$|p(a)|>r\delta$$it follows that $p$ has no zero in $Q$; in that case reject $Q$, otherwise accept it.

Now if $K_k$ is the union of the squares accepted at the $k$-th stage it follows that $K_k$ converges to the zero set of $p$.

Proof: Say $S=\bigcap_k K_k$. Since each rejected square contains no zeroes of $p$ it follows that the zero set is contained in $S$. Conversely, say $z\in S$. Then there exists a sequence $Q_k$ such that $Q_k$ is one of the squares accepted at the $k$-th stage, $Q_{k+1}\subset Q_k$, and $\{z\}=\bigcap_k Q_k$. Now since $Q_k$ passed the test, $$|p(a_k)|\le\delta r_k\to0;$$since $a_k\to z$ this shows that $p(z)=0$.

Multiple roots: Things blow up in the presence of multiple roots. Why: Imagine first that $p(z)=z$. I believe it's easy to see that at each stage there will be exactly four accepted squares, no problem. But now assume $p(z)=z^2$. The problem is that there will be a huge number of squares where $p$ is very small, although the squares contain no zeroes of $p$. In fact it seems to me that for large $k$, at the $k$-th stage there wiil be $c2^k$ accepted squares - very bad.

Fix: Implement the division algorithm, giving quotient plus remainder for the quotient of two polynomials. Use that with Euclid's algorithm to find $d=gcd(p,p')$. Let $q=p/d$. Then $q$ has the same complex zeroes as $p$, except that all the zeroes of $q$ are simple.

Some good news: At first I thought that the fix wouldn't work in floating-point, so you'd need to implement exact calculations in $\Bbb Q[x]$, which would be a slight pain in an environment with no built-in rational data type. But on second thought I believe that you can do it using floating point (if you do have exact rational arithmetic handy it's not clear to me which would be more efficient.)

The (non?)-problem: Say $\tilde d$ is your floating-point approximation to $d=gcd(p,p')$. Then $p$ is not exactly divisible by $\tilde d$, so you're sunk. No, I think this is ok. The division algorithm gives you $p=\tilde d\tilde q+r$. Now, $r$ would vanish if not for roundoff error, so $r$ will be small. Hence $\tilde q$ gives a good approximation to $q=p/d$, and it follows that the zero set of $q$ is appproximated by the zero set of $\tilde q$ (and also that the zereos of $\tilde q$ are simple, if the roundoff error is small enough).

• Other keywords are Weyl scheme, bisection-exclusion method, inner root radius computation via Dandelin-Graeffe iterations, see for instance perso.math.univ-toulouse.fr/yak/files/2012/04/… – LutzL Aug 8 '18 at 15:05
• Other methods that follow the domain subdivision path are the Lehmer-Schur algorithm and Schönhage's Splitting circle method (see also the referenced papers for inspired methods). – LutzL Aug 9 '18 at 18:18
• There are all-at-once iteration methods like the Durand-Kerner and Aberth-Ehrlich methods. Multiple or close-by roots result in linear, straight paths of the root approximations so that a cluster analysis can detect them and accelerate convergence. There are papers by V. Pan about this and other algorithmic improvements. – LutzL Aug 9 '18 at 18:18