Tangent 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

Revised version, August 26, 1999

Abstract

Graeffe iteration was the choice algorithm for solving univariate polynomials in the XIX-th and early XX-th century. In this paper, a new variation of Graeffe iteration is given, suitable to IEEE floating-point arithmetics of modern digital computers.

We prove that under a certain generic assumption the proposed algorithm converges. We also estimate the error after N iterations and the running cost.

The main ideas from which this algorithm is built are: classical Graeffe iteration and Newton Diagrams, changes of scale (renormalization), and replacement of a difference technique by a differentiation one.

The algorithm was implemented successfully and a number of numerical experiments are displayed.

Contents

1 Introduction 2 1.1 Main Result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 What the Graeffe Iteration is; Its historical weaknesses . . . . . . . . 4

2 Renormalizing Graeffe 6 2.1 The Renormalized Graeffe Iteration . . . . . . . . . . . . . . . . . . 6 2.2 The Renormalized Newton Diagram . . . . . . . . . . . . . . . . . . 7

3 Computation of the Newton Diagram 12 3.1 Algorithm and Main Statements . . . . . . . . . . . . . . . . . . . . 12 3.2 Some Estimates about Symmetric Functions . . . . . . . . . . . . . . 13 3.3 A Decision Criterion . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.4 Proof of Proposition 2 . . . . . . . . . . . . . . . . . . . . . . . . . 18

4 Tangent Graeffe Iteration 19 4.1 Perturbation Methods, Infinitesimals, 1-Jets of Polynomials . . . . . . 19 4.2 The Iteration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 4.3 Convergence Results . . . . . . . . . . . . . . . . . . . . . . . . . . 20 4.4 The Main Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . 23 4.5 Proof of Lemmas 5 and 6 . . . . . . . . . . . . . . . . . . . . . . . 24 4.6 ‘Deterioration of Condition’ and Stability Properties . . . . . . . . . 27

5 Numerical Results and Final Remarks 27 5.1 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 5.2 Further Practical Remarks . . . . . . . . . . . . . . . . . . . . . . . 28

1 Introduction

Many present day numerical algorithms have originated in highly acclaimed methods dating from last century or even earlier. Such is certainly the case of Euler’s method or of Newton’s method, whose numerical and theoretical consequences still impact us today [15, 33, 34, 35, 36, 37, 39].

Graeffe’s classical method for finding simultaneously all roots of a polynomial was introduced independently by Graeffe, Dandelin and Lobatchevsky [11]. Its simplicity, as well as importance throughout last century indicate its potential as an effective numerical algorithm.

Surprisingly, Graeffe’s method has not received much attention in present day numerical computations. Very few modern discussions about it or its applications can be found. See the review by V. Pan [28], and also [2, 5, 6, 8, 16, 21, 22, 24, 27, 29, 32].

One of the main reasons for Graeffe’s lack of popularity stems from the fact that its traditional form leads to exponents that easily exceed the maximum allowed by floating-point arithmetic. Other reasons, such as the “chaotic” behavior of the arguments of the roots of the iterates contribute to such stigma.

Also, Graeffe iteration is a many-to-one map. It can map well-conditioned polynomials into ill-conditioned ones, as pointed out by Wilkinson in [40]. We shall refer to this as ‘Wilkinson’s Deterioration of Condition’.

In this work we present a version of Graeffe’s algorithm, which is well suited for floating-point arithmetic computations. Furthermore, it has excellent complexity and memory allocation characteristics. Our method computes both the moduli and the argument of all the roots, provided that certain generic conditions are satisfied. These claims are backed by our theoretical results presented in the next section, and proved throughout the paper, as well as the numerical experiments presented in the end of the paper.

The main ingredients in our approach are the following:

• The idea of renormalizing the relevant operations at each iteration step, akin to what is done in dynamical systems and physics [19, 23].

• The idea of using the differential of our Renormalized Graeffe iteration as a way of keeping the information concerning the argument of the roots. This will allow us to avoid the harmful effects of Wilkinson’s ‘deterioration of condition’, as discussed in Section 4.6.

• A renormalized version of Newton’s diagram that allows us to recognize and locate pairs of conjugate roots, as well as roots of higher multiplicity.

The first idea mentioned above was developed in our earlier work [22], which in a certain sense laid the conceptual framework for our present approach. It is not however essential in understanding the proofs presented herein.

The second ingredient mentioned above is explained and motivated in Section 4. In rather vague terms it could be compared to the advantage of using derivatives, when those are available, as compared to using differences. This idea can be traced back to Brodetsky and Smeal [4] in 1924, in a more ad-hoc fashion. We are not aware of recent applications of that method in modern literature.

Finally, the concept of Newton’s diagram, as well as the power of Graeffe’s method was present throughout Ostrowski’s masterpiece [25]. While writing the present paper we could not help but wonder what would have been the outcome of that research if he had available at that point the present day technology of high speed computers.

We wish to thank two anonymous referees for their comments and for suggesting some extra references such as [4], which we were not aware of in the first draft.

1.1 Main Result

We will introduce an algorithm for solving real and complex univariate polynomials. The following genericity condition will be required at input:

Definition 1. A real polynomial f will be called circle free if, and only if, for any couple ζ, ξ of distinct roots of f, one has either |ζ| ≠ |ξ|, or ζ = ξ̄.

Definition 2. A complex polynomial f will be called circle free if, and only if, for any couple ζ, ξ of distinct roots of f one has |ζ| ≠ |ξ|.

It is obvious that given any real polynomial f, one can obtain a circle free polynomial by (pre)composing it with a conformal transform of the form:

φ : ℂ̄ → ℂ̄
x ↦ φ(x) = \frac{x cos θ−sin θ}{x sin θ+cos θ}

and then clearing denominators; for all but a finite number of θ ∈ (−π, π], the resulting polynomial

f̃(x) = (x sin θ + cos θ)^{d}f (\frac{x cos θ − sin θ}{x sin θ + cos θ} )

is circle free, where d is the degree of f.

Tangent Graeffe Iteration will be shown to converge for all circle free polynomials; Given an arbitrary polynomial, one can first find all zero roots (in the obvious way), apply a random conformal transform, then Tangent Graeffe Iterations, and finally recover the roots of the original polynomial.

When counted with multiplicity, the roots of a circle free polynomial can be canonically ordered by:

1. |ζ_{1}| ≤ ⋯ ≤ |ζ_{d}|

2. In the real case, 2.1. If |ζ_{i}| = |ζ_{i+1}| then ζ_{i} = \overline{ζ_{i+1}}. 2.2. If i = 1 or |ζ_{i−1}| < |ζ_{i}| then Im ζ_{i} ≥ 0.

If we assume that all the arithmetical operations are performed exactly (including transcendental), the mathematical properties of the algorithm can be summarized by:

Theorem 1. Let f be a real (resp. complex) circle free degree d polynomial, not vanishing at 0. Denote by ζ the vector of all the roots of f with multiplicity canonically ordered as above.

Then, a total of N iterations of Renormalized Tangent Graeffe (Algorithm 6) produces ζ^{(N)} ∈ ℂ^{d}, such that

ζ^{(N)} ⟶ ζ (as N → ∞).

The running time for each iteration is O(d^{2}) exact arithmetic operations (including transcendental operations). The relative truncation error bound in each coordinate after N iterations is 2^{−2^{N−C}}, where C depends on f.

1.2 What the Graeffe Iteration is; Its historical weaknesses

In this section we shall briefly review the main ideas behind the method and describe also some of its weaknesses.

Graeffe iteration maps a degree d polynomial f(x) into the degree d polynomial

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

If ζ_{1}, ζ_{2}, … ζ_{d} are the roots of f, then the roots of Gf are ζ_{1}^{2}, ζ_{2}^{2}, … ζ_{d}^{2}

Assume that g = G^{N}f is the N-th iterate of f. Then, assuming that f is monic, the coefficients of g(x) = g_{0} + g_{1}x + ⋯ + g_{d}x^{d} satisfy:

g_{0} = (−1)^{d} σ_{d} (ζ_{1}^{2^{N}} , ζ_{2}^{2^{N}} , ⋯, ζ_{d}^{2^{N}} )
g_{1} = (−1)^{d−1} σ_{d−1} (ζ_{1}^{2^{N}} , ζ_{2}^{2^{N}} , ⋯, ζ_{d}^{2^{N}} )
⋮
g_{j} = (−1)^{d−j} σ_{d−j} (ζ_{1}^{2^{N}} , ζ_{2}^{2^{N}} , ⋯, ζ_{d}^{2^{N}} )
⋮
g_{d} = σ_{0} = 1

where σ_{j} is the j-th elementary symmetric function. In the particular case that |ζ_{1}| < |ζ_{2}| < ⋯ < |ζ_{d}|, we can further approximate

g_{0} = (−1)^{d}ζ_{1}^{2^{N}} ζ_{2}^{2^{N}} … ζ_{d}^{2^{N}}
g_{1} ≃ (−1)^{d−1}ζ_{2}^{2^{N}} … ζ_{d}^{2^{N}}
⋮
g_{j} ≃ (−1)^{d−j}ζ_{j+1}^{2^{N}} … ζ_{d}^{2^{N}}
⋮
g_{d} = σ_{0} = 1

Hence, it is possible to determine

ζ_{j}^{2^{N}} ≃ −\frac{g_{j}}{g_{j+1}} .

We stress two main weaknesses. The first big weakness of classical Graeffe iteration is coefficient growth. As the coefficients of g_{j} grow doubly exponentially in the number of iterations, the exponent (not the mantissa) of the floating-point system gets overflowed:

Example 1. Let f have roots 1, 2, 3, 4. Then the N-th Graeffe iterate of f has roots 1, 2^{2^{N}}, 3^{2^{N}}, 4^{2^{N}}. The coefficient g_{0} is 24^{2^{N}}. If N = 8, then g_{0} is approximately 1.68 × 2^{1173}, while IEEE double precision numbers used in most modern computers cannot contain floating point values more than 2^{1024}, since the exponent is represented by 11 bits (sign included) [10]. (As a matter of fact, the representation is a little more complicated, as it allows for ‘subnormal’ numbers [7]). Therefore, we would have an overflow when computing the 8-th Graeffe iterate of f.

Example 2. On the example above, assume that f would have an additional root 1.01. Namely,

\begin{aligned}f(x) & = & (x − 1)(x − 1.01)(x − 2)(x − 3)(x − 4) \\ & = & −24.24 + 74.5x − 85.35x^{2} + 45.1x^{3} − 11.01x^{4} + x^{5} \end{aligned}

We will show that 8 Graeffe iterations are not enough to compute the first root (namely 1) to an accuracy of 10^{−4}. However, as shown in Example 1, 8 iterates are enough to overflow the IEEE double precision number system.

Indeed, the first root is obtained as:

ζ^{2^{N}} ≃ −\frac{g_{0}}{g_{1}} = −\frac{1.01^{2^{N}} 24^{2^{N}}}{(1^{2^{N}} +1.01^{2^{N}})24^{2^{N}}}
= \frac{1}{1+1.01^{−2^{N}}}
≃ 0.927 ( for N = 8).

Thus,

ζ ≃ 1 − 2.9 × 10^{−4}

The error obtained is therefore larger than 10^{−4}.

The introduction of the idea of renormalization allows us to avoid coefficient growth, and replace a diverging algorithm by a convergent one (See Section 2 and also [22]). Alternative approaches for the number range growth are suggested in [9] and in [8].

A certain geometrical invariant of the polynomial, the limiting Newton diagram, appears naturally in the context of Renormalized Graeffe Iteration. It allows to recover the information about multiple roots and pairs of roots. (See [25]).

In Section 3, we give a procedure to obtain the limiting Newton diagram of a given polynomial. It is effective in the sense that, if we can bound the separation

max_{|ζ_{i}|>|ζ_{j}|} \frac{|ζ_{i}|}{|ζ_{j}|} ,

then we can effectively identify the multiple roots and pairs of roots. It will converge, and eventually provide the list of multiple roots and pairs for any circle-free polynomial, in a finite (but unknown, not effective) number of iterations.

The other big weakness of classical Graeffe iteration is the fact that it returns the moduli of the roots, but not the actual roots. As a matter of fact, information about the argument of the roots is lost, and should be recovered by other means:

Example 3. Consider polynomials f(x) = x^{2}−2x+1, g(x) = x^{2}−1, h(x) = x^{2}+1. After two Graeffe iterations, all the three polynomials are mapped into f(x).

Many algorithms have been proposed to recover the arguments [28]. In this paper, we will differentiate the Graeffe iteration operator, and obtain an iteration defined on the appropriate tangent bundle. This new operator will define a mapping between 1-jets of polynomials. By the latter we mean expressions of the form f(x) + ϵḟ(x), where ϵ is a formal parameter. This procedure is discussed in section 4. In the end of the same section, we shall discuss the stability properties of this process.

In section 5, we compare the numerical behavior of Renormalized Tangent Graeffe Iteration to other publicly available algorithms.

2 Renormalizing Graeffe

2.1 The Renormalized Graeffe Iteration

Example 2 shows a typical behavior of classical Graeffe iteration performed by digital computers [9]. In order to avoid that sort of overflow, the authors introduced in [22] the Renormalized Graeffe Iteration. Although the details and the mathematical foundations of the algorithm are described in [22], to keep the present work self-contained, we give below a very short description of the main ideas:

One should consider the computation of g = G^{N}f as divided in several renormalization levels.

Level 0 Coefficients of f
Level 1 Coefficients of Gf
Level 2 Coefficients of G^{2}f
⋮
Level N Coefficients of G^{N}f

At renormalization level N, all coefficients g_{j} of g = G^{N}f should be represented in coordinates

r_{j}^{(N)} = −2^{−N} log |g_{j}| ,

and

α_{j}^{(N)} = g_{j}/|g_{j}| .

Therefore, we shall obtain convergence of the radial coordinates r_{j}^{(N)} (at least in the case of roots of different moduli). The dynamics of the angular coordinates α_{j}^{(N)} is typically chaotic.

In order to pass from level N to level N +1, a Renormalized Graeffe Operator was defined in [22]. Intermediate computations were performed in coordinates r_{j}^{(N)} and α_{j}^{(N)} by means of renormalized arithmetic operations. For instance, the renormalized sum (r, α) of (r_{1}, α_{1}) and (r_{2}, α_{2}) can be defined (in renormalization level N) by

r = −2^{−N} log |α_{1}e^{−2^{N} r_{1}} + α_{2}e^{−2^{N} r_{2}}|

α = \frac{α_{1}e^{−2^{N} r_{1}} + α_{2}e^{−2^{N} r_{2}}}{|α_{1}e^{−2^{N} r_{1}} + α_{2}e^{−2^{N} r_{2}}|}

Renormalized sum can be computed without overflow by the formula in Algorithm 1. This is a simplified, non-optimal version of renormalized sum. Notice that one or two of the inputs can be the renormalization of 0, i.e., ∞. Under the usual conventions, ∞ is greater than any real number. Therefore, if only one of the arguments is ∞, the correct result will be returned.

A few extra mathematical ideas related to the renormalized Graeffe operators, as well as some other mathematical results can be found in [22].

2.2 The Renormalized Newton Diagram

The first goal of this section is to introduce the concept of Renormalized Newton diagram, which is going to play a key role in the practical implementation of the algorithm discussed in this paper. The second is to prove a convergence result based on such idea using some earlier results of Ostrowski’s.

We start by reviewing the concept of Newton diagram, which has been used extensively by Ostrowski, Puiseux and Dumas, among others.

Algorithm 1 RenSum ( r_{1}, α_{1}, r_{2}, α_{2}, p )

{ It is assumed that r_{1} and r_{2} are real numbers or +∞, and that |α_{1}| = |α_{2}| = 1. The number p should be equal to 2^{N}, where N is the renormalization level. This routine computes (in renormalized coordinates !) the sum of α_{1}e^{−pr_{1}} and α_{2}e^{−pr_{2}}. } if r_{1} = r_{2} = +∞ then

return +∞, 1 Δ ← r_{2} − r_{1} if Δ ≥ 0 then

t ← α_{1} + α_{2}e^{−pΔ}

return r_{1} − \frac{log(|t|)}{p}, t/|t| else

t ← α_{2} − α_{1}e^{pΔ}

return r_{2} − \frac{log(|t|)}{p}, t/|t|

Let

f = ∑_{i=0}^{d} f_{i}x^{i} ,

be a degree d polynomial. As before, we denote by g = G^{N}f the N-th Graeffe iterate of f.

We order the roots of f in nondecreasing order of their moduli, to wit:

|ζ_{1}| ≤ ⋯ ≤ |ζ_{d}| .  (1)

If the above inequalities are all strict, then

lim_{N→∞} 2^{−N} log (\frac{|g_{i}|}{|g_{i+1}|}) = log(|ζ_{i+1}|) .

For each N, consider the piecewise linear function r^{(N)} : [0, d] → ℝ ∪ {+∞} satisfying

r^{(N)}(i) =^{def} −2^{−N} log |g_{i}| .

Notice that under the above assumptions, r^{(N)} is convex for sufficiently large N. (See figure 1). Indeed, since |ζ_{i+1}| > |ζ_{i}|, the inclinations satisfy (for large N)

r^{(N)}(i + 1) − r^{(N)}(i) ≥ r^{(N)}(i) − r^{(N)}(i − 1) .

It is easy to see that if two consecutive roots, say ζ_{i} and ζ_{i+1}, have approximately the same absolute value, then the three corresponding points

(i − 1, r^{(N)}(i − 1)) , (i, r^{(N)}(i)) , (i + 1, r^{(N)}(i + 1))

Figure 1: The function r^{(N)}(i), for N = 0, 1, 2

will be approximately aligned. Furthermore, the functions r^{(N)} converge to a piecewise linear convex function.

However, if the inequalities in (1) are not strict, the functions r^{(N)} may fail to converge.

Example 4. Let f(x) = (x − 1)(x − e^{iθ}). Then its N-th Graeffe iterate is

g(x) = x^{2} − (1 + e^{2^{N} iθ})x + e^{2^{N} iθ}

Therefore, we have r^{(N)}(0) = r^{(N)}(2) = 0, but we also have

r^{(N)}(1) = −2^{−N} log |1 + e^{2^{N} iθ}| = −2^{−N} log |2 cos 2^{N−1}θ| .

Depending on the choice of θ, this last value can range anywhere from −2^{−N} log 2 to +∞.

This is one of the reasons for introducing the convex hull of each r^{(N)}, that will be subsequently called the Renormalized Newton Diagram. (See figure 2).

Our approach has the advantage of simplifying some of the arguments by Ostrowski in [25] by providing plain convergence of Renormalized Newton Diagrams; however, we will quote several of the results by Ostrowski in the sequel.

One of the major goals of Ostrowski in [25] was to obtain effective bounds for the moduli of roots. This was possible by introducing of the majorant of a given polynomial:

Definition 3. A majorant of a given polynomial is any other polynomial, of same degree, with nonnegative coefficients greater than or equal to the given polynomial’s coefficients.

Figure 2: The function r^{(N)}(i) and its Convex Hull

The first step of Ostrowski’s construction is Newton’s majorant:

Definition 4. A polynomial A = ∑_{i=0}^{d} A_{i}x^{i} with nonnegative coefficients is called normal if the following conditions hold:

1. If A_{i} > 0 and A_{j} > 0 for i > j, then A_{l} > 0 for all i < l < j.

2. For l = 1, …, d − 1,

A_{l}^{2} ≥ A_{l−1}A_{l+1} .

A normal majorant T = ∑ T_{i}x^{i} of f is called minimal if for any other majorant T′ of f we have

T_{j} ≤ T′_{j} , j = 0, 1, …, d.  (2)

Notice that Condition 2 above means that the graph of the points of the form (l, − log(r_{l})), for l = 0, …, d is convex.

In the language of majorants,

Proposition 1 (Ostrowski[25]). Any polynomial

f = ∑_{j=0}^{d} f_{j}x^{j}

possesses a unique minimal normal majorant,

M_{f} = ∑_{j=0}^{d} T_{j}x^{j} .

The polynomial M_{f} will be called Newton Majorant of the polynomial f.

The result can be proved by using the convex hull ϕ of the function − log |f_{i}|. and constructing the polynomial M_{f} as the polynomial with positive coefficients (M_{f})_{j} = e^{−ϕ(j)}. We refer the interested reader to Ostrowski’s work [25, 26].

We remark that if the polynomial has roots of strictly increasing moduli, then the coefficients T_{i} of the Newton’s majorant of g = G^{N}f coincide with |g_{i}| for N sufficiently large and i = 0, 1, …, d.

However, the introduction of the Newton Diagram allows us to consider the general situation of possibly many roots of same moduli. As before, we order them in nondecreasing order and consider the indices

i_{0} = 0 < i_{1} < i_{2} < ⋯ < i_{l} < i_{l+1} = d

as 0, d and exactly those integers i between 1 and d − 1 such that |ζ_{i}| < |ζ_{i+1}|. This way, we have that

|ζ_{i_{j−1}}| < |ζ_{i_{j−1}+1}| = ⋯ = |ζ_{i_{j}}| < |ζ_{i_{j}+1}| .

The fundamental result, in this case is

lim_{N→∞} 2^{−N} log \frac{|g_{i_{j}}|}{|g_{i_{j+1}}|} = (i_{j+1} − i_{j}) log |ζ_{i}| ,

for i_{j} < i ≤ i_{j+1} and 0 ≤ j ≤ l. (c.f. equation (79.8) of [26]). In the language of Renormalized Newton Diagrams, that very same equation can be written as:

log |ζ_{i}| = lim_{N→∞} \frac{r^{(N)}(i_{j+1}) − r^{(N)}(i_{j})}{i_{j+1} − i_{j}} , for i_{j} < i ≤ i_{j+1}.

As remarked by Ostrowski, the above formulae are only useful in the determination of the moduli of the roots if we know “a priori” the values i_{1} < i_{2} < ⋯ < i_{l}. This is obviously not the case in most applications. Instead, Ostrowski’s results provide effective bounds for the convergence of the r^{(N)}, and thus for the values of the |ζ_{i}|’ s.

Theorem 2 (Ostrowski[25],Theorem IX.3). Let

ϱ(ν) =^{def} 1 − 2^{−1/ν} ,

and

R_{ν}^{(N)} =^{def} \frac{T_{ν−1}^{(N)}}{T_{ν}^{(N)}} ,

then

ϱ(ν) < \frac{|ζ_{ν}|^{2^{N}}}{R_{ν}^{(N)}} < \frac{1}{ϱ(d − ν + 1)}  ν = 1, …, d.  (3)

As a consequence of the above estimate, Ostrowski gets the following bound

(2d)^{−2^{−N}} < \frac{|ζ_{ν}|}{(R_{ν}^{(N)})^{2^{−N}}} < (2d)^{2^{−N}}

Corollary 1. If r^{(N)}(i) denotes the i-th ordinate of the N-th Renormalized Newton Diagram, then

lim_{N→∞} (r^{(N)}(i) − r^{(N)}(i − 1)) = log |ζ_{i}|

for i = 1, …, d. Furthermore, the error is bounded from above by

2^{−N} log(2d) .

This is indeed a strong result, since nothing is assumed on the coefficients or the roots of the original polynomial. However, it is possible to get a better error bound by assuming a minimal separation on the moduli:

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

for some ϵ > 0.

We note that in the above formula, if ϵ is well defined (i.e. there are at least two roots of different modulus) then it is non-zero.

3 Computation of the Newton Diagram

3.1 Algorithm and Main Statements

The main issue in this section is the following:

We are given a certain polynomial g, obtained after a few Graeffe iterations of a polynomial f. The roots of g are Z_{1}, …, Z_{d}, and we order them so that:

|Z_{1}| ≤ |Z_{2}| ≤ ⋯ ≤ |Z_{d}|

We want to know which of the inequalities are strict. We do not know the actual value of the Z_{i}’s, we know only the coefficients of g, i.e., the symmetric functions of the Z_{i}’s.

We are also ready to assume that

R =^{def} min_{|Z_{i+1}|>|Z_{i}|} \frac{|Z_{i+1}|}{|Z_{i}|}

is a large real number. Indeed, if ζ_{1}, …, ζ_{d} are the roots of f, always ordered such as |ζ_{1}| ≤ ⋯ ≤ |ζ_{d}|, then

ρ =^{def} min_{|ζ_{i+1}|>|ζ_{i}|} \frac{|ζ_{i+1}|}{|ζ_{i}|}

is always strictly greater than one. Hence, given any A > 0, by performing N ≥ log_{2} \frac{log A}{log ρ} iterations, we can assume that R = ρ^{2^{N}} ≥ A.

Recall that the Renormalized Newton Diagram of g is the convex hull of the function i ↦ −2^{−N} log g_{i}. As N grows, the Renormalized Newton Diagram of g converges to the convex hull of

i ↦ − log |f_{0}| + ∑_{j≤i} log |ζ_{i}| .

However, we want to be able to decide in finite time what are the sharp corners of the convex hull of i ↦ log |f_{d}| + ∑_{j≤i} log |ζ_{j}|). As in the preceding section, we write:

r_{i} = −2^{−N} log |g_{i}|

where g is the N-th iterate of f. Notice that we dropped the superscript N of r_{i}^{(N)}.

Proposition 2. Let g = G^{N}f, where f is a degree d polynomial and G denotes the Graeffe iteration. Let Z, ζ and ρ be as above. Assume that

N > 3 + log_{2} \frac{d log 2}{log ρ} .

Then, Algorithm 2 below with input N, d, r, ρ produces the list of the sharp edges of the convex hull of i ↦ ∑_{j≤i} log |ζ_{i}|).

Remark: Algorithm 2 has running time O(d).

3.2 Some Estimates about Symmetric Functions

In order to prove Proposition 2, we need a few estimates about symmetric functions. First of all, let I = {i : |Z_{i}| < |Z_{i+1}|} ∪ {0, d} be the set of sharp corners of the limiting Renormalized Newton Diagram. As before, let σ_{k} denote the k-th elementary symmetric function,

σ_{k}(Z) = ∑_{j_{1}<···<j_{k}} Z_{j_{1}} … Z_{j_{k}}

Then,

Lemma 1. For i ∈ I we have

σ_{d−i}(Z) = Z_{i+1}Z_{i+2} … Z_{d}(1 + c) ,

where |c| ≤ (\binom{d}{i} − 1) R^{−1} ≤ 2^{d}R^{−1}

Proof of Lemma 1: Write

σ_{d−i}(Z) = ∑_{j_{1}<···<j_{d−i}} Z_{j_{1}} … Z_{j_{d−i}}

Algorithm 2 Strict Convex Hull ( N, d, r, ρ )

{ Create a list Λ, containing initially the element Λ_{0} = 0} j ← 0; Λ_{j} ← 0 ;

{ The error bound below will follow from Lemma 4.} R ← ρ^{2^{N}} E ← \frac{1}{2} (2^{−N+1} log(2^{d} + 2^{d}R^{−1}) − 2^{−N+2} log(1 − 2^{d}R^{−1}) + \frac{log ρ}{2})

{ Now, we will try to add more points to the list Λ. At each step, we want to ensure that we have always a convex set.}

for i ← 1 to d do

{ We discard all the points in Λ that are external to the convex hull of Λ and the new point. Let Λ_{j} be the last element of Λ}

while j > 0 and \frac{r_{Λ_{j}} −r_{Λ_{j−1}}}{Λ_{j}−Λ_{j−1}} > \frac{r_{i}−r_{Λ_{j}}}{i−Λ_{j}} − E do j ← j − 1

{ Now, we append the point i}

j ← j + 1

Λ_{j} ← i;

Return (Λ_{0}, ⋯, Λ_{j})

In the sum above, |Z_{j_{1}} … Z_{j_{d−i}}| ≤ R^{−1}|Z_{i+1}Z_{i+2} … Z_{d}| for any choice of j_{1}, … j_{d−i} except i + 1, … d. Since there are \binom{d}{i} − 1 other terms, we obtain that:

|σ_{d−i}(Z) − Z_{i+1}Z_{i+2} … Z_{d}| ≤ (\binom{d}{i} − 1) R^{−1}|Z_{i+1}Z_{i+2} … Z_{d}|

Definition 5. We will say that i_{1} and i_{2} are successive elements of I if and only if:

1. i_{1} ∈ I

2. i_{1} < i < i_{2} ⇒ i ∉ I

3. i_{2} ∈ I

Lemma 2. Let i_{1} and i_{2} be successive elements of I, and let i_{1} < l < i_{2}. Then

|σ_{d−l}(Z)| ≤ (\binom{i_{2} − i_{1}}{i_{2} − l} + c′) |Z_{i_{2}}|^{i_{2}−l}|Z_{i_{2}+1}| … |Z_{d}|

with c′ ≤ (\binom{d}{i} − \binom{i_{2}−i_{1}}{i_{2}−l}) R^{−1} < 2^{d}R^{−1}

Proof of Lemma 2: Write

σ_{d−l}(Z) = ∑_{j_{1}<···<j_{d−l}} Z_{j_{1}} … Z_{j_{d−l}}

= ∑′Z_{j_{1}} … Z_{j_{d−l}} + ∑′′Z_{j_{1}} … Z_{j_{d−l}}

where ∑′ ranges over the j such that i_{1} < j_{1} < ⋯ < j_{i_{2}−l} < i_{2} + 1 and j_{i_{2}−l+1} = i_{2} + 1, …, j_{d−l} = d. Of course, ∑′′ ranges over all the other terms.

We can rewrite ∑′ as:

∑′ = Z_{i_{2}+1}Z_{i_{2}+2} … Z_{d} ( ∑_{i_{1}<j_{1}<···<j_{i_{2}−l}≤i_{2}} Z_{j_{1}} … Z_{j_{i_{2}−l}} )

Hence,

|∑′| ≤ \binom{i_{2} − i_{1}}{i_{2} − l}|Z_{i_{2}}|^{i_{2}−l}|Z_{i_{2}+1}| … |Z_{d}|

The terms in ∑′′ are all smaller than R^{−1}|Z_{i_{2}}|^{i_{2}−l}|Z_{i_{2}+1}| … |Z_{d}|. Since there are \binom{d}{i} − \binom{i_{2}−i_{1}}{i_{2}−l} of them,

|∑′′| < (\binom{d}{i} − \binom{i_{2} − i_{1}}{i_{2} − l}) R^{−1}|Z_{i_{2}}|^{i_{2}−l}|Z_{i_{2}+1}| … |Z_{d}|

Adding those two bounds, we obtain indeed:

|σ_{d−l}(Z)| ≤ (\binom{i_{2} − i_{1}}{i_{2} − l} + (\binom{d}{i} − \binom{i_{2} − i_{1}}{i_{2} − l}) R^{−1}) |Z_{i_{2}}|^{i_{2}−l}|Z_{i_{2}+1}| … |Z_{d}|

The estimates above can be converted into ‘logscale’ estimates:

Lemma 3. Let i_{1}, i_{2} be successive elements of I, and let i_{1} < l < i_{2}. Then the following three equations are true:

1.

\frac{r(i_{2}) − r(i_{1})}{i_{2} − i_{1}} = log |ζ_{i_{2}}| + 2^{−N+1} log(1 + c)

with |c| ≤ (max (\binom{d}{i_{1}}, \binom{d}{i_{2}}) − 1) R^{−1} < 2^{d}R^{−1}.

2.

\frac{r(i_{2}) − r(l)}{i_{2} − l} ≤ log |ζ_{i_{2}}| + 2^{−N} log (\binom{i_{2} − i_{1}}{i_{2} − l} + c′) + 2^{−N} log |1 + c|

where |c| ≤ (\binom{d}{i_{2}} − 1) R^{−1} ≤ 2^{d}R^{−1} and c′ ≤ (\binom{d}{i} − \binom{i_{2}−i_{1}}{i_{2}−l}) R^{−1} <

2^{d}R^{−1}.

\frac{r(l) − r(i_{1})}{l − i_{1}} ≥ log |ζ_{i_{2}}| − 2^{−N} log (\binom{i_{2} − i_{1}}{i_{2} − l} + c′) + 2^{−N} log |1 + c|

where |c| ≤ (\binom{d}{i_{1}} − 1) R^{−1} ≤ 2^{d}R^{−1} and c′ ≤ (\binom{d}{i} − \binom{i_{2}−i_{1}}{l−i_{1}}) R^{−1} <

2^{d}R^{−1}.

Proof of Lemma 3: By using Lemma 1 with i = i_{2}, we obtain:

−2^{−N} log |σ_{d−i_{2}}(Z)| = −2^{−N} (log |Z_{i_{2}+1}| + … log |Z_{d}|) − 2^{−N} log(|1 + c′′|) (4)

Using the same lemma with i = i_{1}, we get:

−2^{−N} log |σ_{d−i_{1}}(Z)| = −2^{−N} (log |Z_{i_{1}+1}| + … log |Z_{d}|) − 2^{−N} log(|1 + c′′|) (5)

Subtracting the two previous expressions and dividing by i_{2} − i_{1} we get:

\frac{r(i_{2}) − r(i_{1})}{i_{2} − i_{1}} = 2^{−N} log |Z_{i_{2}}| + 2^{−N+1} log(1 + c′′)
= log |ζ_{i_{2}}| + 2^{−N+1} log(1 + c′′)

This shows the first part of the Lemma.

By using Lemma 2, we can also bound:

−2^{−N} log |σ_{d−l}(Z)| ≥ −2^{−N}(i_{2} − l) log |Z_{i_{2}}|
−2^{−N} (log |Z_{i_{2}+1}| + … log |Z_{d}|)
−2^{−N} log(|\binom{i_{2}−i_{1}}{i_{2}−l} + c′|)

(6)

where c′ is as in Lemma 2.

We can now estimate equation (4) minus equation (6), altogether divided by i_{2} − l:

\frac{r(i_{2}) − r(l)}{i_{2} − l} ≤ log |ζ_{i_{2}}| + 2^{−N} log (\binom{i_{2} − i_{1}}{i_{2} − l} + c′) + 2^{−N} log |1 + c|

We can also estimate equation (6) minus equation (5), altogether divided by l − i_{1}:

\frac{r(l) − r(i_{1})}{l − i_{1}} ≥ log |ζ_{i_{2}}| − 2^{−N} log (\binom{i_{2} − i_{1}}{i_{2} − l} + c′) + 2^{−N} log |1 + c|

3.3 A Decision Criterion

Lemma 3 can be used do decide if a point in the convex hull of g is converging to a sharp corner of the limiting convex hull or not. Lemma 4. Assume that a. m ≥ max(i_{2} − i_{1}) when i_{1} and i_{2} are successive elements of I. b. 2^{−N+1} log(2^{m} + 2^{d}R^{−1}) − 2^{−N+2} log(1 − 2^{d}R^{−1}) < E < \frac{log ρ}{2} c. i < j < k

Then,

1. If i and j are successive elements of I and there is no other element of I between j and k, then \frac{r(j)−r(i)}{j−i} < \frac{r(k)−r(j)}{k−j} − E

2. If i and k are successive elements of I then \frac{r(j)−r(i)}{j−i} > \frac{r(k)−r(j)}{k−j} − E Proof of Lemma 4: Part 1: Assume that i, j are successive elements of I. Then, part 1 of Lemma 3 implies:

\frac{r(j) − r(i)}{j − i} ≤ log |ζ_{j}| − 2^{−N+1} log(1 − 2^{d}R^{−1})

For the evaluation of \frac{r(k)−r(j)}{k−j}, we have to distinguish two cases: If k ∈ I, then

\frac{r(k) − r(j)}{k − j} ≥ log |ζ_{k}| + 2^{−N+1} log(1 − 2^{d}R^{−1})

If k ∉ I, let m be such that j and m are successive elements of I. Recall that j < k < m by hypothesis. Using part 3 of Lemma 3, we get:

\frac{r(k) − r(j)}{k − j} ≥ log |ζ_{m}| − 2^{−N} log(2^{m} + 2^{d}R^{−1}) + 2^{−N} log(1 − 2^{d}R^{−1})

In any case,

\frac{r(j) − r(i)}{j − i} ≤ \frac{r(k) − r(j)}{k − j} + log |ζ_{j}| − log |ζ_{k}|
+2^{−N} log(2^{m} + 2^{d}R^{−1}) − 2^{−N+2} log(1 − 2^{d}R^{−1})

We use the hypothesis E < \frac{log ρ}{2} to deduce that log |ζ_{j}| − log |ζ_{k}| + E < −E, and:

\frac{r(j) − r(i)}{j − i} < \frac{r(k) − r(j)}{k − j} − E

Part 2: Using Lemma 3, we have:

\frac{r(j) − r(i)}{j − i} ≥ log |ζ_{k}| − 2^{−N} log(2^{m} + 2^{d}R^{−1}) + 2^{−N} log(1 − 2^{d}R^{−1})

\frac{r(k) − r(j)}{k − j} ≤ log |ζ_{k}| + 2^{−N} log(2^{m} + 2^{d}R^{−1}) − 2^{−N} log(1 − 2^{d}R^{−1})

Subtracting, we obtain:

\frac{r(j) − r(i)}{j − i} ≥ \frac{r(k) − r(j)}{k − j} − 2^{−N+1} log(2^{m} + 2^{d}R^{−1})
+2^{−N+1} log(1 − 2^{d}R^{−1})
> \frac{r(k) − r(j)}{k − j} − E

3.4 Proof of Proposition 2

Proof of Proposition 2: Let N > 3 + log_{2} \frac{d log 2}{log ρ}. It is easy to check that R > 2^{8d}, hence 2^{d}R^{−1} < 2^{−7d}. So we can bound:

2^{−N+1} log(2^{m} + 2^{d}R^{−1}) − 2^{−N+2} log(1 − 2^{d}R^{−1}) <

< \frac{2(m + 1) log ρ log 2}{8d log 2} + \frac{4 log ρ}{8d log 2} \frac{1}{2^{8} − 1} < \frac{log ρ}{2}

Therefore, we are in the conditions 1 and 2 of Lemma 4, with m = d. Correctness of the algorithm 2 can be proved now by induction.

Induction Hypothesis . At step i, the list Λ contains Λ_{0}, …, Λ_{s}, Λ_{s+1}, Λ_{j} where Λ_{0}, … Λ_{s} are all successive elements of I and Λ_{s+1}, …, Λ_{j} are not in I. (Possibly, we can have s = j).

The induction hypothesis is true at step 1, with j = 0, and 0 ∈ I. At each step, there are two possibilities:

Case 1: i ∉ I. In that case a few of the Λ_{s+1}, … Λ_{j} may be discarded; but part 1 of Lemma 4 prevents the algorithm from discarding elements of I.

Case 2: i ∈ I. In that case, part 2 of Lemma 4 guarantees that all the Λ_{s+1}, … Λ_{j} will be discarded.

Hence, the induction hypothesis is true at step i + 1. At step d, the last point d is added to Λ. Since d ∈ I, Λ = I.

A note on the running time: although the usual complexity of a convex hull algorithm is O(d log d) for d points in the plane, the complexity is smaller when those points are ‘ordered’ like ours: (i, r(i)). (Compare with Theorem 4.12 in [30]). Algorithm 2 has a running time of O(d) operations (including a fixed number of transcendental operations). Indeed, each point is added to the list Λ precisely one time. It can be discarded only once, so the interior ‘while’ loop is executed at most d − 1 times in one execution of the algorithm.

4 Tangent Graeffe Iteration

4.1 Perturbation Methods, Infinitesimals, 1-Jets of Polynomials

Graeffe iteration provides the absolute values of each root in the case such roots are all of different moduli. Recovering the actual value of each root, and recovering pairs of conjugate roots or multiple roots require further work.

Many algorithms have been proposed to recover the actual roots, such as reverse Graeffe iteration, splitting algorithms. See [28] and references therein.

A possibility of theoretical interest would be to consider a perturbation of f; assume first that f is a polynomial with roots ζ_{1}, …, ζ_{d} such that |ζ_{1}| < |ζ_{2}| < ⋯ < |ζ_{d}|. Then, consider also the iterates of

f(x + ϵ)

Graeffe iteration of f(x) will provide |ζ_{1}|, …, |ζ_{d}|, while Graeffe iteration of f(x + ϵ) will provide |ζ_{1} − ϵ|, …, |ζ_{d} − ϵ|. Therefore, we will be able to compute:

|ζ_{i}|^{2} − |ζ_{i} − ϵ|^{2} = −2ϵRe (ζ_{i}) + ϵ^{2}

thus recovering ζ_{i}.

As mentioned before, this is a possibility of theoretical interest only. The perturbation method above would lose half of the working precision in any reasonable implementation. Therefore, we will prefer to compute the derivative of |ζ − ϵ|^{2} with respect to ϵ.

The value ϵ will be treated as an infinitesimal; therefore, instead of storing in memory a certain value z + ϵż, we will store z and ż separately. When computing some differentiable function G(z + ϵż), we will obtain a result G(x) + ϵDG(z)ż. So we will compute G(z) and DG(z)ż, but we will never need to assign an actual value to ϵ.

A quantity of the form z + ϵż is called a 1-jet. It can also be interpreted as an element of the tangent bundle of the manifold where z is supposed to live.

In this paper, we are mainly concerned with 1-jets of polynomials. We will represent degree d polynomials as points in ℝ^{d+1} or ℂ^{d+1}. Therefore, a 1-Jet of polynomials can be represented as a point of ℝ^{2d+2} or ℂ^{2d+2}, since we are working with a linear space.

The dot notation (such as in ż) will be reserved in this paper to the ‘tangent’ coordinate of a 1-jet z + ϵż. We reserve the notation f′ to the derivative \frac{∂}{∂x}f of a univariate function f = f(x), and the notation DF to the derivative of a multivariate function F. We need the following construction from Calculus on Manifolds [1, 18]: Let G be a differentiable function from manifold X into manifold Y. Its tangent map can be written, in our 1-Jet notation, as:

\begin{aligned}T G : & T X & → & T Y \\ & f + ϵḟ & ↦ & G(f) + ϵDG_{f} ḟ, \end{aligned}

where as usual T X and T Y denote the tangent bundle of X and Y respectively.

The iteration of the 1-jet of polynomials f + ϵf′ can be used to recover the actual value of the roots of f. For instance:

Example 5. Let f be a real circle free polynomial, not vanishing at 0. Consider the 1-jet f(x + ϵ) = f(x) + ϵf′(x); its solutions are ζ_{j} − ϵ, where ζ_{j} are the roots of f. Let g + ϵġ = T G^{N}(f + ϵf′). Lemma 6 below will imply, in the particular case ζ_{j} is a real isolated root, that:

ζ_{j} = lim_{N→∞} 2^{−N} (\frac{|g_{j}|}{|g_{j−1}|})^{\frac{1}{2^{N−1}}} (\frac{ġ_{j}}{g_{j}} − \frac{ġ_{j−1}}{g_{j−1}} )

In case ζ_{j} and ζ_{j+1} = ζ̄_{j} are an isolated pair of conjugate roots, the limit will be:

Re ζ_{j} = lim_{N→∞} 2^{−N−1} (\frac{|g_{j+1}|}{|g_{j−1}|})^{\frac{1}{2^{N}}} (\frac{ġ_{j+1}}{g_{j+1}} − \frac{ġ_{j−1}}{g_{j−1}} )

In the following section, we compute the tangent map of the Graeffe operator in usual and renormalized coordinates.

4.2 The Iteration

Let f + ϵḟ be a 1-Jet of polynomials. Then its Tangent Graeffe Iterate is:

T G(f(x) + ϵḟ(x)) = (−1)^{d} (f(\sqrt{x}) + ϵḟ(\sqrt{x})) (f(−\sqrt{x}) + ϵḟ(−\sqrt{x}))

This can be rewritten as:

T G(f(x) + ϵḟ(x)) = G(f) + (−1)^{d}ϵ (f(\sqrt{x})ḟ(−\sqrt{x}) + f(−\sqrt{x})ḟ(\sqrt{x}))

Precise formulae for computing g + ϵġ = T G(f + ϵḟ) are:

\{ \begin{matrix}g_{i} & = & (−1)^{d+i}f_{i}^{2} + 2 ∑_{j≥1}(−1)^{d+i+j}f_{i+j}f_{i−j} \\ ġ_{i} & = & 2 ∑_{j}(−1)^{d+i+j}f_{i−j}ḟ_{i+j} \end{matrix}

For an efficient root-finding algorithm the equations above need to be renormalized. At each step, this is done by replacing products and sums by their renormalized counterparts. An adjustment is necessary to pass from one renormalization level to another (division of the coordinates r by 2). Those adjustments are summarized in Algorithm 3

4.3 Convergence Results

It is now time to show convergence of the (Renormalized) Tangent Graeffe Operator. Assume one is given a circle free polynomial f. Its roots will be ordered as |ζ_{1}| ≤ |ζ_{2}| ≤ ⋯ |ζ_{d}|. One can use the (Renormalized) Newton diagram to collect together the roots with same moduli. Those will represent single roots, multiple roots, or (in the case of real polynomials) pairs of conjugate roots or pairs of multiple conjugate roots.

Algorithm 3 TangentGraeffe (N, d, r, α, r̂, α̂)

{ N (Renormalization level) and d (degree) are integers; r and r̂ should be real arrays, and α and α̂ should be array of modulus one complex numbers. This routine computes, in renormalized coordinates, the Tangent Graeffe Iterate of the 1-jet: ∑_{i} (α_{i}e^{−2^{N} r_{i}} + ϵα̂_{i}e^{−2^{N} r̂_{i}}) x^{i}. The coordinates of the result are given in renormalization level N + 1.} p ← 2^{N+1} for i ← 0 to d do

(s_{i}, β_{i}) ← (r_{i}, (−1)^{i}α_{i}^{2})

(ŝ_{i}, β̂_{i}) ← ((r_{i} + r̂_{i})/2 − log(2)/p, (−1)^{i}α_{i}α̂_{i})

for j ← 1 to min(d − i, i) do (s_{i}, β_{i}) ← RenSum(s_{i}, β_{i}, (r_{i+j} + r_{i−j})/2 + log(2)/p, (−1)^{i+j}α_{i+j}α_{i−j}, p) (ŝ_{i}, β̂_{i}) ← RenSum(ŝ_{i}, β̂_{i}, (r_{i+j} + r̂_{i−j})/2 + log(2)/p, (−1)^{i+j}α_{i+j}α̂_{i−j}, p) (ŝ_{i}, β̂_{i}) ← RenSum(ŝ_{i}, β̂_{i}, (r_{i−j} + r̂_{i+j})/2 + log(2)/p, (−1)^{i−j}α_{i−j}α̂_{i+j}, p)

return (s, β, ŝ, β̂)

Lemma 5 (Complex case). Let f be a complex circle-free polynomial with roots ζ_{1} ≠ 0, …, ζ_{d} ordered as in Theorem 1. Let

ρ =^{def} min_{|ζ_{i+1}|>|ζ_{i}|} \frac{|ζ_{i+1}|}{|ζ_{i}|}

and let g + ϵġ = (T G)^{N}(f + ϵf′). Suppose that j and j + d′ are successive elements of I = {i : |ζ_{i}| < |ζ_{i+1}|} ∪ {0; d}. Then,

lim_{N→∞} −\frac{2^{−N}}{d′} \overline{(\frac{ġ_{j+d′}}{g_{j+d′}} − \frac{ġ_{j}}{g_{j}} )} = \frac{ζ_{j+d′}}{|ζ_{j+d′}|^{2}} .

Furthermore, the error is bounded by:

2^{d+3} \frac{d}{d′} \frac{|ζ_{d}|}{|ζ_{1}|}ρ^{−2^{N}} |ζ_{j+d′}|^{−1}

Lemma 5 will be proved in Subsection 4.5.

Also, real polynomials have usually pairs of conjugate roots; they may have pairs with multiplicity. In that case, we can show that:

Lemma 6 (Real case). Let f be a real circle-free polynomial with roots ζ_{1} ≠ 0, …, ζ_{d} ordered as in Theorem 1. Let

ρ =^{def} min_{|ζ_{i+1}|>|ζ_{i}|} \frac{|ζ_{i+1}|}{|ζ_{i}|}

and let g + ϵġ = (T G)^{N}f + ϵf′. Suppose that j and j + d′ are successive elements of I = {i : |ζ_{i}| < |ζ_{i+1}|} ∪ {0; d}. Then,

lim_{N→∞} −\frac{2^{−N}}{d′} (\frac{ġ_{j+d′}}{g_{j+d′}} − \frac{ġ_{j}}{g_{j}} ) = \frac{Reζ_{j+d′}}{|ζ_{j+d′}|^{2}} .

Furthermore, the error is bounded by:

2^{d+3} \frac{d}{d′} \frac{|ζ_{d}|}{|ζ_{1}|}ρ^{−2^{N}} |ζ_{j+d′}|^{−1}

The proof of Lemma 6 is also postponed to subsection 4.5. Lemmas 5 and 6 can be used to recover the roots of a polynomial from the Tangent Graeffe iterates of its 1-jet:

Algorithm 4 RealRecover ( N, d, I, r, α, r̂, α̂ )

{ This procedure attempts to recover the roots of the degree d real polynomial ∑_{i} α_{i}e^{−2^{N} r_{i}}x^{i}. The list of sharp corners of its Newton Diagram is supposed given in I = (I_{0}, ⋯, I_{1+size(I)}). See Lemma 6 for a justification } for k ← 0 to Size(I) do

d′ ← I_{k+1} − I_{k} (b, β) ← RenSum(r̂_{I_{k+1}} − r_{I_{k+1}}, \frac{α̂_{I_{k+1}}}{α_{I_{k+1}}} , r̂_{I_{k}} − r_{I_{k}}, −\frac{α̂_{I_{k}}}{α_{I_{k}}} , 2^{N})

m ← exp (2\frac{r_{I_{k+1}} −r_{I_{k}}}{d′} )

x ← −β \frac{2^{−N}}{d′} m exp −2^{N}b if I_{k+1} − I_{k} is even and m > |x|^{2} then

y ← \sqrt{m − |x|^{2}} else

x ← m\frac{x}{|x|}

y ← 0 for j ← 0 to I_{k+1} − I_{k} − 1 do

ζ_{I_{k}+j+1} = x + (−1)^{j}y return ζ

Algorithm 5 ComplexRecover ( N, d, I, r, α, r̂, α̂ )

{ This procedure attempts to recover the roots of the degree d complex polynomial ∑_{i} α_{i}e^{−2^{N} r_{i}}x^{i}. The list of sharp corners of its Newton Diagram is supposed given in I. See Lemma 5 for a justification } for k ← 0 to Size(I) do

d′ ← I_{k+1} − I_{k} (b, β) ← RenSum(r̂_{I_{k+1}} − r_{I_{k+1}}, \frac{α̂_{I_{k+1}}}{α_{I_{k+1}}} , r̂_{I_{k}} − r_{I_{k}}, −\frac{α̂_{I_{k}}}{α_{I_{k}}} 2^{N}, )

m ← exp (2\frac{r_{I_{k+1}} −r_{I_{k}}}{d′} )

x ← −β̄\frac{2^{−N}}{d′} m exp −2^{N}b for j ← 0 to I_{k+1} − I_{k} − 1 do

ζ_{I_{k}+j+1} = x return ζ

4.4 The Main Algorithm

We can now state the algorithm of Theorem 1. We start with a fixed, arbitrary value for ρ(f) = max_{|ζ_{i}|>|ζ_{j}|} \frac{|ζ_{i}|}{|ζ_{j}|}. Proposition 2 guarantees that if

N > 3 + log_{2} \frac{d log 2}{log ρ(f)} ,

then after the N-th iterate the convex hull of the Newton Diagram of f is computed correctly.

Algorithm 6 Solve (d, f, isreal)

{ It is assumed here that f is a degree d, circle-free real or complex polynomial. In the general case, one should first find and output the trivial (0 and ∞) roots of f, then deflate f. After that, one should perform a random real (resp. complex) conformal transform on f so it becomes circle-free }

for i ← 0 to d do

if f_{i} ≠ 0 then

α_{i} ← f_{i}/|f_{i}| else

α_{i} ← 1 r_{i} ← − log |f_{i}|

for i ← 0 to d − 1 do f′_{i} ← (i + 1)f_{i+1} if f′_{i} ≠ 0 then

α̂_{i} ← f′_{i}/|f′_{i}| else

α̂_{i} ← 1 r̂_{i} ← − log |f′_{i}|

N ← 0 ρ = 2 loop

r, α, r̂, α̂ ← TangentGraeffe (N, d, r, α, r̂, α̂) N ← N + 1 I ← Convex Hull (N, d, r, ρ) if isreal then

ζ ← RealRecover ( N, d, I, r, α, r̂, α̂ ) else

ζ ← ComplexRecover ( N, d, I, r, α, r̂, α̂ ) Output ζ_{1}, ⋯ ζ_{d} if N > 3 + log_{2} \frac{d log 2}{log ρ} then

{ Proposition 2 implies that at this point, I is indeed correct for all the polynomials with separation ration ≥ ρ. Therefore, it is time to decrease ρ. } ρ ← \sqrt{ρ}

After the N-th iteration, convergence is guaranteed by the following bounds: According to Lemma 3, at the execution of algorithm Complex Recover (resp. Real Recover),

|exp (\frac{2}{i_{k+1} − i_{k}} (g_{i_{k+1}} − g_{i_{k}}))| = |ζ_{i_{k+1}}|(1 + δ_{1})

where |δ_{1}| ≤ e^{2^{−N+1} log 1+2^{d}ρ^{−2^{N}}}.

Introducing the error bound of Lemma 5, one gets in the complex case:

−\frac{2^{−N}}{d′} \overline{(a − b)} = \frac{ζ_{i_{k+1}}}{|ζ_{i_{k+1}}|^{2}} (1 + δ_{2})

with |δ_{2}| < 2^{d+3} \frac{d}{d′} \frac{|ζ_{d}|}{|ζ_{1}|}ρ^{−2^{N}}, and where a and b are as in the Algorithms. Therefore,

|ζ_{i_{k}+j+1} − x̄| ≤ δ|ζ_{i_{k}+j+1}|

where δ < (1 + δ_{1})^{2}(1 + δ_{2}) − 1.

The real case is analogous. According to Lemma 6,

−\frac{2^{−N}}{d′} (a − b) = \frac{Re ζ_{i_{k+1}}}{|ζ_{i_{k+1}}|^{2}} (1 + δ_{2})

so that

|Re ζ_{i_{k}+j+1} − x| < δ|ζ_{i_{k}+j+1}|

and

|Im ζ_{i_{k}+j+1} − \sqrt{1 − x^{2}}| < δ′|ζ_{i_{k}+j+1}|

where δ′ = δ + O(δ^{2})

In both cases, |δ| and eventually |δ′| are dominated by ρ^{−2^{N}}. Since

Aρ^{−2^{N}} = 2^{log_{2} A−2^{N} log_{2} ρ}2^{−2^{N−C}} ,

the bound in Theorem 1 follows.

4.5 Proof of Lemmas 5 and 6

Consider the 1-jet of degree d polynomials f + ϵḟ, with solutions ζ_{i} + ϵζ̇_{i}. After N steps of Tangent Graeffe Iteration, we obtain a 1-jet of polynomials

g + ϵġ = (T G)^{N} (f + ϵḟ) .

Since the differential of the transformation P^{N} : ζ ⟼ z = ζ^{2^{N}} gives

DP^{N} : ζ + ϵζ̇ ⟼ (ζ_{j})^{2^{N}} + ϵ2^{N} (ζ_{j})^{2^{N}} \frac{ζ̇_{j}}{ζ_{j}} ,

it is clear that the roots Z_{j} + ϵŻ_{j} of g + ϵġ will be

(ζ_{j})^{2^{N}} + ϵ2^{N} (ζ_{j})^{2^{N}} \frac{ζ̇_{j}}{ζ_{j}} .

We can now compute the derivative at ϵ = 0 of g_{i} + ϵg_{i} = σ_{d−i}(Z + ϵŻ). Let’s denote by Zĵ the vector (Z_{1}, ⋯, Ẑ_{j}, ⋯, Z_{d}) ∈ ℂ^{d−1}. Take i ∈ I and notice that

\frac{∂}{∂ϵ}|_{ϵ=0}σ_{d−i}(Z + ϵŻ) = ∑_{j} Ż_{j}σ_{d−i−1}(Zĵ)

= ∑_{j} \frac{Ż_{j}}{Z_{j}} Z_{j}σ_{d−i−1}(Zĵ)

Thus,

\frac{∂}{∂ϵ}|_{ϵ=0} \frac{σ_{d−i}(Z + ϵŻ)}{σ_{d−i}(Z)} = ∑_{j} \frac{Ż_{j}}{Z_{j}} \frac{Z_{j}σ_{d−i−1}(Zĵ)}{σ_{d−i}(Z)}

= (∑_{j>i} + ∑_{j≤i}) \frac{Ż_{j}}{Z_{j}} \frac{Z_{j}σ_{d−i−1}(Zĵ)}{σ_{d−i}(Z)} .

Due to Lemma 1

∑_{j>i} \frac{Ż_{j}}{Z_{j}} \frac{Z_{j}σ_{d−i−1}(Zĵ)}{σ_{d−i}(Z)} = ∑_{j>i} \frac{Ż_{j}}{Z_{j}} (1 + η_{j})

and

∑_{j≤i} \frac{Ż_{j}}{Z_{j}} \frac{Z_{j}σ_{d−i−1}(Zĵ)}{σ_{d−i}(Z)} = ∑_{j≤i} \frac{Ż_{j}}{Z_{j}} \frac{Z_{j}}{Z_{i+1}} (1 + η_{j}) ,

where

|η_{j}| ≤ 4\binom{d}{i}R^{−1} .

Although this estimate is by no means sharp, it suffices for our purposes. Using that i ∈ I and hence |z_{j}/z_{i+1}| < R^{−1} for i ≤ j, we get

| \frac{∂}{∂ϵ} |_{ϵ=0} \frac{σ_{d−i}(z + ϵż)}{σ_{d−i}(z)} − ∑_{j>i} \frac{ż_{j}}{z_{j}} | ≤ max | \frac{ż_{j}}{z_{j}} | \frac{4\binom{d}{i}(d − i) + i (1 + 4\binom{d}{i}R^{−1})}{R}
< max | \frac{ż_{j}}{z_{j}} |d2^{d+2}R^{−1}

Therefore, if we take the logarithmic derivative of the expression

(g + ϵġ) = (T G)^{N}(f + ϵḟ)

and evaluate at a successive couple of elements i_{1}, i_{2} ∈ I we get

| (\frac{ġ_{i_{1}}}{g_{i_{1}}} − \frac{ġ_{i_{2}}}{g_{i_{2}}} ) + ∑_{i_{1}<j≤i_{2}} 2^{N} \frac{1}{ζ_{j}} | < \frac{2^{N+d+3}d}{R} max_{l} |\frac{1}{ζ_{l}} | .

Now let’s assume we are under the hypothesis of Lemma 5. Then, since i_{1} and i_{2} are consecutive, it follows that

∑_{i_{1}<j≤i_{2}} ζ_{j}^{−1} = ∑_{i_{1}<j≤i_{2}} \frac{\overline{ζ_{j}}}{|ζ_{j}|^{2}} = (i_{2} − i_{1})\frac{\overline{ζ_{i_{2}}}}{|ζ_{i_{2}}|^{2}}

and so

|\frac{2^{−N}}{i_{2} − i_{1}} (\frac{ġ_{i_{2}}}{g_{i_{2}}} − \frac{ġ_{i_{1}}}{g_{i_{1}}} ) |ζ_{i_{2}}|^{2} − \overline{ζ_{i_{2}}} ||ζ_{i_{2}}|^{−1} < \frac{2^{d+3}d}{R(i_{2} − i_{1})} max_{r,s} \frac{|ζ_{r}|}{|ζ_{s}|} .  (7)

On the other hand, if we are under the hypothesis of Lemma 6, we get

∑_{i_{1}<j≤i_{2}} ζ_{j}^{−1} = (i_{2} − i_{1})\frac{Re ζ_{i_{2}}}{|ζ_{i_{2}}|^{2}} .

Thus,

|\frac{2^{−N}}{i_{2} − i_{1}} (\frac{ġ_{i_{2}}}{g_{i_{2}}} − \frac{ġ_{i_{1}}}{g_{i_{1}}} ) |ζ_{i_{2}}|^{2} − Re ζ_{i_{2}} ||ζ_{i_{2}}|^{−1} < \frac{2^{d+3}d}{R(i_{2} − i_{1})} max_{r,s} \frac{|ζ_{r}|}{|ζ_{s}|} . (8)

In either case, we have that each right hand side of equations (7) and (8) is bounded by

2^{d+3} \frac{d}{R(i_{2} − i_{1})} \frac{|ζ_{d}|}{|ζ_{1}|} .

Now, using the fact that the separation radius at the N-th step

R = ρ^{2^{N}} ,

where ρ is the separation radius of the original roots of f before applying Graeffe, we get that

\frac{2^{d+3}d}{R(i_{2} − i_{1})} \frac{|ζ_{d}|}{|ζ_{1}|} = 2^{d+3} \frac{d}{i_{2} − i_{1}} \frac{|ζ_{d}|}{|ζ_{1}|}ρ^{−2^{N}} ⟶ 0 , as N ⟶ ∞ .

4.6 ‘Deterioration of Condition’ and Stability Properties

It is important to understand that we will never have to solve g(x) = (G^{N}f)(x). Therefore, the actual condition number of g does not matter at all. In order to determine the roots of f, we will be using the extra information provided by ġ, where (g, ġ) = T G^{N}(f, ḟ).

A valid source of concern is the propagation of rounding-off error. In the Tangent Graeffe algorithm, that error would typically double at each step (it actually doubles at each step ‘in the limit’).

However, Lemmas 5 and 6 guarantee that the truncation error decreases as ρ^{−2^{N}}. Hence, in order to obtain a truncation error smaller than a certain δ > 0, we need N ≈ c + log_{2} log_{2} δ^{−1}, where c = − log_{2} log_{2} ρ is a constant depending only on the original polynomial f.

In order to reduce the accumulated rounding-off error to the same order, one would need that

Cϵ_{m}2^{N} < δ

where ϵ_{m} is the ‘machine epsilon’ and C is a constant depending on f. Thus, we just need

ϵ_{m} < \frac{δ}{C} 2^{−N} ≈ \frac{δ}{C2^{c} log_{2} δ^{−1}}

What is a reasonable value for δ ? The strength of Renormalized Tangent Graeffe Iteration is its capacity to solve the ‘global’ problem: given a polynomial, approximate all its roots. Once a suitable approximation of each root is found, local iterative algorithms (such as Newton Iteration) cheaply provide better refinements of the roots. Such a two-step procedure entails a reasonable range of values for δ. Namely, δ should be smaller (but not much smaller) than the radius of quadratic convergence of Newton’s iteration.

This radius is of the order of the reciprocal condition number of the original polynomial. (See [3] Theorem 1 and Remark 1 p. 263 for a precise statement).

5 Numerical Results and Final Remarks

5.1 Numerical Results

A polynomial solver based on the algorithms above was implemented and tested under IEEE 754 double and double-extended arithmetic. The results below are intended to make a case in favor of the stability and the practical feasibility of Renormalized Tangent Graeffe Iteration.

The first set of tests was designed to measure the performance of our algorithm for large degree polynomials. The test polynomials are pseudo-random real (Table 1 and Figure 3) and complex (Table 2 and Figure 4) polynomials, under the U(2)-invariant probability measure [17, 34, 22]. Under this probability measure, random polynomials are well-conditioned on the average.

Figure 3: Real pseudo-random polynomials

The results were certified using alpha-theory [38, 33, 20]. The running time (certification excluded) was compared to the code of Jenkins and Traub [14, 13] for the values where this code succeeds.

Running time is measured in user-time seconds of a Pentium-133 computer running Linux and the gcc compiler.

In the second set of experiments, we tried to check the behavior of Renormalized Tangent Graeffe Iteration in the presence of very badly conditioned polynomials. The test polynomials are Wilkinson’ s ‘perfidious’ polynomials [41]

p_{d}(x) = (x − 1)(x − 2) ⋯ (x − d)

and Chebyshev polynomials (Table 3).

T_{d} = ∏_{0≤m<d} (x − cos(\frac{π}{2d} + \frac{π}{d} m))

The error of the solutions of the perfidious polynomials is measured as max |ζ − round ζ|, ζ a ‘solution’ found by the program. Similarly, the error in Chebyshev polynomials is measured as max |m− round m| where m = \frac{d arccos ζ− \frac{π}{2}}{π} and ζ a ‘solution’ found by the program. Again, those results are compared to the ones provided by the software by Jenkins and Traub.

5.2 Further Practical Remarks

• Graeffe process (and hence our algorithm) is known to be parallelizable [12, 31].

• The algorithm presented here needs only O(d) memory storage; therefore, all intermediate computations for a reasonable degree d fit into the ‘cache’ memory of modern computers.

Running time (s)

Degree  Algorithm  Seed
0  1  2  3  4  5  6  7  8  9
50  Graeffe  0.08  0.08  0.08  0.07  0.09  0.08  0.08  0.08  0.08  0.07
J-T  0.03  0.02  0.04  0.02  0.05  0.05  0.07  0.04  0.04  0.03
100  Graeffe  0.20  0.20  0.20  0.21  0.21  0.21  0.19  0.20  0.19  0.20
J-T  0.12  0.11  0.12  0.09  0.12  0.12  0.14  0.09  0.10  0.10
150  Graeffe  0.36  0.36  0.36  0.36  0.37  0.35  0.36  0.37  0.35  0.36
J-T  0.26  0.26  0.28  0.24  0.23  0.24  0.27  0.20  0.26  0.24
200  Graeffe  0.56  0.54  0.56  0.56  0.55  0.55  0.56  0.55  0.56  0.55
J-T  0.53  0.42  0.51  0.35  0.43  0.38  0.54  0.40  0.47  0.36
250  Graeffe  0.77  0.78  0.77  0.78  0.78  0.78  0.78  0.78  0.77  0.76
300  Graeffe  1.02  1.02  1.01  1.02  1.02  1.00  1.02  1.01  1.01  1.03
350  Graeffe  1.29  1.29  1.29  1.29  1.30  1.29  1.29  1.29  1.28  1.29
400  Graeffe  1.59  1.59  1.57  1.58  1.58  1.60  1.57  1.59  1.59  1.58
450  Graeffe  1.90  1.90  1.89  1.91  1.89  1.89  1.90  1.90  1.91  1.91
500  Graeffe  2.32  2.27  2.31  2.30  2.31  2.31  2.30  2.29  2.32  2.34
550  Graeffe  2.57  2.58  2.58  2.58  2.59  2.59  2.58  2.59  2.57  2.59
600  Graeffe  2.96  2.98  2.97  2.96  2.95  2.97  2.96  2.97  2.98  2.96
650  Graeffe  3.38  3.39  3.37  3.36  3.37  3.34  3.38  3.37  3.37  3.37
700  Graeffe  3.78  3.78  3.78  3.80  3.79  3.77  3.78  3.80  3.78  3.77
750  Graeffe  4.21  4.22  4.21  4.21  4.19  4.20  4.19  4.22  4.20  4.21
800  Graeffe  4.66  4.66  4.65  4.65  4.64  4.64  4.63  4.65  4.65  4.65
850  Graeffe  5.18  5.20  5.18  5.17  5.23  5.18  5.19  5.20  5.19  5.17
900  Graeffe  5.61  5.61  5.58  5.59  5.60  5.57  5.56  5.60  5.60  5.60
950  Graeffe  6.16  6.16  6.15  6.16  6.13  6.14  6.16  6.18  6.14  6.18
1000  Graeffe  6.60  6.58  6.58  6.58  6.59  6.61  6.59  6.60  6.58  6.59

Table 1: Real pseudo-random polynomials Figure 4: Complex pseudo-random polynomials

Running time (s)

Degree  Algorithm  Seed
0  1  2  3  4  5  6  7  8  9
50  Graeffe  0.13  0.12  0.12  0.13  0.12  0.10  0.11  0.10  0.11  0.11
J-T  0.07  0.08  0.09  0.09  0.07  0.08  0.07  0.08  0.07  0.09
100  Graeffe  0.34  0.32  0.32  0.32  0.33  0.33  0.30  0.33  0.31  0.32
J-T  0.25  0.28  0.26  0.25  0.28  0.26  0.26  0.27  0.26  0.25
150  Graeffe  0.60  0.60  0.61  0.60  0.61  0.60  0.60  0.61  0.60  0.60
J-T  0.55  0.60  0.58  0.66  0.61  0.60  0.60  0.59  0.59  0.59
200  Graeffe  0.95  0.95  0.97  0.94  0.96  0.94  0.92  0.94  0.94  0.93
J-T  1.13  1.17  1.07  1.04  1.09  1.06  1.05  1.06  1.03  1.05
250  Graeffe  1.33  1.34  1.30  1.32  1.35  1.32  1.32  1.32  1.35  1.33
J-T  1.63  1.68  1.64  1.68  1.66  1.70  1.66  1.63  1.63  1.66
300  Graeffe  1.77  1.76  1.77  1.78  1.77  1.76  1.77  1.75  1.79  1.76
J-T  2.34  2.49  2.43  2.44  2.50  2.46  2.48  2.45  2.36  2.51
350  Graeffe  2.26  2.30  2.26  2.25  2.28  2.34  2.08  2.36  2.27  2.26
J-T  3.30  3.41  3.39  3.34  3.36  3.72  3.51  3.45  3.35  3.52
400  Graeffe  2.74  2.75  2.77  2.76  2.77  2.76  2.75  2.77  2.74  2.74
450  Graeffe  3.27  3.28  3.28  3.30  3.38  3.30  3.30  3.28  3.34  3.30
500  Graeffe  3.92  3.87  3.89  3.89  3.88  3.89  3.91  3.87  3.90  3.91
550  Graeffe  4.49  4.48  4.51  4.52  4.50  4.52  4.51  4.49  4.50  4.48
600  Graeffe  5.16  5.14  5.17  5.18  5.15  5.14  5.15  5.15  5.18  5.13
650  Graeffe  5.89  5.88  5.87  5.83  5.83  5.86  5.89  5.84  5.89  5.86
700  Graeffe  6.62  6.59  6.59  6.58  6.59  6.61  6.59  6.59  6.58  6.62
750  Graeffe  7.37  7.42  7.40  7.32  7.48  7.27  7.30  7.30  7.43  7.31
800  Graeffe  8.10  8.06  8.07  8.03  8.10  8.07  8.10  8.06  8.11  8.09
850  Graeffe  8.96  8.94  8.92  8.95  8.97  8.94  8.92  8.92  8.95  8.96
900  Graeffe  9.72  9.70  9.70  9.70  9.72  9.72  9.71  9.70  9.67  9.71
950  Graeffe  10.63  10.67  10.65  10.64  10.67  10.65  10.67  10.61  10.65  10.60
1000  Graeffe  11.49  11.47  11.54  11.46  11.50  11.49  11.54  11.50  11.51  11.50

Table 2: Complex pseudo-random polynomials

Perfidious polynomials

Degree  Algorithm  Error
10  Graeffe  5.123013 × 10−12
J-T  4.859935 × 10−11
15 Graeffe 3.968295 × 10−08
J-T  5.508868 × 10−09
20 Graeffe 1.780775 × 10−03
J-T  1.275754 × 10−04

Chebyshev polynomials

Degree  Algorithm  Error
10  Graeffe  8.790711 × 10−16
J-T  8.790711 × 10−16
15 Graeffe 2.169163 × 10−15
J-T  2.169163 × 10−15
20 Graeffe 1.903848 × 10−14
J-T  1.278977 × 10−13
25 Graeffe 1.266375 × 10−11
J-T  1.663025 × 10−11
30 Graeffe 5.511325 × 10−11
J-T  3.301608 × 10−10
35 Graeffe 5.708941 × 10−09
J-T  3.840000 × 10−08

Table 3: Perfidious and Chebyshev polynomials • Use of higher precision may be required to handle very badly conditioned poly-

nomials; as a matter of fact, polynomials of that sort are only meaningful if their coefficients are known with sufficiently high accuracy. For example, when they are obtained by symbolic manipulation.

Acknowledgements

JPZ’s work was supported in part by CNPq, through grant MA 521.329/94-9, and by PRONEX grant 76.97.1008-00. GM’s work was supported by CNPq through grant MA 520.423/96-8 and FAPERJ E-26/170.027/98.

References

[1] V. Arnold, Méthodes Mathématiques de la Mécanique Classique, Mir, Moscow, 1976.

[2] D. Bini and V. Y. Pan, Graeffe’s, Chebyshev-like, and Cardinal’s processes for splitting a polynomial into factors, J. Complexity, 12 (1996), pp. 492–511. Special issue for the Foundations of Computational Mathematics Conference (Rio de Janeiro, 1997).

[3] L. Blum, F. Cucker, M. Shub, and S. Smale, Complexity and Real Computation. Springer, 1998.

[4] Mr S. Brodetsky and G. Smeal, On Graeffe’s method for complex roots of algebraic equations, Proceedings of the Cambridge Philosophical Society XXII Part II (1924), pp. 83–87

[5] J.-P. Dedieu, À propos de la méthode de Dandelin-Graeffe, C. R. Acad. Sci. Paris Sér. I Math., 309 (1989), pp. 1019–1022.

[6] J. P. Dedieu, X. Gourdon, and J. C. Yakoubsohn, Computing the distance from a point to an algebraic hypersurface, in The mathematics of numerical analysis (Park City, UT, 1995), vol. 32 of Lectures in Appl. Math., Amer. Math. Soc., Providence, RI, 1996, pp. 285–293.

[7] J. W. Demmel, Applied numerical linear algebra, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1997.

[8] A. A. Grau, On the reduction of number range in the use of the Graeffe process, J. Assoc. Comput. Mach., 10 (1963), pp. 538–544.

[9] P. Henrici, Applied and computational complex analysis. Vol. 3, Pure and Applied Mathematics, John Wiley & Sons Inc., New York, 1986. Discrete Fourier analysis—Cauchy integrals—construction of conformal maps—univalent functions, A Wiley-Interscience Publication.

[10] N. J. Higham, Accuracy and stability of numerical algorithms, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1996.

[11] A. S. Householder, Dandelin, Lobačevskiĭ, or Graeffe?, Amer. Math. Monthly, 66 (1959), pp. 464–466.

[12] P. Jana and B. Sinha, Fast parallel algorithms for Graeffe’ s root squaring, Comput. Math. Appl., 35 (1998), pp. 71–80.

[13] M. A. Jenkins, Algorithm 493 – Zeros of a real polynomial [C2], ACM Transactions on Mathematical Software, 1 (1975), pp. 178–189.

[14] M. A. Jenkins and J. F. Traub, Algorithm 419 – Zeros of a complex polynomial [C2], Comm. of the ACM, 15 (1972), pp. 97–99.

[15] L. V. Kantorovich and G. P. Akilov, Functional analysis, Pergamon Press, Oxford, second ed., 1982. Translated from the Russian by Howard L. Silcock.

[16] P. Kirrinis, Partial fraction decomposition in ℂ(z) and simultaneous newton iteration for factorization in ℂ[z]. Preprint, Bonn, (1995).

[17] E. Kostlan, Random polynomials and the statistical fundamental theorem of algebra, 1987. Preprint, MSRI, 1987.

[18] S. Lang, Differentiable Manifolds, Addison-Wesley Pub. Co., Reading, Mass., 1972.

[19] R. S. MacKay, Renormalisation in area-preserving maps, vol. 6 of Advanced Series in Nonlinear Dynamics, World Scientific Publishing Co. Inc., River Edge, NJ, 1993.

[20] G. Malajovich, On generalized Newton algorithms: quadratic convergence, path-following and error analysis, Theoret. Comput. Sci., 133 (1994), pp. 65– 84. Selected papers of the Workshop on Continuous Algorithms and Complexity (Barcelona, 1993).

[21] G. Malajovich and J. P. Zubelli, A fast and stable algorithm for splitting polynomials, Comput. Math. Appl., 33 (1997), pp. 1–23.

[22] G. Malajovich and J. P. Zubelli, On the geometry of Graeffe iteration, 1997. Informes de Matemática Série B-118, IMPA.

[23] C. T. McMullen, Complex dynamics and renormalization, vol. 135 of Annals of Mathematics Studies, Princeton University Press, Princeton, NJ, 1994.

[24] C. A. Neff and J. H. Reif, An efficient algorithm for the complex roots problem, J. Complexity, 12 (1996), pp. 81–115.

[25] A. Ostrowski, Recherches sur la méthode de Graeffe et les zéros des polynomes et des séries de Laurent, Acta Math., 72 (1940), pp. 99–155.

[26] , Recherches sur la méthode de Graeffe et les zéros des polynomes et des séries de Laurent. Chapitres III et IV, Acta Math., 72 (1940), pp. 157–257.

[27] V. Y. Pan, Optimal and nearly optimal algorithms for approximating polynomial zeros, Comput. Math. Appl., 31 (1996), pp. 97–138.

[28] V. Y. Pan, Solving a polynomial equation: some history and recent progress, SIAM Rev., 39 (1997), pp. 187–220.

[29] V. Y. Pan, M.-H. Kim, A. Sadikou, X. Huang, and A. Zheng, On isolation of real and nearly real zeros of a univariate polynomial and its splitting into factors, J. Complexity, 12 (1996), pp. 572–594. Special issue for the Foundations of Computational Mathematics Conference (Rio de Janeiro, 1997).

[30] F. P. Preparata and M. I. Shamos, Computational geometry, Texts and Monographs in Computer Science, Springer-Verlag, New York, 1985. An introduction.

[31] T. A. Rice and L. H. Jamieson, A highly parallel algorithm for root extraction, IEEE Trans. Comput., 38 (1989), pp. 443–449.

[32] A. Schönhage, Equation solving in terms of computational complexity, in Proceedings of the International Congress of Mathematicians, Vol. 1, 2 (Berkeley, Calif., 1986), Providence, RI, 1987, Amer. Math. Soc., pp. 131–153.

[33] M. Shub and S. Smale, Complexity of Bézout’s theorem. I. Geometric aspects, J. Amer. Math. Soc., 6 (1993), pp. 459–501.

[34] M. Shub and S. Smale, Complexity of Bezout’s theorem. II. Volumes and probabilities, in Computational algebraic geometry (Nice, 1992), vol. 109 of Progr. Math., Birkhäuser Boston, Boston, MA, 1993, pp. 267–285.

[35] M. Shub and S. Smale, Complexity of Bezout’s theorem. III. Condition number and packing, J. Complexity, 9 (1993), pp. 4–14. Festschrift for Joseph F. Traub, Part I.

[36] M. Shub and S. Smale, Complexity of Bezout’s theorem. V. Polynomial time, Theoret. Comput. Sci., 133 (1994), pp. 141–164. Selected papers of the Workshop on Continuous Algorithms and Complexity (Barcelona, 1993).

[37] M. Shub and S. Smale, Complexity of Bezout’s theorem. IV. Probability of success; extensions, SIAM J. Numer. Anal., 33 (1996), pp. 128–148.

[38] S. Smale, Newton’s method estimates from data at one point, in The merging of disciplines: new directions in pure, applied, and computational mathematics (Laramie, Wyo., 1985), Springer, New York, 1986, pp. 185–196.

[39] M. M. Vainberg, Variational methods for the study of nonlinear operators, Holden-Day Inc., San Francisco, Calif., 1964. With a chapter on Newton’s method by L. V. Kantorovich and G. P. Akilov. Translated and supplemented by Amiel Feinstein.

[40] J. H. Wilkinson, Rounding errors in algebraic processes, Prentice Hall, Englewood Cliffs NJ (1963).

[41] J. H. Wilkinson, The perfidious polynomial, in Studies in numerical analysis, vol. 24 of MAA Stud. Math., Math. Assoc. America, Washington, DC, 1984, pp. 1–28.