Algorithm for Approximating Complex Polynomial Zeros
Victor Pan
Lehman College, CUNY
Algorithms Seminar
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).
Abstract
An algorithm for approximating complex polynomial zeros is presented. Its
complexity is optimal up to polylogarithmic factors and holds the
current record.
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
following theorem.
Theorem 1 Let p(x) be a monic polynomial of degree n
and z_{1},...,z_{n} its zeros, with z_{i}£1, i=1,...,n. For a
fixed positive b, approximations z_{i}^{*} satisfying

z_{i}z 

<2^{b}, i=1,...,n
(1) 
can be computed at a cost bounded by O^{~}(n) arithmetic
operations and O^{~}(n^{2}(b+n)) boolean operations. The
notation O^{~} means that factors log n, log b or smaller
are neglected.
Much more precise statements, proofs and parallel complexity
estimates can be found in [5] and a pedagogical introduction
to this area is [6].
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^{~}(n^{2}(b+n)) is optimal in
the very frequent case n=O(b).
Actually, O(n^{2}b) is even a lower
bound for the computation of one root of polynomials of
degree n.
This bound follows from the high susceptibility of the
roots of a polynomial with respect to the coefficients. For instance,
the polynomial x^{n}a with a small a>0 has for
root a^{1/n}. If this root is of order 2^{b},
changing a to 0 is a change of the nbth bit of a
coefficient that changes the bth bit of the root. This reasoning
extends to other coefficients: let p=O(n) and consider x^{n}a
x^{p}. Then again a change of a bit at position O(nb) modifies
the bth bit of the solution. Thus b bits of the
solution depend on O(nb) bits of each of O(n) coefficients, whence
the O(n^{2}b) 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 nk 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:

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.
Each of these steps has to be performed in O(n^{2}b+n^{3}) boolean
operations to yield the theorem.
The factors p_{k} and p_{nk} 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
inequality:

p(x)p_{k}(x)p_{nk}(x)£e 

p(x),
(2) 
where q(x) denotes the sum of the
moduli of the coefficients of a polynomial q.
Lemma 1 [8] If
½½ ½½ ½½ ½½ 
p(x) 


½½ ½½ ½½ ½½ 
<ep(x), 
with log_{2}e³ bn+n+2, the inequalities (1)
are satisfied.
Lemma 2 [8] Let p(x), f_{1}(x),...,f_{k}(x)
and f(x),g(x) be polynomials such that

p(x)f_{1}(x)··· f_{k}(x) 

(3) 

f_{1}(x)f(x)g(x) 
£e_{k}f_{1}(x), 
(4) 
then
p(x)f(x)g(x)f_{2}(x)···
f_{k}(x)£e 

p(x) 
holds, provided
From these lemmas follows that it is sufficient to compute the
splitting with e^{*}£e/(n2^{n})
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 p_{k}(x) be the monic polynomial
whose k roots are those roots of p lying inside the circle. The
computation of p_{k}(x) relies on the following integral
representation of the power sums s_{j} of its zeros:
This idea originates in [2] and was refined
by [8] to produce error bounds, i.e., to bound Q such
that the s_{j}'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 [8]. If the closest root to the circle is
at distance O(1/n), a value of Q of order O(n^{2}) is
used^{1} and the
corresponding p'(w^{q}) and p(w^{q}) are computed by a discrete
Fourier transform. From there, the sums s_{j}^{*} for j=0,...,K
are computed by DFT, K being the smallest power of 2 larger
than k=s_{0}. An approximation of the factor p_{k}(x) can then be
recovered efficiently by a
variant of NewtonHensel'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
p(x)p_{k}^{(0)}(x)p_{nk}^{(0)}(x)£e,
where p_{k}^{(0)} has degree k, the aim is to find a
refinement p_{k}^{(1)}=p_{k}^{(0)}+q_{k},
p_{nk}^{(1)}=p_{nk}^{(0)}+q_{nk} with deg q_{i}<i, improving
the error. Since
pp_{k}^{(1)}p_{nk}^{(1)}=(pp_{k}^{(0)}p_{nk}^{(0)})p_{k}^{(1)}p_{nk}^{(0)}
p_{k}^{(0)}p_{nk}^{(1)}p_{k}^{(1)}p_{nk}^{(1)},
the Newton iteration is obtained by satisfying
(
p
p_{k}^{(0)}p_{nk}^{(0)})=
p_{k}^{(1)}p_{nk}^{(0)}+
p_{k}^{(0)}p_{nk}^{(1)},
(5)
which determines p_{k}^{(1)} and p_{nk}^{(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 p_{nk}^{(i)}
modulo p_{k}^{(i)} by a second, parallel, Newton iteration and
then p_{k}^{(i)} is given
by q^{(i)}p=q^{(i)}(pp_{k}^{(i)}p_{nk}^{(i)})mod p_{k}^{(i)}. A
similar formula gives p_{nk}^{(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 [3]). Starting from p(x) of degree n, one performs the
following iteration:
p_{i+1}(x^{2})=(1)^{n}p_{i}(x)p_{i}(x),
which transforms the polynomial p_{i}(x) into a
polynomial p_{i+1}(x) whose roots are the squares of the roots
of p_{i}(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.
Lemma 3
Let z_{1},...,z_{n} be the roots of p(x), satisfying
z_{1}£...£z_{n}£1. Given c>0 and d³0, it is possible
to
compute r_{1},r_{1},...,r_{n},r_{n}
such that r_{k}£z_{k}£r_{k}=(1+c/n^{d})r_{k},
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,1r), with some r<1 such that
z_{k+1}/z_{k}³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'=2r_{11n/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 [5], in particular it
is shown that it is not necessary to compute all the zeros
of p^{(ë n/2û +1)}.
Conclusion
This summary is a very rough sketch of a very detailed study given
in [5]. 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 [4].
References
 [1]

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.
 [2]

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. 543560.
 [3]

Henrici (Peter). 
Applied and computational complex analysis. 
WileyInterscience, New York, 1974, Pure and
Applied Mathematics, vol. Volume 1: Power series, integration, conformal
mapping, location of zeros, xv+682p.
 [4]

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. 488496. 
ACM Press, 1998.
 [5]

Pan (V. Y.). 
Optimal and nearly optimal algorithms for approximating polynomial
zeros. Computers & Mathematics with Applications, vol. 31, n°12,
1996, pp. 97138.
 [6]

Pan (Victor). 
Solving polynomials with computers. American Scientist,
vol. 86, 1998, pp. 6269.
 [7]

Pan (Victor Y.). 
Optimal (up to polylog factors) sequential and parallel algorithms
for approximating complex polynomial zeros. In SIAM Journal on
Computing. pp. 741750. 
ACM Press, 1995.
 [8]

Schönhage (Arnold). 
The fundamental theorem of algebra in terms of computational
complexity. 
Technical report, Mathematisches Institut der Universität
Tübingen, 1982. Preliminary report.
 [9]

Schönhage (Arnold). 
Equation solving in terms of computational complexity. In Proceedings of the International Congress of Mathematicians,
pp. 131153. 
1987. Berkeley, California, 1986.
 1
 More precise values are given
in [8, p. 35].
This document was translated from L^{A}T_{E}X by
H^{E}V^{E}A.