On the Geometry of Graeffe Iteration

Gregorio Malajovich

Dep. de Matemática Aplicada, Universidade Federal do Rio de Janeiro

Caixa Postal 68530, Rio de Janeiro, RJ, 21945 – BRASIL

gregorio@labma.ufrj.br, http://www.labma.ufrj.br/˜ gregorio

Jorge P. Zubelli

IMPA

Estrada Dona Castorina 110, Rio de Janeiro, RJ, 22460-320, BRASIL

zubelli@impa.br, http://www.impa.br/˜ zubelli

Second revision, Aug 26, 1999

Abstract

A new version of the Graeffe algorithm for finding all the roots of univariate complex polynomials is proposed. It is obtained from the classical algorithm by a process analogous to renormalization of dynamical systems.

This iteration is called Renormalized Graeffe Iteration. It is globally convergent, with probability 1. All quantities involved in the computation are bounded, once the initial polynomial is given (with probability 1). This implies remarkable stability properties for the new algorithm, thus overcoming known limitations of the classical Graeffe algorithm.

If we start with a degree-d polynomial, each renormalized Graeffe iteration costs O(d^{2}) arithmetic operations, with memory O(d).

A probabilistic global complexity bound is given. The case of univariate real polynomials is briefly discussed.

A numerical implementation of the algorithm presented herein allowed us to solve random polynomials of degree up to 1000.

Contents

1 Introduction 3 2 Iterative Algorithms and Renormalization 7 3 Recurrence Relations and the Renormalization of Graeffe 17 4 Probability of Success 23 5 Newton Diagram Revisited 28 6 Numerical Results 30 7 Acknowledgments 32

1 Introduction

Graeffe Iteration was one of the most prestigious XIX century algorithm for finding roots of polynomials. At that time, computations were performed by hand, by people payed specifically to perform those computations. They were called calculateurs [34] or computers [21].

Let f be a univariate polynomial, of degree d. Its Graeffe iterate is defined as:

Gf(x) = (−1)^{d}f(\sqrt{x})f(−\sqrt{x})  (1)

This defines a many-to-one mapping in the space of all degree-d polynomials (real or complex, as wish). The effect of this mapping is to square each root of f.

After a few Graeffe iterations, the roots of G^{k}f have (hopefully) incommensurate moduli. This is not true for complex-conjugate roots, which can be worked out in a different way.

Assume, for simplicity, that f is a complex polynomial, with no two roots of the same modulus. Then, G^{k}f can be written

G^{k}f(x) = ∑_{i=0}^{d} a_{i}^{(k)}x^{i}  (2)

where a_{i}^{(k)} is given by the (d − i)-th symmetric function of the roots of G^{k}f. Therefore, a_{i}^{(k)} is dominated by:

(−1)^{d−i}ζ_{1}^{2^{k}}ζ_{2}^{2^{k}} … ζ_{d−i}^{2^{k}}  (3)

where ζ_{1}, …, ζ_{d} are the roots of f ordered with decreasing modulus. (A more rigorous statement of this will appear in Section 5)

Therefore, −a_{i−1}^{(k)}/a_{i}^{(k)} is a good approximation for ζ_{d−i+1}^{2^{k}}.

Hence it is computationally easy to approximate |ζ_{i}| for all i. Although we also obtain arg ζ_{i} mod 2^{1−k}π, we will discard this information in this paper to avoid additional complications. There are many classical algorithms to recover the actual value of ζ_{i}. See Pan [38] for a discussion. ^{1}

In this note, we apply a suitable non-uniform change of coordinates (renormalization) to the Graeffe iteration operator, to make it ‘convergent’ with probability 1.

1Added in the revised version: We deal with this issue in [28].

The algorithm obtained by this change of coordinates will be called Renormalized Graeffe Iteration. We use the following systems of coordinates for each iterate G^{k}f of the Graeffe method applied to f: The coefficients of G^{k}(f) and all related intermediate computations will be represented in scaled polar coordinates, where a complex number w is represented by ‘magnitude’ 2^{−k} log_{2} |w| and ‘argument’ arg w ∈ [−π, π]. Calculations will always be performed ‘in coordinates’. The ‘magnitude’ variables of G^{k}f ‘in coordinates’ will converge with probability 1 (Theorem 1 below).

The precise construction of the Renormalized Graeffe Operator is postponed to Section 2.

We also claim that the Renormalized Graeffe algorithm compares well with available numerical software or theoretical algorithms.

However, in this paper, we are considering a modified problem: our algorithm is designed to find the absolute values of the roots, not the actual roots. Therefore we will compare its complexity to the complexity of finding the absolute value of the roots by other existing algorithms.

In [28], we explain how to modify this algorithm to obtain the actual roots, without endangering the complexity estimates.

The algorithm presented here has arithmetic complexity O(d^{2}) for each iteration, and memory size O(d), where d is the degree of the polynomial.

The number of iterations will be bounded also in terms of a probability of failure (Theorem 2 below). The authors personally believe that this is only possible due to the clean mathematical structure of Graeffe iteration. Our bound improves previous probabilistic bounds (in the sense of probability of success) on the complexity of solving polynomials. (See Renegar [41] and Shub-Smale [45]).

Also, our algorithm compares well with practical software, like for instance the algorithm in Matlab (running time O(d^{3}) and memory O(d^{2})). Instead, Renormalized Graeffe Iteration seems to run in time O(d^{2}). It also seems much more stable, for most (in a probabilistic sense) complex polynomials of degree, say, 1000. (See Section 6 for a discussion of preliminary experimental results).

In order to compare with deterministic algorithms one should bear in mind that the algorithm presented here is probabilistic and may be quite slow on a set of non-zero measure. Also, the complexity of deterministic algorithms such as [32, 38, 48, 49] can be given in terms of number of arithmetic operations or in number of bit operations. We believe that estimates of the order of O(d^{1+α}), α > 0 in the number of operations are made at the expense of a large increase in the working precision, of the order of O(d^{2+α}) bits. However, this is not a lower bound and it may as well happen that those algorithms can work with substantially smaller precision on a large set of inputs.

In comparison, our Renormalized Graeffe Iteration was designed for floatingpoint arithmetic (our results assume an arbitrary, but fixed floating-point precision). However, the details of the rounding-off analysis of renormalized Graeffe iteration will not be discussed in this paper.

Instead, some preliminary numerical results are presented in Section 6. They support the empirical fact that typical random degree 1000 polynomials can be solved within a precision of 64 bits of mantissa (IEEE 854 long double). We expect a factor of a polynomial in log d bits of mantissa to be sufficient in general, with probability 1. (This is a conjecture).

The result in Theorem 1 holds for “reasonable” probability distributions on the space of polynomials. By reasonable, we mean all probabilities with bounded Radon-Nikodym derivative with respect to Lebesgue probability in the projectivization of coefficient space.

Theorem 1. There is a renormalization of the Graeffe iteration, such that if f is a degree-d polynomial (in a measure theoretical sense) then with probability 1 this renormalized Graeffe iteration produces d + 1 sequences, each one converging to some h_{i}, s.t. log |ζ_{i}| = h_{i} − h_{i+1}, and ζ_{1}, ⋯, ζ_{d} are roots of f. Moreover, each iteration can be performed in O(d^{2}) arithmetical operations and all iterations can be performed with memory O(d).

This theorem is constructive, in the sense that an explicit construction of the renormalized Graeffe iteration will be given.

In Section 2, we discuss the precise meaning of renormalization in our context. Its main consequence, will be to produce an algorithm operating on a bounded set of numbers. This solves the main stability problem of classical Graeffe iteration that prevented it from finding all roots at once. See for example Henrici’s comments on FFT-based Graeffe iteration [18] (Vol III, last paragraph of p. 69). The definition of renormalization will outlaw FFT-based Graeffe iteration. Indeed, although FFT is known to be stable with respect to vector norms, it is not component-wise stable with respect to the relative error. This means that some of the coefficients of the k-th Graeffe iterate may have a large relative error, and hence some of the roots will be extremely inaccurate. This may be disastrous if one wants to retrieve all the roots at the same time.

Theorem 2. Let f be a random complex polynomial of degree d. Let b ≥ 1 + log_{2} d. Then, with probability 1 − δ, k steps of the Renormalized Graeffe Iteration will approximate the log |ζ_{i}|’s with relative precision 2^{−b}, where

k ≥ c_{1} + c_{2} log_{2} b − c_{3} log_{2} δ

and c_{2}, c_{3} are universal constants. The constant c_{1} depends on the choice of the probability distribution, and on d.

Whenever speaking of random polynomials, we like to consider the normally invariant probability density introduced by Kostlan [25] (See also Section 4). However, the above mentioned result is true for any reasonable probability distribution.

The experimental results in figure 6 support the conjecture that, under Kostlan’s probability distribution, we can fix

c_{1} = c_{4} log_{2} d

where c_{4} ≈ 2 is a universal constant.

We will briefly discuss the real case, and how to deal with complexconjugate roots or roots with same modulus in Section 5.

Historical remarks. The Graeffe iteration was developed independently by Dandelin (1826), by Graeffe (1837) and by Lobachevsky (1834). We call it Graeffe iteration to conform with most of the literature. See Householder [20] for early references and priority questions. See Dedieu [9] for an application of Graeffe ’s algorithm.

Important theoretical results were obtained by Ostrowskii [34] in 1940. Also, by that time, numerical analysis books mentioned Graeffe iteration as the preferred algorithm for zero-finding (see e.g. Uspensky [51] Page 318. For another early computer implemented algorithm, see Bareiss [1, 2] and also Blish and Curry [6]).

With the advent of digital computing, the practical use of Graeffe iteration seems to have been forgotten.

Most popular zero-finding algorithms seem to be based now on QR iteration (Matlab) or in a several steps, root-finding plus deflation scheme. (e.g.

Jenkins and Traub [22]). See, however, Cardinal [8], Edelman and Murakami [13], Emiris, Galligo and Lombardi [14], and Toh and Trefethen [50].

In a more theoretical perspective, Graeffe iteration is considered as a sort of pre-conditioning for polynomial splitting. Splitting a polynomial means factorizing it into one factor with large roots, and another with small roots. Splitting is used to obtain extremely fast theoretical algorithms (see Schönhage [48, 49], Kirrinnis [23], Neff and Reif [32], Bini and Pan [5], Mourrain and Pan [31], Pan [36, 37, 38], Pan et al. [39], Malajovich and Zubelli [27]). The main practical difficulty for those algorithms seems to be the large precision required by Graeffe iteration. Also, those algorithms are quite close to known lower bounds on topological complexity (See Vassiliev [52]). For a related lower bound see Novak and Woźniakowski [33].

An important paper by Grau in 1963 [17] laid some of the bases for a version of Graeffe iteration adapted to digital computers. He identified the problem of the increasing numerical range. During Graeffe iteration, some of the coefficients can become so large that the floating point system cannot accommodate them anymore.

While most of the literature suggests to find one root at a time and then use deflation, disregarding some stability problems (See e.g. Henrici [18]), Grau proposed a globally convergent algorithm. Grau’s algorithm would involve only bounded quantities.

As far as we know, that paper was completely forgotten. The algorithm suggested by Grau has complexity O(d^{2}) and memory usage of O(d^{2}). It may be considered as the precursor of the one we shall introduce below.

2 Iterative Algorithms and Renormalization

In this paper, we will produce a version of Graeffe iteration that has bounded numerical range, for most input polynomials. The crucial concept in the construction of this algorithm is the idea of renormalization.

Renormalization is a tool used in understanding the qualitative behavior of iterative phenomena that range over different scales. A rich theory of renormalization exists for one-dimensional dynamical systems. See Feigenbaum [15], McMullen [29], and De Melo-Strien [30]. As for the multidimensional case see Palis-Takens [35].

We shall illustrate what we mean by renormalization in the setting of iterative algorithms by a well-known example and then proceed with the definitions

Example 1. Let f : ℂ → ℂ be a smooth function. Newton iteration associates to a point x the point N(x) = x − f(x)/f′(x). The sequence of Newton iterates converges to a zero of f, provided that x is picked close enough to a non-degenerate zero.

We assume now that f(x) = 0 + f′(0)x + h.o.t. , and that f′(0) ≠ 0. We will consider now the renormalization of the iterative algorithm (or mapping) x ↦ N(x). Although we are dealing here with a simple example (the answer is always 0), its renormalization will have some of the main features of the renormalized Graeffe iteration, yet to be defined.

The basic idea, therefore, is to look at Newton iteration with a variable microscope of variable lens. We look at the mapping:

N_{ϵ} = h_{ϵ^{−2}} ∘ N ∘ h_{ϵ}  (4)

where h_{ϵ} means the homothety of ratio ϵ.

When ϵ tends to zero, N_{ϵ} tends to the map: y ↦ γy^{2}, defined from the disk D(|γ^{−1}|) of radius |γ^{−1}| onto itself. Here, γ is a parameter chosen equal to f^{(2)}(0)/2f′(0). (Proof of this fact: Taylor’s Theorem. The choice of the radius |γ|^{−1} makes the limiting map surjective).

This process (which we call renormalization) gives us qualitative information on the dynamics of Newton iteration near a root. We can summarize this information in an eventually commutative diagram (commutative in the limit):

\begin{aligned}D(|γ^{−1}|ϵ) & ∋ & x & ⟶^{(h_{ϵ})^{−1}} & y & ∈ & D(|γ^{−1}|) \\ ^{N}↓ & & ↓^{y↦γy^{2}} \\ D(|γ^{−1}|ϵ^{−2}) & ∋ & N(x) & ⟶^{(h_{ϵ})^{−2}} & γy^{2} & ∈ & D(|γ^{−1}|) \end{aligned}  (5)

In the example above, the homothety (h_{ϵ})^{−1} is generating the renormalization group. In general, we want to consider renormalizations that lead to a commutative diagram, on the limit:

\begin{matrix}y & ⟶^{R} & R(y) \\ ^{G}↓ & & ↓_{G^{R}} \\ G(y) & ⟶_{R^{2}} & R^{2}(G(y)) \end{matrix}  (6)

Above, G is the original algorithm (possibly after some uniform change of coordinates). We denote by R the renormalization map, and by G^{R} the renormalized version of our iteration. Usually, G^{R} depends on a renormalization parameter. However, we want G^{R} to converge to a limiting map, with a simple dynamics. Moreover, computation of G^{R} should be ‘stable’, in a precise sense.

This suggests the following heuristics, in the case of Graeffe iteration: we would like to consider a polynomial f represented by its roots. More precisely, a polynomial f can be represented by the vector (log |ζ_{i}|, arg(ζ_{i})) ∈ ℝ^{d} × 𝕋^{d}, where ζ_{i} is the i-th root of f, and 𝕋 denotes the additive group ℝ/2πℤ. In that case, Graeffe iteration is just multiplication by 2.

Renormalization would be a division by 2 of the log of the radii of the roots. We would have:

ℝ^{d} × 𝕋^{d}∋  (log |ζ_{i}|, arg ζ_{i})  ⟶^{R^{k}}  (r, θ)  ∈ℝ^{d} × 𝕋^{d}
×2↓  ↓^{(r,θ)↦(r,2θ)}
ℝ^{d} × 𝕋^{d}∋(2 log |ζ_{i}|, 2 arg ζ_{i} mod 2π) ⟶_{R^{k+1}} (r, 2θ mod 2π)∈ℝ^{d} × 𝕋^{d}

(7)

Therefore, the limit of the renormalized map should be the map: (r, θ) ↦ (r, 2θ). Of course, we do not know in advance the roots of the polynomial. Instead, we will produce a chain of commutative diagrams ‘converging’ to diagram (7). In order to do that, let us assume that a polynomial

a_{0} + a_{1}x + a_{2}x^{2} + ⋯ + a_{d}x^{d}  (8)

is represented by the vector

(â_{i}, α_{i}) = (log |a_{i}| , arg(a_{i}))  (9)

In this paper, for clarity of exposition, we will ignore the case where some of the a_{i} is zero. We can do that because most of our results hold ‘with probability 1’. In practice, polynomials with a zero coefficient do arise. See Section 6 for implementation comments.

If the operator G represents Graeffe iteration in this new system of coordinates, we will have:

\begin{aligned}& ↓ & & ↓ \\ ℝ^{d+1} × 𝕋^{d+1} & ∋ & G^{k}(f) & ⟶^{R^{k}} & (â^{k}, α^{k}) & ∈ & ℝ^{d+1} × 𝕋^{d+1} \\ & ^{G}↓ & & ↓^{(G_{k})} \\ ℝ^{d+1} × 𝕋^{d+1} & ∋ & G^{k+1}(f) & ⟶_{R^{k+1}} & (â^{k+1}, α^{k+1}) & ∈ & ℝ^{d+1} × 𝕋^{d+1} \\ & ↓ & & ↓ \end{aligned}  (10)

Above, â^{k} = (â_{0}^{k}, ⋯, â_{d}^{k}) and α^{k} = (α_{0}^{k}, ⋯, α_{d}^{k}). Also, k is a superscript, not an exponent. However, R^{k} denotes the k-th iterate of R.

In this diagram, R maps (r, θ) into (r/2, θ). Although this defines G_{k} as a mapping, the algorithmic construction of G_{k} is postponed to Section 3. On the limit, G_{k} ‘converges’ to the mapping: (a, θ) ↦ (a, 2θ). This will ensure that lim â^{k} exists. We will show that this limit satisfies with probability 1

lim_{k→∞} â_{j}^{k} − lim_{k→∞} â_{j+1}^{k} = log |ζ_{d−j}|

where ζ_{d−j} is the (d − j)-th root of the original polynomial, the roots being ordered by decreasing modulus..

Also, by using the renormalized Graeffe iteration G_{k} instead of the classical Graeffe iteration G, we will be able to bound all the intermediate calculations to a compact, depending on the input f. This will imply several stability results.

In order to define formally what we mean by a renormalized algorithm, we need to state the conditions that we expect G_{k} to satisfy. These conditions will define a class of algorithms, that we call renormalized iterative algorithms. Some examples: usual Newton, renormalized Newton, all reasonable algorithms based on a contraction principle and, of course, Renormalized Graeffe Iteration (To be constructed).

Before defining what we mean by a renormalized algorithm, we will briefly discuss the notion of algorithm. There are many possible definitions (See Blum, Cucker, Shub and Smale [7]).

Definition 1. An iterative algorithm M is a Blum-Shub-Smale (BSS) machine over ℝ, modified as follows:

1. If the input is in ℝ^{l}×𝕋^{m} for a pair (l, m) , then the output is in ℝ^{l}×𝕋^{m}. The integer l + m is called the input size. If the input is denoted by x, the output is denoted M(x), and M^{k} means the composition M ∘ M ∘ ⋯ ∘M

Also, we say that the iterative algorithm ‘computes’ the function x ↦ lim_{k→∞} M^{k}(x)

2. Computation nodes are allowed to perform the following elementary functions: polynomial evaluation, absolute value, (real) logarithms, (real) exponential, sine, cosine, rational functions, and compositions of those.

3. M_{ϵ}(x) will denote the result of approximating M(x) by allowing each operation with non-integer parameters to be performed with relative precision ϵ. (This is also known in numerical analysis as the (1 + ϵ) property.) The parameter ϵ is allowed to vary in (0, \frac{1}{2}). The machine is supposed to “know” ϵ. This means that the value of ϵ can be used in intermediate calculations. Also, the approximation is performed by some prescribed algorithm. (e.g., IEEE arithmetic with − log_{2} ϵ bits of mantissa).

4. For all ϵ, the arithmetic complexity of M applied to the input x is the number of elementary function evaluations (from the list of item 1) and branchings performed with input x. We require the arithmetic complexity of the approximate algorithm M_{ϵ} applied to the input x to be the same as the arithmetic complexity of the original algorithm M applied to input x.

Remark 1. Item 1 will allow us to distinguish between real and angular variables, the latter one being defined modulo 2π. This will simplify notation when we speak of the distance between points in ℝ^{l} × 𝕋^{m}. Let d be that distance.

Remark 2. It has been argued that the outcome of a branching node would not be well-defined in the presence of numerical error. This is not true in the definition above. The branching nodes of the machine M can be assumed, without loss of generality, to branch on queries of the form y > 0 or y ≥ 0 or y = 0. When the machine M gets replaced by M_{ϵ}, the branching nodes stay formally the same. The value of y, however, is contaminated with a certain numerical error. It is still a perfectly defined real number, and it may be compared to zero. Thus, the branching nodes branch ‘correctly’, for a slightly perturbed input.

This would lead to disastrous results if an approximate machine enters a ‘loop’ due to numerical errors. Item 4 in the definition above is there to preclude that sort of loop, by ensuring that the arithmetic complexity of each iteration does not depend on the working precision.

Definition 2. A function φ can be computed in finite time if and only if there is a BSS machine over ℝ modified as above items 2, 3 and 4, that computes φ.

One consequence of the previous definition is the following: suppose a certain function φ can be computed in finite time. Then, its branching set is given by a finite set of equations. Outside the branching set, φ can be written (locally) as a composition of elementary functions. Moreover, only a finite number of such compositions may appear, one corresponding to each set of possible branchings.

A few definitions are in order now. In the sequel we will need to use algorithms depending effectively on a parameter k ∈ ℕ. Those will be given by a machine M with two inputs, say k and f. However, we will denote the output as M_{k}(f) and we will write M_{k} for each of the inputoutput mappings obtained by restricting that machine to some fixed value of k. Also, in that situation we will speak explicitly of the algorithm (M_{k})_{k∈ℕ} or M_{k} for short.

The sequence (M_{k})_{k=1}^{∞} can be considered as a sequence of mappings from ℝ^{l} × 𝕋^{m} into itself. The orbit of f by the sequence (M_{k})_{k=1}^{∞} is the set

orbf =^{def} (f, M_{1}(f), M_{2} ∘ M_{1}(f), …) ⊂ ℝ^{l} × 𝕋^{m} .

This set should be understood as the orbit of f by the non-autonomous dynamical system M_{k} (Not the semi-group !) The closure of orbf will be denoted by \overline{orbf}.

We shall say that a subset of ℝ^{l} × 𝕋^{m} has full-measure if its complement is contained in a set of null measure. We say that some property is true almost everywhere (a.e.) if that property is true in a full-measure set. We can now define a renormalized algorithm and make precise our concept of renormalization:

Definition 3. The algorithm (M_{k})_{k∈ℕ} (or M_{k} for short) is said to be a renormalized iterative algorithm to compute φ : ℝ^{l} × 𝕋^{m} → ℝ^{s}, f ↦ φ(f), where 0 ≤ s ≤ l, if, and only if, it satisfies the Axioms 1 through 4 below.

Axiom 1 (Consistency). For almost every f ∈ ℝ^{l} × 𝕋^{m} we have

lim_{k→∞} (π ∘ M_{k} ∘ M_{k−1} ∘ ⋯ ∘ M_{1}) (f) = φ(f)

where π is the projection of ℝ^{l} × 𝕋^{m} onto the first s coordinates of ℝ^{l}.

Axiom 2 (Arithmetic Complexity). The arithmetic complexity of M_{k} with input f is bounded in terms of the size l + m of f, independently of k and the coefficients of f.

Axiom 3 (Propagation). For almost every f ∈ ℝ^{l} × 𝕋^{m}, there exists a compact neighborhood V ⊂ ℝ^{l} × 𝕋^{m} of orbf (under {M_{k}}) and C such that M_{k}|_{V} is Lipschitz with constant C_{k} and eventually C_{k} < C.

Axiom 4 (Stability). For almost every f ∈ ℝ^{l} × 𝕋^{m}, there exists a compact neighborhood V ⊂ ℝ^{l} × 𝕋^{m} of orbf (under {M_{k}}) and B ∈ ℝ such that ∀g ∈ V and ∀k,

d (M_{k,ϵ}(g) , M_{k}(g)) < ϵB .

Our concept of renormalized iterative algorithm subsumes several ‘reasonable’ properties of iterative algorithms. Axiom 1 allows our algorithm to carry more information than what is actually required at output. Yet, we want our algorithm to produce a sequence converging to the expected result, for almost every input. Axioms 3 and 4 will rule out unstable algorithms. The idea behind Axiom 2 is that any honest iterative algorithm should have bounded arithmetic complexity for each iteration. This prevents the use of multiple precision arithmetic to obtain stability at the expense of (possibly exponentially many) extra arithmetic operations.

We shall now explore some consequences of the definition of renormalized iterative algorithm we proposed above. More specifically, our first goal is to obtain a (non-uniform) error bound on the result of iterating k times a renormalized algorithm with precision ϵ.

Lemma 1. Let M be a renormalized iterative algorithm . Then, for almost every f and for each k, we have that for sufficiently small ϵ > 0

d (M_{k,ϵ} ∘ ⋯ ∘ M_{1,ϵ}(f), M_{k} ∘ ⋯ ∘ M_{1}(f)) < kA^{k}Bϵ .

Here, A and B, depend on f, but not on k. In particular, that will be true for values of ϵ of the form

ϵ ≤ \frac{ρ}{kA^{k}B} ,

where ρ is defined below.

As an immediate consequence of Lemma 1 we have that

∥ (π ∘ M_{k,ϵ} ∘ ⋯ ∘ M_{1,ϵ})(f) − (π ∘ M_{k} ∘ ⋯ ∘ M_{1})(f)) ∥_{2} < kA^{k}Bϵ .

Proof of Lemma 1. We now describe the construction of the constants in the statement of Lemma 1. Combining Axioms 3 and 4, there exists a compact neighborhood W of orbf such that:

1. Every mapping M_{k} is Lipschitz on W. Moreover, there exists a constant A such that for every k, the norm of the Lipschitz constant of M_{k} is uniformly bounded by A = max (1, sup_{k} C_{k}), C_{k} as in Axiom 3.

2. There exists ϵ_{0} such that ∀ϵ < ϵ_{0} we have

d(M_{k,ϵ}(g), M_{k}(g)) < ϵB , ∀g ∈ W .

We then define ρ as the minimum of ϵ_{0} and the distance of orbf to the boundary of W.

We may now conclude the proof of Lemma 1. Set

g = (M_{k−1} ∘ ⋯ ∘ M_{1})(f) ,

and

h = (M_{k−1,ϵ} ∘ ⋯ ∘ M_{1,ϵ})(f) .

In that case, we want to bound

d(M_{k,ϵ}(h), M_{k}(g)) ≤ d(M_{k,ϵ}(h), M_{k}(h)) + d(M_{k}(h), M_{k}(g))

If h ∈ W, one can bound

d(M_{k}(h), M_{k,ϵ}(h)) ≤ Bϵ

and

d(M_{k}(h), M_{k}(g)) ≤ A d(g, h) .

By induction,

d(g, h) ≤ (k − 1)A^{k−1}Bϵ ,

and hence

d(M_{k,ϵ}(h), M_{k}(g)) ≤ kA^{k}Bϵ .

The condition

ϵ < \frac{ρ}{kA^{k}B}

guarantees that h ∈ W.

If we use a Turing machine model (or any other classical discrete complexity model), we can perform all the operations of M in finite precision, and obtain:

Lemma 2. Let M be a renormalized iterative algorithm. For almost every f, assume that the truncation error is bounded by

∥ (π ∘ M_{k} ∘ ⋯ ∘ M_{1}) (f) − φ(f)∥_{2} < E^{2^{k}} ,

where E = E(f) ∈ (0, 1). Then, the complexity of approximating φ(f) with precision δ is O((log_{2} \frac{1}{δ})^{1+α}) where α > 0 is arbitrarily small.

It is assumed above that the cost of arithmetic with l bits of mantissa is O(l^{1+α}), for all α > 0. See [4] pp 78-79 for a sharper bound on the complexity of long integer multiplication.

Proof of Lemma 2. Let

k = ⌈log_{2} \frac{− log_{2} \frac{δ}{2}}{− log_{2} E} ⌉ ,

so that

E^{2^{k}} ≤ \frac{δ}{2} .

Choose ϵ so that Lemma 1 holds, i.e., ϵ < ρ/kA^{k}B, and such that kA^{k}Bϵ < δ/2. This can be done for

ϵ < c_{1}δ(log δ^{−1})^{c_{2}}

for some constants c_{1} and c_{2} dependent of f. The total cost of computing M_{k} is O (a(log_{2} ϵ^{−1})^{1+α}), where a is the arithmetic complexity of the algorithm and α is arbitrarily small. Since f is fixed, a is a constant and k < log_{2} log_{2} δ^{−1} the total cost is

O ((log_{2} δ^{−1})^{1+α})

for all α arbitrarily small.

The complexity bound of Lemma 2 is non-uniform in f. It is also noneffective, in the sense that we give no procedure to estimate k and ϵ without the knowledge of ρ, B and A. Indeed, those quantities may depend on f. It would be of some interest to bound those quantities in a probabilistic setting, similar to Theorem 2.

Definition 4. Let M be an iterative algorithm to compute Φ : F ↦ Φ(F), i.e., lim_{n→∞} M^{n}(F) = Φ(F)). We say that M_{k} is a renormalization of M if

1. M_{k} is a renormalized iterative algorithm to compute φ(f).

2. There are functions ψ and η, defined almost everywhere and computable in finite time, such that the diagram

F ⟶^{ψ} ψ(F)
Φ↓  ↓^{φ}
Φ(F) ⟵^{η} φ(ψ(F))

commutes.

3. There is a function R, computable in finite time, so that the diagrams:

F ⟶^{R^{k}∘ψ}  R^{k}(ψ(F))
M↓  ↓^{M_{k}}
M(F) ⟶^{R^{k+1}∘ψ} R^{k+1}(ψ(M(F)))

commute.

Theorem 1 can be stated now in a more concise way: Let ζ : f ↦ ζ(f) be the function that associates, to any univariate degree d polynomial f, its roots ζ_{1}, …, ζ_{d} ordered by decreasing modulus. We have the following:

Theorem. There is a renormalization of the Graeffe iteration to compute |ζ(f)| = (|ζ_{1}(f)|, …, |ζ_{d}(f)|). Each iteration has arithmetic complexity O(d^{2}) and needs memory O(d).

3 Recurrence Relations and the Renormalization of Graeffe

It is time to construct the renormalized Graeffe iteration. Let f(x) = f_{0} + f_{1}x + ⋯ + f_{d}x^{d}. Let h = Gf be its Graeffe iterate. The coefficients of h can be written as :

h_{i} = (−1)^{d} ∑_{\begin{matrix}0≤i−j≤d \\ 0≤i+j≤d \end{matrix}}(−1)^{i−j}f_{i−j}f_{i+j}  (11)

For convenience, we rewrite (11) as:

h_{i} = (−1)^{d+i}f_{i}^{2} + 2 ∑_{1≤j≤min(i,d−i)}(−1)^{d+i−j}f_{i−j}f_{i+j}  (12)

The next step is to write those equations in terms of the log of the coefficients. More precisely, we will have to deal with the following two quantities:

f_{i}^{log} =^{def} log |f_{i}|  (13)

f_{i}^{arg} =^{def} arg f_{i}  (14)

It is possible now to construct the renormalized Graeffe iteration G_{k}. Recall from diagram (10) that this iteration will map 2^{−k}f^{log} and f^{arg} into 2^{−k−1}h^{log} and h^{arg}.

For that purpose, we introduce the notation:

\{ \begin{matrix}f^{k} & =^{def} & (2^{−k}f^{log}, f^{arg}) \\ h^{k+1} & =^{def} & (2^{−k−1}h^{log}, h^{arg}) \end{matrix}  (15)

We also introduce operators:

(x, α) ⊠ (y, β) =^{def} (x + y, α + β)  (16)
(x, α) ^{λ}  =^{def} (λx, λα)  (17)
z _{k}(x, α) =^{def} (x + 2^{−k} log |z|, α + arg z)  (18)

and

(x, α)⊞_{k}(y, β) =^{def} (2^{−k} log |e^{iα+2^{k}x} + e^{iβ+2^{k}y}| , arg (e^{iα+2^{k}x} + e^{iβ+2^{k}y}))(19)

We remark that the purpose of the sub-index k in the above formulae is to keep track of the degree of the renormalization. For operations which do not change, we omitted the sub-index. The operator z _{k} stands for the multiplication of a renormalized value by a (non-renormalized) constant z

Also, binary renormalized operations are defined for operands with the same renormalization index. Therefore, one should first convert f_{i}^{k} to f_{i}^{k+1} before attempting to ‘multiply’ it with a factor of renormalization index of order k + 1. This conversion will be implicit in the formulae below.

Equation (12) becomes:

h_{i}^{k+1} = ( s_{i0} _{k+1} (f_{i}^{k}) ^{2} ) ⊞_{k+1} ( 2 _{k+1} + _{k+1}_{j=1}^{min(i,d−i)} ( s_{ij} _{k+1}f_{i−j}^{k} ⊠ f_{i+j}^{k})) , (20)

where s_{ij} = (−1)^{d+i−j}. Recall that above, k is a superscript, not an exponent. The ‘renormalized operations’ above are easy to implement in terms of the classical ones. The most delicate being the renormalized sum, so we give here our preferred algorithm Example 2. How to compute the ‘Renormalized sum’:

(c, γ) := (a, α)⊞_{k}(b, β)

    If a > b, do:

s = exp iα + exp (iβ + 2^{k}(b − a))
c = a + 2^{−k} ln |s|
γ = arg(s)

    else

s = exp iβ + exp (iα + 2^{k}(a − b))
c = b + 2^{−k} ln |s|
γ = arg(s)

    endif

We remark that in the above formula, the complex arithmetic operations can be performed in terms of real elementary ones. Moreover, in numerical implementations if k is large enough, as compared to ϵ, it may be faster to approximate (c, γ) with (a, α) or (b, β), whichever is larger.

In order to finish the proof of Theorem 1, we still need a few remarks about the renormalized Graeffe iteration, which we just constructed. Clearly, in order to compute formula (20), we only need memory space of the order O(d) and time of the order O(d^{2}).

When k → ∞, the quantities |f_{i}^{k}| are all convergent.

If G_{k} is the renormalized Graeffe iteration, we define G_{∞} as the limit of G_{k} when k goes to infinity. This means that renormalized sums ⊞_{k} are replaced by their limit ⊞_{∞}, where:

(a, α)⊞_{∞}(b, β) = \{ \begin{matrix}(a, α) & if a > b \\ (b, β) & if a < b \\ undefined & if a = b \end{matrix}

It is easy to see that G_{k} → G_{∞} almost everywhere, pointwise and in the C^{1} topology. In order to avoid a rather tedious calculation, we can establish this fact from the pointwise C^{1} convergence almost everywhere of (a, α)⊞_{k}(b, β) to (a, α)⊞_{∞}(b, β) and similarly for the other renormalized operations.

We define:

ψ(f_{0}, ⋯, f_{d}) = (log |f_{0}|, ⋯, log |f_{d}|; arg |f_{0}|, ⋯, arg |f_{d}|)
η(r_{0}, ⋯, r_{d}; a_{0}, ⋯, a_{d}) = (exp (r_{0} − r_{1}), ⋯, exp (r_{d−1} − r_{d}))
R(r_{0}, ⋯, r_{d}; a_{0}, ⋯, a_{d}) = (\frac{r_{0}}{2} , ⋯, \frac{r_{d}}{2} ; a_{0}, ⋯, a_{d})

We define Φ as the function that associates to a polynomial f the values |ζ_{1}|, …, |ζ_{d}| where ζ_{i} are the roots of f, ordered by decreasing modulus.

Then, we set

φ = η^{−1} ∘ Φ ∘ ψ^{−1}

With the definitions above, G_{k} is indeed a renormalization of the classical Graeffe algorithm to compute Φ:

Proposition 1. Renormalized Graeffe Iteration is a renormalized iterative algorithm (in the sense of Definition 3) to compute Φ.

The proof of this proposition will conclude the proof of Theorem 1. In order to establish Proposition 1, we should check that our algorithm, as described in equation (20), satisfies Axioms 1 to 4.

Axiom 1 is verified by construction. Axiom 2 follows from the recurrence formula (20).

The proof that our algorithm satisfies Axiom 3 will require a technical lemma. Before stating it, we recall some notation:

orb(f) = {f; G_{1}(f); G_{2} ∘ G_{1}(f); …; G_{k} ∘ ⋯ ∘ G_{1}(f); … } is the orbit of f under the sequence (G_{k}). Its closure is denoted by \overline{orb(f)}.

We will show the following lemma, which implies satisfaction of Axiom 3.

Lemma 3. For almost every f and any δ > 0, there exist a compact neighborhood W ⊆ ℝ^{l} × 𝕋^{m} of orb(f) and an integer k_{0} such that, G_{k} is a local diffeomorphism W → G_{k}(W), and such that for k_{0} ≤ k ≤ ∞, the derivative of G_{k} is bounded by 2 + δ

Proof of Lemma 3. The proof is divided in several steps.

Step 1: For any j ∈ ℕ, there is a full-measure set U_{j} such that G_{j} ∘ ⋯ ∘ G_{1} is well-defined, and a local diffeomorphism in U_{j}. Hence, there is a full-measure set U_{j,k} such that G_{k} ∘ (G_{j} ∘ ⋯ ∘ G_{1}) is defined on U_{j,k} and for all f ∈ U_{j,k} there is a neighborhood V_{f} of G_{j} ∘ ⋯ ∘ G_{1}(f) such that G_{k} is a diffeomorphism V_{f} → G_{k}(V_{f}).

Step 2: Let U_{∞} be the set of all f such that G_{k} is a diffeomorphism near (π^{−1} ∘φ)(f) for all values of k that are large enough. Then U_{∞} contains the set of complex polynomials without roots of the same modulus. Hence U_{∞} has full measure.

Step 3: Let U = U_{∞} ∩ (⋂_{j} U_{j}) ∩ (⋂_{jk} U_{jk}). Then U has full measure. Moreover, let f ∈ U. Then G_{k} is a local diffeomorphism with derivative of norm < 2 + δ/2 in an open neighborhood V_{k}(g) of every g ∈ \overline{orbf}, for all k ≥ k_{0}(g). Indeed, if we write G_{k} = G_{∞} + (G_{k} − G_{∞}), we can make the C^{1} norm of the second term arbitrarily small, namely less than δ/2. We know that the norm of the derivative of G_{∞} is precisely 2, hence the bound 2 + δ/2 In the particular case δ = ∞ we can set k_{0} = 1.

Step 4: Since G_{k} → G_{∞} pointwise in the C^{1} topology and for g almost everywhere, we can assume that ⋂_{k≥k_{0}} V_{k}(g) contains an open ball V (g) of center g, where G_{k} is a local diffeomorphism with derivative bounded by 2 + δ/2.

Step 5: Since \overline{orb(f)} is compact, the union ⋃_{g∈\overline{orb(f)}} V (g) has a finite sub-cover ⋃_{g∈Γ} V (g), and we set W = ⋃_{g∈Γ} \overline{V (g)}. Then we set k_{0} = max_{g∈Γ} k_{0}(g), and we obtain that for any k ≥ k_{0}, G_{k} is a local diffeomorphism in W, with derivative of norm bounded by 2 + δ.

The proof that our algorithm satisfies Axiom 4 is divided in two parts , dealing (respectively) with small and large values of k.

Lemma 4. For almost every f and for all k, there exist an open neighborhood U containing f and B > 0 such that for all g ∈ U and for all ϵ small enough

d (G_{k,ϵ}(g) − G_{k}(g)) < Bϵ .

Lemma 5. For almost every f, there exist k_{0}, B > 0, and an open neighborhood U containing f, such that for all g ∈ U and k ≥ k_{0}, for all ϵ small enough,

d (G_{k,ϵ}(g), G_{k}(g)) < Bϵ

Lemmas 4 and 5 together imply that for almost all f, there is a neighborhood U containing f such that for all g ∈ U,

d (G_{k,ϵ}(g), G_{k}(g)) < Bϵ

independently of k.

Since orb(f) is defined for almost every f and admits a compact neighborhood V, we can select a finite subcover of those neighborhoods, and hence find a finite B valid for all g ∈ V. Therefore, our algorithm satisfies Axiom 4.

Proof of Lemma 4: Our iteration G_{k} can be written in terms of the following real operations: +, −, cos, sin, arctan, multiplication by 2^{k}, by 2^{−k}, absolute value, exp, log.

The set of inputs such that a ‘log of zero’ or an ‘absolute value of zero’ occurs has zero measure. Therefore, for almost every f, G_{k}(f) is computed by a composition of analytic functions. Also, for almost every f, none of the output values is zero.

Therefore, for every intermediate quantity x_{l}, the derivative of any output y_{m} with respect to x_{l} is finite (say ≤ D_{lm}). Therefore, a relative perturbation of ϵ in x_{l} leads to a perturbation of size |x_{l}|D_{lm}ϵ in y_{m}. Thus, we set B = \frac{1}{2} ∑_{l,m} D_{lm}|x_{l}|. By continuity, a perturbation of ϵ in each intermediate value x_{l} leads to a perturbation smaller than Bϵ, for input g in a certain neighborhood of f.

Proof of Lemma 5: Let W be the set of all f such that G_{∞}(f) is well-defined. Recall that

\begin{aligned}⊞_{∞} : & (a, α), (b, β) & ↦ & (a, α) & if a > b \\ & & & (b, β) & if b > a \end{aligned}

and that the operator ⊞_{∞} is not defined for a = b. Therefore, W is open and has full measure. Let f ∈ W and U be a small connected neighborhood containing f. Let g ∈ U, then by taking U small enough and k large enough, we can guarantee that G_{k}(g) and G_{k,ϵ}(g) are well-defined. Since U is connected, all the branching outcomes in the computation of G_{k}(g) and G_{k,ϵ} are the same, hence G_{k} restricted to U is a composition of locally analytic functions. We can assume without loss of generality that all derivatives are bounded, hence there is a constant B such that

d (G_{k,ϵ}(g), G_{k}(g)) < Bϵ

for ϵ small enough, but still independent of the choice of g.

This concludes the proof of Proposition 1 and hence of Theorem 1

4 Probability of Success

Through this section, ∥.∥_{d} will denote Weyl’s unitary invariant norm ([53] III–7, pp 137–140) in the space P_{d} of complex polynomials of degree at most d (See [7]). If f(x) = ∑_{i=0}^{d} f_{i}x^{i}, then

∥f∥_{d} = \sqrt{∑_{i=0}^{d} \frac{|f_{i}|}{( \begin{matrix}d \\ i \end{matrix} )}}

This norm is invariant under the following action of the group U(2) of

unitary 2 × 2 matrices: if φ = [ \begin{matrix}α & β \\ γ & δ \end{matrix} ] is unitary, define

f^{φ}(y) = (γy + δ)^{d}f (\frac{αy + β}{γy + δ} )

Thus, ∥f^{φ}∥_{d} = ∥f∥. (See Ch. 12, Th. 1 of [7].)

This is in some sense the most ‘natural’ norm in the space of all polynomials. More information about that norm, its associated probability distribution and its applications can be found in [3, 10, 11, 25, 26, 27, 42, 43, 44, 45, 46].

Let ℙP_{d} be the projectivization of normed complex vector space (P_{d}, ∥.∥_{d}). We can define the ‘sine’ distance in ℙP_{d} by:

d_{ℙ}(f, g) = min_{λ∈ℂ} \frac{∥f − λg∥}{∥f∥} = sin ϱ(f, g)

where ϱ is the usual Riemann metric in ℙP_{d}. Also, there is a natural volume form in ℙP_{d}. Normal invariant distributed polynomials in P_{d} correspond to uniformly distributed ‘polynomials’ in ℙP_{d}.

Let f be a random polynomial (in any of the two equivalent senses above). Then we may consider its roots ζ_{1}, ⋯, ζ_{d} as random variables. The joint distribution of ζ_{1}, ⋯, ζ_{d} was studied by Kostlan in [25]. However, the result below will rely on elementary estimates:

Lemma 6. Let f be random in the sense above. The probability that

min_{|ζ_{i}|>|ζ_{j}|} \frac{|ζ_{i}|}{|ζ_{j}|} > 1 + ϵ

is larger than 1 − Mϵ, where M is a positive constant depending on d.

This result is also true if we chose f random with respect to any other probability distribution with bounded Radon-Nikodym derivative with respect to the volume form in ℙP_{d}. The constant M will have to be multiplied by the maximum of the Radon-Nikodym derivative of the new distribution with respect to the volume form.

Lemma 6 will be a consequence of the ‘Condition Number Theorem’ below. Let

ρ(f) = min_{|ζ_{i}|<|ζ_{j}|} 1 − \frac{|ζ_{i}|}{|ζ_{j}|}

We will interpret ρ(f)^{−1} as a condition number. Let Σ_{G} be the locus of ill-posed problems, i.e., the set of polynomials such that |ζ_{i}| = |ζ_{j}| for some i ≠ j. Then,

Theorem 3 (Condition Number Theorem for Graeffe Iteration).

ρ(f) ≥ \frac{d_{ℙ}(f, Σ_{G})}{\sqrt{d}}

Therefore, the probability that ρ(f) > ϵ is no less than 1 − Vol V_{ϵ\sqrt{d}}Σ_{G} where V_{ϵ\sqrt{d}} denotes an ϵ\sqrt{d}-neighborhood. This volume is of O(ϵ\sqrt{d}) ≤ Mϵ.

Lemma 7. Let f(x) = (x − ζ_{1})g(x) ∈ P_{d}, where g ∈ P_{d−1}. Then

∥g∥_{d−1} ≤ \sqrt{\frac{d}{1 + |ζ_{1}|^{2}}} ∥f∥_{d}

Proof of Lemma 7. We start with the easy case and assume that ζ_{1} = 0. Set g(x) = ∑_{i=0}^{d−1} g_{i}x^{i}. Then, f(x) = ∑_{i=1}^{d} g_{i−1}x^{i}. Now,

∥g∥_{d−1}^{2} = ∑_{i=0}^{d−1} \frac{|g_{i}|^{2}}{( \begin{matrix}d − 1 \\ i \end{matrix} )} = ∑_{i=1}^{d} \frac{|g_{i−1}|^{2}}{( \begin{matrix}d − 1 \\ i − 1 \end{matrix} )} = ∑_{i=1}^{d} \frac{d}{i} \frac{|g_{i−1}|^{2}}{( \begin{matrix}d \\ i \end{matrix} )} ≤ d∥f∥_{d}^{2}

and hence ∥g∥_{d−1} ≤ \sqrt{d}∥f∥_{d}.

For the general case, we will use U(2)-invariance of ∥.∥_{d} and ∥.∥_{d−1}. Let φ be a convenient unitary matrix:

φ = \frac{1}{\sqrt{1 + |ζ_{1}|^{2}}} [ \begin{matrix}1 & ζ_{1} \\ −ζ̄_{1} & 1 \end{matrix} ]

Set f^{φ} = f ∘ φ and g^{φ} = g ∘ φ. The choice of φ has the particularity that f^{φ}(0) = f(ζ_{1}) = 0. We can compute f^{φ} in terms of g^{φ}:

f^{φ}(y) = \sqrt{1 + |ζ_{1}|^{2}} y g^{φ}(y)

Using the easy case,

∥\sqrt{1 + |ζ_{1}|^{2}}g^{φ}∥_{d−1} ≤ \sqrt{d}∥f^{φ}∥_{d}

By U(2)-invariance,

∥g∥_{d−1} = ∥g^{φ}∥_{d−1}
= \frac{1}{\sqrt{1 + |ζ_{1}|^{2}}} ∥\sqrt{1 + |ζ_{1}|^{2}}g^{φ}∥_{d−1}

≤ \begin{matrix}√ & \overline{d} \\ √ & \overline{1 + |ζ_{1}|^{2}} \end{matrix} ∥f^{φ}∥_{d}

= \begin{matrix}√ & \overline{d} \\ √ & \overline{1 + |ζ_{1}|^{2}} \end{matrix} ∥f∥_{d}

Lemma 8. Let g ∈ P_{d−1}. Then ∥g∥_{d} ≤ ∥g∥_{d−1}.

Proof of Lemma 8.

∥g∥_{d}^{2} = ∑_{i=0}^{d−1} \frac{|g_{i}|^{2}}{( \begin{matrix}d \\ i \end{matrix} )} = ∑_{i=0}^{d−1} \frac{d − i}{d} \frac{|g_{i}|^{2}}{( \begin{matrix}d − 1 \\ i \end{matrix} )} ≤ ∥g∥_{d−1}^{2}

Putting Lemma 7 and Lemma 8 together, Lemma 9. Let f(x) = (x − ζ_{1})g(x) ∈ P_{d}, where g ∈ P_{d−1}. Then

∥g∥_{d} ≤ \begin{matrix}√ & \overline{d} \\ √ & \overline{1 + |ζ_{1}|^{2}} \end{matrix} ∥f∥_{d}

Proof of Theorem 3. Let f(x) = (x − ζ_{1})(x − ζ_{2}) ⋯ (x − ζ_{d}) and order the ζ_{i}’s such that

ρ(f) = min_{|ζ_{i}|<|ζ_{j}|} 1 − \frac{|ζ_{i}|}{|ζ_{j}|} = 1 − \frac{|ζ_{2}|}{|ζ_{1}|}

Define h(z) = (x − ζ_{1}(1 − ρ))(x − ζ_{2}) ⋯ (x − ζ_{d}) ∈ Σ_{G}. Then

∥f − h∥_{d} = ∥ζ_{1}ρ(x − ζ_{2}) ⋯ (x − ζ_{d})∥_{d}
= |ζ_{1}|ρ∥(x − ζ_{2}) ⋯ (x − ζ_{d})∥_{d}
≤ |ζ_{1}|ρ\begin{matrix}√ & \overline{d} \\ √ & \overline{1 + |ζ_{1}|^{2}} \end{matrix} ∥f∥_{d}

Hence,

\frac{∥f − h∥_{d}}{∥f∥_{d}} ≤ \frac{|ζ_{1}|}{\sqrt{1 + |ζ_{1}|^{2}}} ρ\sqrt{d} ≤ ρ\sqrt{d}

Hence,

d_{ℙ}(f, Σ_{G}) ≤ ρ\sqrt{d}

Proof of Theorem 2. We set δ = Mϵ, where M is the constant such that the volume of an ϵ\sqrt{d} neighborhood of Σ_{G} is less than Mϵ. With probability larger than 1 − δ,

max_{|ζ_{i}|<|ζ_{j}|} \frac{|ζ_{i}|}{|ζ_{j}|} > 1 − ϵ .

We now use the fact that, for any N > (log 2)/ϵ, we have (1 − ϵ)^{N} < 1/2. We set k_{1} = 1 + ⌈log_{2} ϵ^{−1}⌉ and using N = 2^{k_{1}} in the previous formula, we obtain:

max_{|ζ_{i}|<|ζ_{j}|} \frac{|ζ_{i}|^{2^{k_{1}}}}{|ζ_{j}|^{2^{k_{1}}}} < (1 − ϵ)^{2^{k_{1}}} < \frac{1}{2} .

An extra 1 + log_{2} b iterations ensures that

max_{|ζ_{i}|<|ζ_{j}|} \frac{|ζ_{i}|^{2^{k}}}{|ζ_{j}|^{2^{k}}} < 2^{−1−b} ,

provided that

k = k_{1} + 1 + log_{2} b .

We claim that we have in this case that the sum

∑_{i_{1}<···<i_{r}} ζ_{i_{1}}^{2^{k}} … ζ_{i_{r}}^{2^{k}} = ζ_{1}^{2^{k}} … ζ_{r}^{2^{k}}(1 + ε) ,

where

|ε| < 2^{−b} .

Indeed, assume we have reordered the roots so that |ζ_{1}| > ⋯ > |ζ_{d}|. It then follows that

|ζ_{i_{1}}^{2^{k}} … ζ_{i_{r}}^{2^{k}}| < |ζ_{1}^{2^{k}} … ζ_{r}^{2^{k}}|2^{−(1+b)(i_{1}+···+i_{r}−r(r+1)/2)} .

Now, for r fixed, let

S_{k} = {(i_{1}, …, i_{r}) | i_{1} < ⋯ < i_{r} and i_{1} + … i_{r} − r(r + 1)/2 = k} .

Note that S_{0} = {(1, …, r)}, and S_{1} = {(1, …, r − 1, r + 1)}. Also, S_{2} = {(1, …, r − 1, r + 2); (1, …, r − 2, r, r + 1)}

Figure 1: Examples of Newton diagrams

In general, every multi-index i_{1} < ⋯ < i_{r} may be obtained by starting from 1 < 2 < ⋯ < r and increasing one of the indices, in such a way not two indices are equal. Then S_{k} is the set of multi-indices obtained after k steps. Therefore, we may bound #S_{k+1} ≤ d #S_{k} to get

∑_{k≥1} 2^{−(b+1)k} #S_{k} ≤ 2^{−b−1}(1 + 2^{−b} + 2^{−b}d + 2^{−2b−1}d^{2} …) ≤ 2^{−b} ,

under the assumption that d ≤ 2^{b}. This concludes the proof of Theorem 2.

5 Newton Diagram Revisited

Following Ostrowskii [34], we introduce Newton diagram, but now with a log scale. The Newton diagram of a polynomial f is the graph of the piecewise linear function defined by g(i) = − log |f_{i}|. (Our construction differs from Ostrowskii’s point of view. He used to define the Newton diagram as the convex hull of that set of points). It makes sense to renormalize Newton diagram, as we did with polynomials G^{k}f.

Therefore, we may speak of a Newton diagram on the limit. Assume (as we did through this point) that all roots of f have different moduli. Then the derivative of this diagram on the piecewise linear part over [i, i + 1] will correspond to log |ζ_{d−i}|, where ζ_{i} is the i-th largest root of f in modulus. It follows that the Newton diagram will be the graph of a convex function.

We should consider now the general case. If f has several roots of the same moduli, some of the ratios f_{i}/f_{i+1} will converge to the ratios of the coefficients of the factor of f containing those roots. The case of three roots of same moduli is illustrative:

Let ζ_{1}, …, ζ_{d} be arranged by non-increasing moduli, and assume |ζ_{i}| = |ζ_{i+1}| = |ζ_{i+2}|. Then,

f_{d−i} = ζ_{1} … ζ_{i−1} (ζ_{i} + ζ_{i+1} + ζ_{i+2}) + …  (21)
f_{d−i−1} = ζ_{1} … ζ_{i−1} (ζ_{i}ζ_{i+1} + +ζ_{i}ζ_{i+2} + ζ_{i+1}ζ_{i+2}) + …  (22)
f_{d−i−2} = ζ_{1} … ζ_{i+2} + …  (23)

Hence, on the limit, we approximate:

\frac{f_{d−i−2}}{f_{d−i−1}} ≃ \frac{ζ_{i}ζ_{i+1}ζ_{i+2}}{ζ_{i}ζ_{i+1} + +ζ_{i}ζ_{i+2} + ζ_{i+1}ζ_{i+2}}  (24)
\frac{f_{d−i−1}}{f_{d−i}} ≃ \frac{ζ_{i}ζ_{i+1} + +ζ_{i}ζ_{i+2} + ζ_{i+1}ζ_{i+2}}{ζ_{i} + ζ_{i+1} + ζ_{i+2}}  (25)
\frac{f_{d−i}}{f_{d−i+1}} ≃ ζ_{i} + ζ_{i+1} + ζ_{i+2}  (26)

The moduli of those roots is therefore given by:

\sqrt[3]{|\frac{f_{d−i−2}}{f_{d−i−1}} | |\frac{f_{d−i−1}}{f_{d−i}} | |\frac{f_{d−i}}{f_{d−i+1}} |}  (27)

The same formula extends for factors of any degree, provided they have roots in a circle and all other roots are far away from this circle.

It is useful to have a decision criterion for the existence of factors of degree greater than 1. We will do that for degree-2 factors, since this is the interesting case for real polynomials. See Ostrowskii [34] for more results.

Clearly, it is enough to consider the case of a polynomial f of degree 2:

f(x) = f_{2}x^{2} + f_{1}x + f_{0} = x^{2} + (−ζ_{1} − ζ_{2})x + ζ_{1}ζ_{2}  (28)

In case |ζ_{1}| ≫ |ζ_{2}|, we have:

− log (\frac{|f_{2}|}{|f_{1}|}) + log (\frac{|f_{1}|}{|f_{0}|}) ≥ log \frac{|ζ_{1}|}{|ζ_{2}|} ≫ 0  (29)

In case R = |ζ_{1}| = |ζ_{2}|, we can bound:

− log (\frac{|f_{2}|}{|f_{1}|}) + log (\frac{|f_{1}|}{|f_{0}|}) < log 4  (30)

However, we should look at the renormalized Newton diagram. In that case, we should divide equations (29) and (30) by 2^{k} , and expect that if |ζ_{1}| ≠ |ζ_{2}| then the following holds:

−2^{−k} log (\frac{|f_{2}|}{|f_{1}|}) + 2^{−k} log (\frac{|f_{1}|}{|f_{0}|}) ≥ σ  (31)

where σ is an a priori bound on root moduli separation. It may be obtained from a probabilistic analysis (Lemma 6), or from any other a priori knowledge on the polynomial; for the choice of σ, notice that if |ζ_{1}| = |ζ_{2}|,

−2^{−k} log (\frac{|f_{2}|}{|f_{1}|}) + 2^{−k} log (\frac{|f_{1}|}{|f_{0}|}) < 2^{−k} log 4  (32)

Thus we may choose k such that 2^{−k} log 4 < σ. In case equation (31) is not satisfied, it is reasonable to assume that the two roots have the same modulus indeed.

6 Numerical Results

The results discussed below are tentative, and our algorithm deserves further experimentation. Moreover, at this moment we are using a tentative algorithm to find the argument of the roots. Those matters will be dealt with in a subsequent paper ([28]). The purpose of this section is to illustrate the numerical properties of Renormalized Graeffe Iteration, when applied to random real polynomials.

Figure 2: Zeros of a random polynomial

A C implementation of our algorithm was tested for pseudo-random real and complex polynomials. The arguments of the solutions were recovered using an algorithm derived from the Renormalized Graeffe Iteration.

IMPLEMENTATION DETAILS: Polynomials with zero coefficients do arise in practice, and we explain now how to deal with the ‘log of zero’ problem. The IEEE floating point arithmetic, implemented by most modern computers, has a few useful features for this sort of situation. When one asks a computer to produce log 0, it returns a special IEEE value called −∞. This is not an error, and further calculations can be carried out (as long as they are defined). If they are not defined (e.g. ∞ − ∞) then another special value, called a ‘not a number’ is returned.

This last case (∞ − ∞) can appear during the computation of a renormalized sum. It can be dealt by testing a and b in Example 2 for finitesness. In case a or b is infinite, then c should be set to ∞. For an introduction of IEEE arithmetic, see [19] pages 45-48 or [12] pages 9-15. The correctness of the results was certified by estimates as in [26].

Experiments were performed in a Pentium 66 system running Linux operating system. Since the objective here was to illustrate the asymptotic behaviour of the algorithm, we did not perform experiments in other systems. Those would be necessary if one wanted to compare with other algorithms with same asymptotic properties. However, this goes far beyond the scope of this paper.

The table below shows the average and median user time in a Pentiumbased computer, using ‘long double’ precision. Time does not include validation time. Ten pseudo-random polynomials were tested for each degree. The actual experimental data is plotted in Figure 6.

Degree  Real polynomials  Complex polynomials
avg time (s) median time(s) avg time (s) median time (s)
100  0.87  0.87  1.19  1.18
200  3.23  3.24  4.35  4.30
300  7.07  7.07  9.99  9.73
400  12.60  12.58  17.28  17.10
500  19.41  19.38  26.75  26.71
600  27.67  27.73  37.02  35.96
700  37.50  37.32  51.30  49.91
800  48.89  48.72  65.39  63.61
900  61.56  61.06  79.28  78.74
1000  75.89  75.47  102.33  101.80

We also computed (approximately) the relative separation of the moduli of the solutions ζ_{i}:

min_{|ζ_{i}|<|ζ_{j}|} \frac{|ζ_{i}| − |ζ_{j}|}{\sqrt{1 + |ζ_{i}|^{2}}}

The values obtained are also plotted in figures 3 to 6.

Further experimentation is necessary to obtain data about polynomials of degree ≫ 1000. Indeed, due to underflow, we cannot represent random high-degree polynomials in the usual floating point representation.

7 Acknowledgments

JPZ was supported by CNPq and GM was partially supported by CNPq. The authors wish to thank the following people for their help: Jean-Pierre Dedieu, Carlos Isnard, Peter Kirrinnis, Welington de Melo, Victor Pan, Mike Shub, Jose Felipe da Silva, Steve Smale, Benar Fux Svaiter, Jean-Claude Yakoubsohn.

Figure 3: Timing for 100 random real polynomials 33 Figure 4: Timing for 100 random complex polynomials

34 Figure 5: Separation for 100 random real polynomials

35 Figure 6: Separation for 100 random complex polynomials

36

The authors gratefully acknowledge a large number of helpful comments by the anonymous referees, which helped to enhance the presentation tremendously.

References

[1] E. H. Bareiss. Resultant Procedure and the Mechanization of the Graeffe Process. J. ACM 7 pp 346-386, 1960.

[2] E. H. Bareiss. The numerical Solution of Polynomial Equations and the Resultant Proceudures. In: A. Ralston and H. S. Wilf, eds, Mathematical Methods for Digital Computers Vol II, pp 185-214 John Wiley and Sons, Inc, New York, 1967.

[3] B. Beauzamy and J. Degot: Differential identities. Transactions of the AMS 347, pp 2607-2619, 1995.

[4] Dario Bini and Victor Pan. Polynomial and Matrix Computations, Vol 1 Progress in Theoretical Computer Science, Birkhäuser, Boston, 1994.

[5] Dario Bini and Victor Pan, Graeffe’s, Chebyshev-like and Cardinal’s Processes for Splitting a Polynomial into Factors. Journal of Complexity 12, pp 492-511, 1996.

[6] Sharon L. Blish and James H. Curry, On the Geometry of Factorization Algorithms. in: Eugene L. Allgower and Kurt Georg, eds, Computational Solution of Nonlinear Systems of Equations. Lectures in Applied Mathematics 26, American Mathematical Society, Rhode Island, 1990.

[7] Lenore Blum, Felipe Cucker, Mike Shub and Steve Smale. Complexity and Real Computation Springer-Verlag, New York, 1998.

[8] Jean-Paul Cardinal, On Two Iterative Methods for Approximating the Roots of a Polynomial. The Mathematics of Numerical Analysis. 1995 AMS-SIAM Summer Seminar in Applied Mathematics, July 17- August 11, 1995, Park City, UT. Lectures in Applied Mathematics 32, pp 165–188. American Mathematical Society, Rhode Island, 1996.

[9] Jean-Pierre Dedieu, A propos de la méthode de Graeffe. C.R. Acad. Sci. Paris, t. 309, pp 1019-1022, 1989.

[10] Jean-Pierre Dedieu, Approximate Solutions of Numerical Problems, Condition Number Analysis and Condition Number Theorem. in: James Renegar, Michael Shub, Steve Smale (eds), The Mathematics of Numerical Analysis. 1995 AMS-SIAM Summer Seminar in Applied Mathematics, July 17- August 11, 1995, Park City, UT. Lectures in Applied Mathematics 32, pp 263-283. American Mathematical Society, Rhode Island, 1996.

[11] Jean-Pierre Dedieu, Condition Number Analysis for Sparse Polynomial Systems. in: Felipe Cucker and Mike Shub (eds), Foundations of Computational Mathematics – Selected Papers of a Conference held at IMPA in Rio de Janeiro, Jan 1997, pp 75–101. Springer-Verlag, Berlin, 1997.

[12] James W. Demmel, Applied Numerical Linear Algebra, SIAM, Philadelphia, 1997.

[13] Alan Edelman and H. Murakami, Polynomial Roots from Companion Matrix Eigenvalues. Mathematics of Computation 64 No 210, pp 763-776, 1995.

[14] Ioannis Z. Emiris, André Galligo, and Henri Lombardi, Numerical Univariate Polynomial GCD. The Mathematics of Numerical Analysis. 1995 AMS-SIAM Summer Seminar in Applied Mathematics, July 17- August 11, 1995, Park City, UT. Lectures in Applied Mathematics 32, pp 323-343. American Mathematical Society, Rhode Island, 1996.

[15] Mitchell J. Feigenbaum Quantitative Universality for a class of Nonlinear Transformations. J. of Statistical Physics 19, no.1, 1978.

[16] S. Goedecker, Remark on Algorithms to Find Roots of Polynomials. SIAM J. Sci. Computing 15 No 5, pp 1059-1063, 1994.

[17] A. A. Grau. On the Reduction of Number Range in the Use of the Graeffe Process. Journal of the ACM, 1963.

[18] P. Henrici. Applied and Computational Complex Analysis. Wiley, New York, 1974.

[19] Nicolas J. Higham, Accuracy and Stability of Numerical Algorithms, SIAM, Philadelphia, 1996.

[20] Alston S. Householder. Dandelin, Lobatchevskii or Graeffe ? American Mathematical Monthly.

[21] Anthony Hyman, Charles Babbage, Pioneer of the Computer. Princeton University Press, 1982.

[22] M.A.Jenkins and J.F.Traub. A Three-Stage variable shift iteration for polynomial zeros and its relation to generalized Rayleigh iteration. Numer. Math. 14, 1970.

[23] Peter Kirrinnis, Partial Fraction Decomposition in ℂ^{n} and simultaneous Newton Iteration for Factorization in ℂ[z]. Preprint, Bonn, 1995.

[24] Peter Kirrinnis, Newton Iteration Towards a Cluster of Polynomial Zeros. in: Felipe Cucker and Mike Shub (eds), Foundations of Computational Mathematics – Selected Papers of a Conference held at IMPA in Rio de Janeiro, Jan 1997, pp 193–215. Springer-Verlag, Berlin, 1997.

[25] Eric Kostlan Random Polynomials and the Statistical Fundamental Theorem of Algebra. Preprint, MSRI, 1987.

[26] Gregorio Malajovich, On generalized Newton algorithms: quadratic convergence, path-following and error analysis. Theoretical Computer Science 133, 65-84 (1994).

[27] Gregorio Malajovich and Jorge P. Zubelli. A Fast and Stable Algorithm for Splitting Polynomials. Computers and Mathematics with Applications, 33:3, 1–23, 1997.

[28] Gregorio Malajovich and Jorge P. Zubelli. Tangent Graeffe Iteration. Informes de Matemática. Serie B. Número 119, IMPA, Junho 1998 pp 1–34.

[29] Curtis T. McMullen. Complex Dynamics and Renormalization Annals of Mathematics Studies 135. Princeton University Press, Princeton, New Jersey, 1994.

[30] Welington de Melo and S. van Strien, One-Dimensional Dynamics. Springer Verlag, 1993.

[31] Bernard Mourrain and Victor Y. Pan, Solving Special Polynomial Systems by Using Structured Matrices and Algebraic Residues. in: Felipe Cucker and Mike Shub (eds), Foundations of Computational Mathematics – Selected Papers of a Conference held at IMPA in Rio de Janeiro, Jan 1997, pp 287– 304. Springer-Verlag, Berlin, 1997.

[32] C. Andrew Neff and John H. Reif. An Efficient Algorithm for the Complex Roots Problem. Journal of Complexity 12, 81-115, 1996.

[33] Erich Novak and Henryk Woźniakowski, Topological Complexity of Zero-Finding. Journal of Complexity 12, pp380-400, 1996.

[34] Alexandre Ostrowskii. Recherches sur la Méthode de Graeffe et les Zéros de Polynomes et des Séries de Laurent. Acta Mathematica 72, 99-257, 1940.

[35] Jacob Palis and Floris Takens. Hyperbolicity and Sensitive Chaotic Dynamics at Homoclinic Bifurcations Cambridge Studies in Advanced Mathematics 35. CUP, Cambridge, 1993.

[36] Victor Y.Pan, Deterministic Improvement of Complex Polynomial Factorization Based on the Properties of the Associated Resultant. Computers Math. Applic. 30 No 2, pp 71-94,1995.

[37] Victor Y.Pan, Optimal and Nearly Optimal Algorithms for Approximating Polynomial Zeros. Computers Math. Applic. 31 No 12, pp 97–138, 1996.

[38] Victor Y.Pan. Solving a Polynomial Equation: Some History and Recent Progress. SIAM Review 39 No 1, 1997.

[39] Victor Y.Pan, Myongh-hi Kim, Akimou Sadikou, Xiaohan Huang and Aliong Zheng, On Isolation of Real and nearly Real Zeros of a Univariate Polynomial and its Splitting into Factors. Journal of Complexity 12, pp 772–594, 1996.

[40] M.G. Pólya, Sur la méthode de Graeffe. Comptes Rendus de l’Académie des Sciences de Paris, 14 avril 1913.

[41] James Renegar, On the Cost of Approximating all Roots of a Complex Polynomial Mathematical Programming 32 (1985), 319-336.

[42] Michael Shub and Steve Smale, On the Complexity of Bezout’s Theorem I - Geometric aspects. Journal of the AMS, 6(2), (1993).

[43] Michael Shub and Steve Smale, On the complexity of Bezout’s Theorem II - Volumes and Probabilities. in: F. Eysette and A. Galligo, eds: Computational Algebraic Geometry. Progress in Mathematics 109 267-285, Birkhauser (1993).

[44] Michael Shub and Steve Smale, Complexity of Bezout’s Theorem III ; Condition number and packing. Journal of Complexity 9, 4-14 (1993).

[45] Michael Shub and Steve Smale, Complexity of Bezout’s Theorem IV ; Probability of success ; Extensions. Siam J. Numerical Analysis 33 No 1, pp 128-148, 1996.

[46] Michael Shub and Steve Smale, Complexity of Bezout’s Theorem V: Polynomial time. Theoretical Computer Science 133(1), 141-164 (1994)

[47] Steve Smale, Complexity Theory and Numerical Analysis. Preprint, Hong Kong, 1996. To appear, Acta Numerica.

[48] Arnold Schönhage, The Fundamental Theorem of Algebra in Terms of Computational Complexity – Preliminary Report Preprint, Tübingen, 1982.

[49] Arnold Schönhage, Equation solving in Terms of Computational Complexity. Proceedings of the International Congress of Mathematicians, Berkeley, 1986. 131-153, American Mathematical Society, 1987

[50] Kim-Chuan Toh and Lloyd N. Trefethen Pseudozeros of Polynomials and Pseudospectra of companion matrices. Numer. Math. 68, pp 403–425 (1994).

[51] J.V.Uspensky. Theory of equations. McGraw Hill, New York, 1948.

[52] Victor Vassiliev Topological Complexity of Root-Finding Algorithms. The Mathematics of Numerical Analysis. 1995 AMS-SIAM Summer Seminar in Applied Mathematics, July 17- August 11, 1995, Park City, UT. Lectures in Applied Mathematics 32, pp 831–856. American Mathematical Society, Rhode Island, 1996.

[53] Hermann Weyl. The Theory of Groups and Quantum Mechanics Dover Publications, Inc. c.1931, 1951.