Algorithm for Approximating Complex Polynomial Zeros
Lehman College, CUNY
June 8, 1998
[summary by Bruno Salvy]
A properly typeset version of this document is available in postscript
and in pdf
If some fonts do not look right on your screen, this might be
fixed by configuring your browser (see the documentation here).
An algorithm for approximating complex polynomial zeros is presented. Its
complexity is optimal up to polylogarithmic factors and holds the
Finding roots of a complex polynomial numerically in a guaranteed way
with a fixed prescribed accuracy is difficult when no approximation
is known in advance. This task cannot be performed in a fixed
precision environment and implementations in computer algebra systems
(where arbitrary precision is available) are seldom able to treat
polynomials of degree a few hundreds. However, polynomials of very high
degree arise frequently when solving a polynomial system by elimination.
The work summarized here provides an algorithm supporting the
Theorem 1 Let p(x) be a monic polynomial of degree n
and z1,...,zn its zeros, with |zi|£1, i=1,...,n. For a
fixed positive b, approximations zi* satisfying
Much more precise statements, proofs and parallel complexity
estimates can be found in  and a pedagogical introduction
to this area is .
can be computed at a cost bounded by O~(n) arithmetic
operations and O~(n2(b+n)) boolean operations. The
notation O~ means that factors log n, log b or smaller
The statement of the theorem can be modified to accommodate polynomials
which are not
monic (by first scaling the coefficients) or with roots of modulus
larger than 1 by computing a bound on the moduli (see below) and then
scaling the polynomial.
1 Lower Bounds
It is clear that the arithmetical complexity O~(n) is
optimal, since n coefficients of the input polynomial have to be
treated. The boolean complexity O~(n2(b+n)) is optimal in
the very frequent case n=O(b).
Actually, O(n2b) is even a lower
bound for the computation of one root of polynomials of
This bound follows from the high susceptibility of the
roots of a polynomial with respect to the coefficients. For instance,
the polynomial xn-a with a small a>0 has for
root a1/n. If this root is of order 2-b,
changing a to 0 is a change of the nb-th bit of a
coefficient that changes the b-th bit of the root. This reasoning
extends to other coefficients: let p=O(n) and consider xn-a
xp. Then again a change of a bit at position O(nb) modifies
the b-th bit of the solution. Thus b bits of the
solution depend on O(nb) bits of each of O(n) coefficients, whence
the O(n2b) lower bound. This example also illustrates why clusters
of zeros defeat many numerical algorithms.
2 Outline of the Algorithm
The algorithm is based on a splitting technique where
the polynomial p is split into factors of
degree k and n-k with k=a n, for
some aÎ(1/2,r), r being fixed. Applying this process
recursively, any polynomial can be
completely factored in O(log n) steps.
The splitting itself is computed in 3 steps:
Each of these steps has to be performed in O(n2b+n3) boolean
operations to yield the theorem.
Find a ``splitting'' circle not ``too close'' to roots of p
and containing a n of them;
- Compute the polynomial vanishing at these a n roots;
- Divide p by this polynomial to obtain the other factor.
The factors pk and pn-k of p are computed numerically. The
following two lemmas show how the
precision with which they are required can be bounded by ensuring
that e* is sufficiently small in the following
where ||q(x)|| denotes the sum of the
moduli of the coefficients of a polynomial q.
|Lemma 1  If
with -log2e³ bn+n+2, the inequalities (1)
|Lemma 2  Let p(x), f1(x),...,fk(x)
and f(x),g(x) be polynomials such that
From these lemmas follows that it is sufficient to compute the
splitting with e*£e/(n2n)
in (2), where e comes from Lemma 1.
The splitting circle method was introduced by
Schönhage [8, 9].
We now review the algorithms used in steps 1 and 2, together with the
recent progress due to Victor Pan.
3 Numerical Factorization
To simplify the notation, assume the unit circle is a splitting circle
for the polynomial p(x). Let pk(x) be the monic polynomial
whose k roots are those roots of p lying inside the circle. The
computation of pk(x) relies on the following integral
representation of the power sums sj of its zeros:
This idea originates in  and was refined
by  to produce error bounds, i.e., to bound Q such
that the sj's can be computed by the discretization
The value of Q depends on a lower bound for |p(z)| on the unit
circle, which in turns is related to a bound on the distance from this
circle to the closest root of p, hence the need for a circle ``not
too close'' to the roots in Step 1 of the algorithm.
Efficiency is attained at the price of quite technical
developments . If the closest root to the circle is
at distance O(1/n), a value of Q of order O(n2) is
used1 and the
corresponding p'(wq) and p(wq) are computed by a discrete
Fourier transform. From there, the sums sj* for j=0,...,K
are computed by DFT, K being the smallest power of 2 larger
than k=s0. An approximation of the factor pk(x) can then be
recovered efficiently by a
variant of Newton-Hensel's lifting (see [1, p. 34]). Then
the other factor is obtained by division. In order to reach the right
level of complexity, it is necessary to compute only O(n) bits
for these steps and then refine the factorization by another Newton
like algorithm as follows. Starting from the approximate factorization
where pk(0) has degree k, the aim is to find a
pn-k(1)=pn-k(0)+qn-k with deg qi<i, improving
the error. Since
the Newton iteration is obtained by satisfying
which determines pk(1) and pn-k(1) uniquely. These
polynomials could be found by Euclid's algorithm, but this is too
expensive. Instead, one also computes an inverse q(i) of pn-k(i)
modulo pk(i) by a second, parallel, Newton iteration and
then pk(i) is given
by q(i)p=q(i)(p-pk(i)pn-k(i))mod pk(i). A
similar formula gives pn-k(i). Then the required precision is
obtained after a few iteration at a cost bounded by O(nloge*).
4 Finding Splitting Circles
The basic technique to find discs containing a known number of roots
of a polynomial is the iteration of Graeffe's
method (see ). Starting from p(x) of degree n, one performs the
which transforms the polynomial pi(x) into a
polynomial pi+1(x) whose roots are the squares of the roots
of pi(x). This process emphasizes the differences of moduli between
the roots. The coefficients of these
iterates are Newton sums from which precise information about the
different moduli of the roots of the original polynomial can be
recovered at a low cost. More precisely, one gets the following lemma.
Let z1,...,zn be the roots of p(x), satisfying
|z1|£...£|zn|£1. Given c>0 and d³0, it is possible
such that rk£|zk|£rk=(1+c/nd)rk,
k=1,...,n with O~(n) arithmetic operations.
This iteration is applied after having first shifted the origin to
the center of gravity of the roots, which is given by the first two
coefficients of the polynomial. When it follows from this computation
that there is a k=a n, a in a fixed interval
(r,1-r), with some r<1 such that
|zk+1|/|zk|³1+c/n for some c fixed in advance, then this
yields a splitting circle and the factoring algorithm of the previous
section can be applied.
It is when no such circle can be found that progress has been made by
Victor Pan recently. In this case, there is an annulus centered at 0 which contains
most of the roots of the polynomial. Now the idea is to
shift the origin to each of r'=2r11n/12 and ir',
and apply the same method. Then either
a good splitting circle is found, or there is a small circle which is
easily computed and contains the intersection of
these three annuli, itself containing an important cluster of zeros
(at least half of the zeros of p if c=1/100). In this case, the
idea is that one of the zeros of a derivative of p of high order
(for instance, one can take p(ë n/2û +1)) is either the
center of a good splitting circle or makes it possible to isolate a
massive cluster of zeros, where more than half of the zeros
of p are at distance less than the desired accuracy 2-b. In
both cases, the polynomial can then be factored numerically and the
computation proceeds on those factors that do not correspond to a
massive cluster. Many refinements are given in , in particular it
is shown that it is not necessary to compute all the zeros
of p(ë n/2û +1).
This summary is a very rough sketch of a very detailed study given
in . For practical polynomial solving, other algorithms
are known to perform extremely well, but their complexity analysis has
yet to be done.
The talk also mentioned extensions to the multivariate case, this is
described in .
Bini (Dario) and Pan (Victor Y.). --
Polynomial and matrix computations. Vol. 1. --
Birkhäuser Boston Inc., Boston, MA, 1994, Progress in Theoretical Computer Science, xvi+415p. Fundamental algorithms.
Delves (L. M.) and Lyness (J. N.). --
A numerical method for locating the zeros of an analytic function.
Mathematics of Computation, vol. 21, 1967, pp. 543--560.
Henrici (Peter). --
Applied and computational complex analysis. --
Wiley-Interscience, New York, 1974, Pure and
Applied Mathematics, vol. Volume 1: Power series, integration, conformal
mapping, location of zeros, xv+682p.
Mourrain (Bernard) and Pan (Victor Y.). --
Asymptotic acceleration of solving multivariate polynomial systems of
equations. In Proceedings of the Thirtieth Annual ACM Symposium on
Theory of Computing. pp. 488--496. --
ACM Press, 1998.
Pan (V. Y.). --
Optimal and nearly optimal algorithms for approximating polynomial
zeros. Computers & Mathematics with Applications, vol. 31, n°12,
1996, pp. 97--138.
Pan (Victor). --
Solving polynomials with computers. American Scientist,
vol. 86, 1998, pp. 62--69.
Pan (Victor Y.). --
Optimal (up to polylog factors) sequential and parallel algorithms
for approximating complex polynomial zeros. In SIAM Journal on
Computing. pp. 741--750. --
ACM Press, 1995.
Schönhage (Arnold). --
The fundamental theorem of algebra in terms of computational
Technical report, Mathematisches Institut der Universität
Tübingen, 1982. Preliminary report.
Schönhage (Arnold). --
Equation solving in terms of computational complexity. In Proceedings of the International Congress of Mathematicians,
pp. 131--153. --
1987. Berkeley, California, 1986.
- More precise values are given
in [8, p. 35].
This document was translated from LATEX by