This is page 356 Printer: Opaque this

A1

Algorithms and Computational Concerns

Overview

Appendix 1 presents some of the basic algorithms for computing with polynomials and rational functions and discusses some of the complexity issues. Included is a discussion of root finding methods. It requires very little background and can essentially be read independently.

Algorithms and Computational Concerns

Algorithms and Computational Concerns 357

Polynomials lend themselves to computation perhaps more than any other object of analysis. Algorithms that involve special functions, differential equations, series, and the like usually must reduce at some point to a finite polynomial or rational approximation or truncation. This often allows analytic problems to be reduced to algebraic ones. This appendix will present, as a series of exercises, some of the principal algorithmic concerns. The reader is encouraged to experiment with the algorithms. With current technology this is most comfortably done in any of the large symbolic manipulation packages available. Code, actual or schematic, is not presented. Indeed, methods rather than algorithms are presented. Current “practical” best methods date quickly in this rapidly evolving area. It is also the authors’ belief that today’s theoretical curiosities may be vital for tomorrow’s algorithms as larger instances are calculated on faster machines. Historically we have already seen this happen repeatedly with algorithms such as the fast Fourier transform algorithm.

One of the most interesting lessons to be learned from the last few decades of revitalized interest in computational mathematics is that many of the most familiar mathematical algorithms, such as how to multiply large numbers, were very poorly understood, and indeed, still are incompletely analyzed. Very many of the familiar processes of mathematics, such as multiplication of large numbers or computation of determinants, can be computed far more expeditiously than allowed by the usual “school” algorithms. See, for example, Aho, Hopcroft, and Ullman [74]; Bini and Pan [92]; Borodin and Munro [75]; Borwein and Borwein [87]; Brent [74]; Knuth [81]; Pan [92]; Smale [85]; and Wilf [86] for the complexity side of the following exercises.

E.1 Complexity and Recursion. We are concerned with measuring the size of an algorithm given an input of a certain length. Unfortunately, there are many different ways of measuring this. (One could, for example, use the length of the tapes of some well-defined instantiation of a Turing machine.) We will settle for less. The measure of input size will usually be chosen to be a natural one, so for polynomials of degree n, the measure will often be n. The complexity measure then depends a bit on the problem. For example, it might count the number of additions of coefficients (we do not distinguish subtraction from addition) and multiplications of coefficients required to evaluate the polynomial at a point. (So, for example, by Horner’s rule O(n) operations suffice.) Care is already required to count naturally. Note that we have not specified the size of the coefficients (this may or may not be reasonable) and so we could cheat on addition of coefficients by doing two additions as one addition of twice the length. (Since a + b and c + d can be decoded from (10^{m}a + c) + (10^{m}b + d), where m is larger than the number of digits in any of a, b, c, or d.) It is more reasonable in this context to fix a precision (or to think of working to infinite precision or over some polynomial ring). Our cases are fairly simple, and the measures should be clear in context.

We adopt the following notations:

f(n) = O(g(n)) means lim sup_{n→∞} \frac{f(n)}{g(n)} < ∞

and

f(n) = Ω(g(n)) means lim sup_{n→∞} \frac{g(n)}{f(n)} < ∞ .

Parts of these exercises are reprinted from Borwein and Borwein [87], with permission from Wiley.

358 A1. Algorithms and Computational Concerns

So the first measure gives an upper bound, while the second gives a lower bound.

Many algorithms are analyzed recursively. For example, addition of two polynomials of degree at most 2n reduces to two additions of polynomials of degree at most n, plus perhaps an “overhead” for reassembling the pieces. In other words, for the complexity of addition, we have

A(2n) ≤ 2A(n) + c ,

from which one can deduce that

A(n) = O(n) .

This general recursive strategy of breaking a problem in half is often called “divide and conquer.”

We introduce the following functions. Here n is the maximum degree of the polynomials p and q. Additions, multiplications, and so on are performed in the underlying field of coefficients (in our case ℂ or ℝ) and are all performed to some predetermined fixed precision (possibly infinite).

A(n) := the maximum number of ± , × , ÷ , to compute p ± q;
M(n) := the maximum number of ± , × , ÷ , to compute pq;
e(n) := the maximum number of ± , × , ÷ , to evaluate

p(α) for an arbitrary fixed α ∈ ℂ;

E(n) := the maximum number of ± , × , ÷ , to evaluate

p(α_{1}), …, p(α_{n}) for arbitrary fixed α_{1}, …, α_{n} ∈ ℂ .

These are the complexity functions for polynomial addition, polynomial multiplication, polynomial evaluation at a single point, and polynomial evaluation at n points, respectively. The input for the computation is the coefficients (and the evaluation points for e(n) and E(n)). So the input may be considered to be in ℂ^{n+1} (or more generally an (n + 1)-dimensional vector space over an infinite field). In the first two cases the output is the sequence of coefficients. In the last two cases, respectively, the output is the evaluation and the sequence of evaluations.

a] Show that usual algorithms give

A(n) ≤ 2n + 2 = O(n) ,
M(n) = O(n^{2}) ,
e(n) = O(n)  (Horner’s rule) ,
E(n) = O(n^{2}) .

E.2 and E.3 of this appendix provide better upper bounds for the last three functions above. b] Show that

A(n) ≥ n + 1 ,
M(n) ≥ n + 1 ,
e(n) ≥ n + 1 ,
E(n) ≥ n + 1 .

Hint: In all cases this is a uniqueness argument. At least one operation must be performed for each coefficient, otherwise the algorithm will not distinguish various different sequences of input. ⊓ ⊔

c] Some Recursive Bounds. Let a, b > 0 and c, d > 1. Suppose that f is monotone on (0, ∞).

If f(n) ≤ af(n/c) + bn, then
f(n) = O(n)  if a < c ,
f(n) = O(n log n) if a = c ,
f(n) = O(n^{log_{c} a}) if a > c .
If f(n) ≤ df(n/d) + bn(log n)^{c−1}, then
f(n) = O(n(log n)^{c}) .

Hint: Analyze the equality case. Then establish the general principle that the equality solution is the maximal solution. ⊓ ⊔

E.2 The (Finite) Fast Fourier Transform (FFT). This is undoubtedly one of the most widely used algorithms. It has, in its various forms, tremendous practical and theoretical applications.

Let w be a primitive (n + 1)th root of unity in either ℂ or a finite field 𝔽_{m}, that is, w^{n+1} = 1 and w^{k} ≠ 1 for k = 1, 2, …, n. In the complex case we may take w := e^{2πi/(n+1)}. Consider the following two problems.

Interpolation Problem. Given n + 1 numbers, α_{0}, α_{1}, …, α_{n}, find the coefficients of the unique polynomial

p(z) := a_{0} + a_{1}z + ⋯ + a_{n}z^{n}

of degree n that satisfies

p(w^{k}) = α_{k} , k = 0, 1, …, n .

Evaluation Problem. Given the coefficients of a polynomial p_{n} of degree at most n, calculate the n + 1 values

p(w^{k}) , k = 0, 1, …, n .

These are the two directions of the (finite) Fourier transform. The classical approaches to either part of the Fourier transform problem have complexity at least cn^{2}. This is the complexity, for example, of evaluating p_{n} at n + 1 points using Horner’s rule. Both directions can, however, be solved with complexity O(n log n).

a] If n + 1 = 2^{m}, then both the interpolation and the evaluation problem have complexity O(n log n). (Here we are counting the number of additions and multiplications in the underlying coefficient field, which for most of our purposes is ℂ.)

Outline. We treat the evaluation first. Suppose

p(x) := a_{0} + a_{1}x + ⋯ + a_{n}x^{n} .

Let

q(x^{2}) := a_{0} + a_{2}x^{2} + a_{4}x^{4} + ⋯ + a_{n−1}x^{n−1}

and

x r(x^{2}) := x(a_{1} + a_{3}x^{2} + ⋯ + a_{n}x^{n−1}) .

Then, with y := x^{2},

p(x) = x r(y) + q(y) ,

where r and q are both polynomials of degree at most 2^{m−1} − 1. The observation that makes the proof work is that for w an (n + 1)th root of unity,

(w^{k})^{2} = (w^{(n+1)/2 + k})^{2} .

Hence, evaluating p(x) at the n + 1 roots of unity reduces to evaluating r and q each at the \frac{1}{2}(n + 1) points (w^{2})^{1}, (w^{2})^{2}, …, (w^{2})^{(n+1)/2} and amalgamating the results. Observe that w^{2} is a primitive (2^{m−1})th root of unity, so we can iterate this process. Let F (2^{m}) be the number of additions and multiplications required to evaluate a polynomial of degree at most 2^{m} − 1 at the 2^{m} points w^{k}, k = 1, 2, …, 2^{m}, where w is a primitive (2^{m})th root of unity. Then

F (1) = 1  and  F (2^{m}) = 2F (2^{m−1}) + 2 · 2^{m} ,  m = 1, 2, ….

The second term comes from the single addition and multiplication required to calculate each p(w^{k}) from r(w^{2k}) and q(w^{2k}). This recursion solves as

F (2^{m}) = 2^{m+1}m ,

and the bound for the evaluation problem is established.

The interpolation problem is equivalent to evaluation. This can be seen as follows. Let w be a primitive (n + 1)th root of unity, and let

W := ( \begin{matrix}1 & 1 & 1 & … & 1 \\ 1 & w & w^{2} & … & w^{n} \\ 1 & w^{2} & w^{4} & … & w^{2n} \\ ⋮ & ⋮ & ⋮ & ... & ⋮ \\ 1 & w^{n} & w^{2n} & … & w^{n^{2}} \end{matrix} ) .

Then

W^{−1} = \frac{1}{n + 1} ( \begin{matrix}1 & 1 & 1 & … & 1 \\ 1 & w^{−1} & w^{−2} & … & w^{−n} \\ 1 & w^{−2} & w^{−4} & … & w^{−2n} \\ ⋮ & ⋮ & ⋮ & ... & ⋮ \\ 1 & w^{−n} & w^{−2n} & … & w^{−n^{2}} \end{matrix} )

and w^{−1} is also a primitive (n+1)th root of unity. The interpolation problem can be formulated as follows. Find (a_{0}, …, a_{n}) so that

W (a_{0}, a_{1}, …, a_{n}) = (α_{0}, α_{1}, …, α_{n}).

However, this can be solved by

W^{−1}(α_{0}, α_{1}, …, α_{n}) = (a_{0}, a_{1}, …, a_{n}) ,

which is exactly the evaluation problem. ⊓ ⊔

See Borodin and Munro [75] and Borwein and Borwein [87]. Versions of FFT exist in a plethora of shapes and sizes. We have just exposed the tip of the iceberg.

As an application we construct a fast polynomial multiplication.

b] Fast Polynomial Multiplication. Suppose the polynomials p, q of degree at most n−1 are given. Compute the coefficients of the product pq as follows:

b1] Use an FFT to evaluate p and q at the primitive (2n)th roots of unity w^{1}, w^{2}, …, w^{2n}.

b2] Form

p(w^{k})q(w^{k}) , k = 1, 2, …, 2n .

b3] Find the coefficients of the product pq by solving the interpolation problem once again by using the FFT. Show that this algorithm requires

O(n log n)

additions, multiplications, and divisions (of complex numbers) and therefore M(n) = O(n log n).

This is the best-known upper bound on the serial complexity of polynomial multiplication. The only known lower bound is the trivial one O(n). In parallel (on a PRAM) polynomial multiplication can be done in O(log n) time on O(N) processors; see Pan [92]. The same bounds apply for the FFT in a].

E.3 Other Elementary Operations.

a] Fast Polynomial Division. For polynomials p of degree n and q of degree m ≤ n, it is possible to find polynomials u and r with deg(r) < deg(q) such that

p(x) = u(x)q(x) + r(x)

in O(n log n) additions and multiplications (of complex numbers).

Outline. Simplify by observing that it suffices to calculate u since r may then be computed by E.2 b]. If we replace x by 1/x then

\frac{p(x^{−1})}{q(x^{−1})} = u(x^{−1}) + \frac{r(x^{−1})}{q(x^{−1})}

and so, for some h ≥ 1,

\frac{p^{∗}(x)}{q^{∗}(x)} = u^{∗}(x) + x^{n−m+h} \frac{r^{∗}(x)}{q^{∗}(x)} ,  where  v^{∗}(x) := x^{deg(v)}v(x^{−1}) .

To calculate u^{∗} (and hence to calculate u) it suffices to calculate the first n−m (= deg(u)) Taylor coefficients of 1/q^{∗}. This can be done by Newton’s method (see the next exercise) as follows: Suppose deg s_{i} = j − 1 and

\frac{1}{q^{∗}(x)} − s_{i}(x) = O(x^{j}) .

Establish that

\frac{1}{q^{∗}(x)} − [2s_{i}(x) − s_{i}^{2}(x)q^{∗}(x)] = \frac{1}{q^{∗}(x)} [1 − s_{i}(x)q^{∗}(x)]^{2} = O(x^{2j}) .

Note that we may assume q^{∗}(0) ≠ 0. Now the computation of

s_{i+1} := 2s_{i} − s_{i}^{2}q^{∗}

can be performed by using an FFT-based polynomial multiplication and it needs only be performed by using the first 2j − 1 coefficients of q^{∗} and s_{i}. By starting with an appropriate first estimate of s_{0} (say, s_{0}(x) := 1/q^{∗}(0)), and proceeding inductively as above (doubling the number of coefficients used at each stage) we can show that the required number of terms of the expansion can be calculated in O(n log n) additions and multiplications (of complex numbers).

b] Fast Reversion of Power Series. Let

f(x) = ∑_{k=0}^{∞} a_{k}x^{k} , a_{0} ≠ 0

be a formal power series with given coefficients. Show that the first n coefficients of the formal Taylor expansion of 1/f can be computed in O(n log n) additions, multiplications, and divisions (of complex numbers).

c] Fast Polynomial Expansion. Given α_{1}, α_{2}, …, α_{n} show that the coefficients of ∏_{i=1}^{n}(x − α_{i}) can all be calculated in O(n(log n)^{2}) additions, multiplications, and divisions (of complex numbers).

Hint: Proceed recursively by dividing the problem into two parts of roughly half the size. Recombine the pieces using part b] of the previous exercise.

⊓ ⊔

d] Fast Polynomial Expansion at Arbitrary Points. Given a polynomial p of degree at most n, and n + 1 distinct points x_{0}, x_{1}, …, x_{n}, show that p(x_{0}), p(x_{1}), …, p(x_{n}) can all be evaluated in O(n(log n)^{2}) multiplications and additions.

Hint: Let

q_{1}(x) := ∏_{i=0}^{⌊n/2⌋−1}(x − x_{i})

and let r_{1} be the remainder on dividing p by q_{1}. Note that r_{1}(x_{i}) = p(x_{i}) for each i < n/2. Similarly, use

q_{2}(x) := ∏_{⌊n/2⌋}^{n} (x − x_{i}).

Thus two divisions reduce the problem to two problems of half the size. ⊓ ⊔

e] Extend d] to rational functions.

f] Evaluation of x^{n}. The S-and-X binary method for calculating x^{n} is the following algorithm. Suppose n has binary representation δ_{0}δ_{1}δ_{2} … δ_{k} with δ_{0} = 1. Given symbols S and X, define

S_{i} := \{ \begin{matrix}SX & if δ_{i} = 1 \\ S & if δ_{i} = 0 \end{matrix}

and construct the rule

S_{1}S_{2} ⋯ S_{k} .

Now let S be the operation of squaring, and let X be the operation of multiplying by x. Let S_{1}S_{2} ⋯ S_{k} operate from left to right beginning with x. For example, for n = 27,

δ_{0}δ_{1}δ_{2}δ_{3}δ_{4} = 11011

and

S_{1}S_{2}S_{3}S_{4} = (SX)(S)(SX)(SX) .

The sequence of calculations of x^{27} is then

x → x^{2} → x^{3} → x^{6} → x^{12} → x^{13} → x^{26} → x^{27}.

f1] Prove that the above method computes x^{n} and observe that it only requires storing x, n, and one partial product.

f2] Show that the number of multiplications in the S-and-X method is less than 2⌊log_{2} n⌋.

f3] Show that the S-and-X method is optimal for computation of x^{2m} (considering only multiplications).

f4] Show that the S-and-X method is not optimal for computing x^{15}.

An extended discussion of this interesting and old problem is presented in Knuth [81].

Much further material on the complexity of polynomial operations and complexity generally may be found in Aho, Hopcroft, and Ulman [74]; Borodin and Munro [75]; Pan [92]; and Wilf [86].

E.4 Newton’s Method. One of the very useful algorithms, both in theory and practice, for zero finding is Newton’s method.

a] Suppose f is analytic in a (complex) neighborhood of z_{0}, and suppose f(z_{0}) = 0 and f′(z_{0}) ≠ 0. Show that the iteration

x_{n+1} := x_{n} − \frac{f(x_{n})}{f′(x_{n})}

converges locally uniformly quadratically, that is, with a constant c independent of n,

|x_{n+1} − z_{0}| < c |x_{n} − z_{0}|^{2}

for initial values x_{0} in some neighborhood of z_{0}. As before, we call this locally quadratic convergence.

Hint: Note that

f(x_{n}) = f(z_{0}) + (x_{n} − z_{0})f′(z_{0}) + O((x_{n} − z_{0})^{2}) ,

which implies

x_{n+1} − z_{0} = (x_{n} − z_{0}) [\frac{f′(x_{n}) − f′(z_{0})}{f′(x_{n})}] + O((x_{n} − z_{0})^{2}) .

⊓ ⊔

b] Show that the iteration

x_{n+1} := x_{n} + \frac{f(x_{n})}{f′(x_{n})}

converges locally quadratically to the simple poles of a meromorphic function f. Note that this is Newton’s method with the sign changed.

c] The iteration

x_{n+1} := x_{n} − \frac{f(x_{n})f′(x_{n})}{f′(x_{n})^{2} − f(x_{n})f′′(x_{n})}

converges locally quadratically to a zero of an analytic f independent of its multiplicity; see Henrici [74].

d] Let g := f^{−1}. The iteration

x_{n+1} := x_{n} + (n + 1) \frac{g^{(n)}(x_{n})}{g^{(n+1)}(x_{n})}

converges locally uniformly to a zero of an analytic function f with (n + 2)th order. Newton’s method is n := 0; Halley’s method is n := 1 (see Householder [70]).

Newton’s method and its variants work tremendously well provided that a good starting value can be found. This is a problem. On a real interval a bisection method can be used initially to localize the zeros. In the plane, life is more complicated as is seen in the next exercise. Another drawback to Newton’s method is the need to compute the derivative. Of course, this is not a problem for polynomials, but in a general setting it is usually replaced by an approximation such as

\frac{f(x_{n}) − f(x_{n−1})}{x_{n} − x_{n−1}}

(which yields the so-called secant method).

e] Consider Newton’s method for computing \sqrt{x} starting with x_{0} := 1. This gives

x_{n+1} := \frac{1}{2} (x_{n} + \frac{x}{x_{n}} ) .

Show that r_{n}(x) := x_{n+1} is a rational function in x with numerator of degree 2^{n} and denominator of degree 2^{n} − 1. Show that r_{n}(x) − \sqrt{x} has a zero of order 2^{n+1} at 1. So r_{n} is in fact the (2^{n}, 2^{n} − 1) Padé approximant to \sqrt{x}. (This implies, though not obviously, that r_{n} has all real negative roots and poles.)

An attractive feature of Newton’s method is that it is “self-correcting.” So, for example, to compute a root to a large precision, one can start at a small precision and double the precision at each iteration. This is a substantial savings both practically and theoretically. The same feature applies to Newton’s method solutions over formal power series, as in E.3 a]. This allows for doubling the number of terms used at each stage. Much additional material on Newton’s method is available in Borwein and Borwein [87], Henrici [74], Householder [70], and Traub [82].

E.5 Newton’s Method in Many Variables.

a] Let f : ℂ^{n} → ℂ^{n}, and suppose f has Jacobian

J := |\frac{∂f_{i}}{∂z_{j}} | ,

where f := (f_{1}, f_{2}, …, f_{n}) with f_{i} : ℂ^{n} → ℂ. Let

x := (x_{1}, x_{2}, …, x_{n}) ∈ ℂ^{n} .

The function J(x) is the Jacobian evaluated at x. Then Newton’s method becomes

x_{k+1} = x_{k} − s_{k}

where s_{k} solves J(x_{k})s_{k} = f(x_{k}). This iteration converges locally uniformly quadratically to a zero z_{0} of f, that is, with a constant c independent of n,

|x_{n+1} − z_{0}| < c |x_{n} − z_{0}|^{2}

for initial values x_{0} in some neighborhood of z_{0}, provided in a neighborhood of z_{0}, f is continuously differentiable, J^{−1} exists and is bounded in norm, and J satisfies a Lipschitz condition. We call this locally quadratic convergence. For a polynomial f, we require only that J^{−1} exists in a neighborhood of the zero z_{0}. (For examples and detail, see Dennis and Schnabel [83].)

b] Finding All Zeros of a Polynomial. Let

p(z) := a_{0} + a_{1}z + ⋯ + a_{n}z^{n} , a_{n} := 1 .

Let f_{i}(x_{1}, …, x_{n}) be the ith coefficient of

g(x) := ∏_{i=1}^{n}(x − x_{i}) − ∑_{i=0}^{n} a_{i}x^{i}

and let

f(x) := (f_{1}, f_{2}, …, f_{n}) .

Then the iteration of part a] applied to f converges locally quadratically to z_{0} := (z_{1}, z_{2}, …, z_{n}), where z_{1}, z_{2}, …, z_{n} are the zeros of p, provided the zeros of p are distinct.

c] Another Approach to Finding All Zeros. Let x_{j}(k) denote the kth approximation to the jth root of p, where p is a polynomial of degree n with n distinct roots. Let

x_{j}(k + 1) = x_{j}(k) − \frac{p(x_{j}(k))}{∏_{i≠j}^{k}(x_{j}(k) − x_{i}(x))} , j = 1, 2, …, n .

Show that for sufficiently good choice of (x_{0}(1), …, x_{0}(k)), the above iteration converges locally quadratically to a sequence of the n distinct zeros of p.

Hint: This is really just the single variable version of Newton’s method for each root, where the derivative is approximated by the derivative of the kth estimate. ⊓ ⊔

d] Observe that the iteration of part c] fails to converge for p(x) := x^{n} −1, n > 2, if the starting values are all taken to be real.

In practice, the iteration of part c] works rather well for reasonably chosen starting values. One can use the techniques of the next exercise to localize the zeros first. With care, an algorithm can be given that computes a zero of a polynomial with an error < 2^{−b} in O(n log b log n) time and all zeros in O(n^{2} log b log n) time; see Pan [92]. On a parallel machine (PRAM) an algorithm requiring O(log^{3}(nb)) time and O(nb)^{O(1)} processors can be given for computing all zeros.

e] Modify the iteration of E.4 c] as given in the previous exercise to get a method that computes all roots even in the presence of repeated roots.

For further discussion, see Aberth [73], Durand [60], Kerner [66], and Werner [82].

E.6 Localizing Zeros.

a] Cauchy Indices. Let r be a real rational function with a real pole α. The Cauchy index of r at α is

1 if lim_{x→α−} r(x) = −∞ and lim_{x→α+} r(x) = ∞
−1 if lim_{x→α−} r(x) = ∞ and lim_{x→α+} r(x) = −∞
0 otherwise .

The Cauchy index of r on an interval [a, b] is the sum of the Cauchy indices of the poles of r in (a, b). (We demand that neither a nor b be poles of r.) We denote this by I_{a}^{b}(r).

b] The Euclidean Algorithm. Let p_{0} and p_{1} be nonzero polynomials. Define polynomials p_{0}, p_{1}, …, p_{m} and q_{1}, q_{2}, …, q_{m} (by the usual division algorithm) so that

p_{0}(z) = p_{1}(z)q_{1}(z) − p_{2}(z) ,  deg(p_{2}) < deg(p_{1}) ;
p_{1}(z) = p_{2}(z)q_{2}(z) − p_{3}(z) ,  deg(p_{3}) < deg(p_{2}) ;
⋮
p_{m−1}(z) = p_{m}(z)q_{m}(z) − 0 ;

the algorithm stops the first time that the remainder is zero. This is called the Euclidean algorithm. Show that p_{m} is the greatest common factor of p_{0} and p_{1}.

c] Let p_{0}, p_{1}, …, p_{m} be the polynomials generated by the Euclidean algorithm in part b]. Let P_{i} := p_{i}/p_{m}, i = 0, 1, …, m. Suppose p_{0}(a)p_{0}(b) ≠ 0. Show that if P_{k}(γ) = 0 for some k and a ≤ γ ≤ β, then

P_{1}(γ) ≠ 0 if k = 0

and

P_{k−1}(γ)P_{k+1}(γ) < 0 if 1 ≤ k ≤ m − 1 .

d] Show that for real polynomials p_{0} and p_{1} with p_{0}(a)p_{0}(b) ≠ 0,

I_{a}^{b}(p_{1}/p_{0}) = v(a) − v(b) , a < b ,

where v(α) is the number of sign changes in the sequence

(p_{0}(α), p_{1}(α), …, p_{m}(α))

and where the polynomials p_{i} are generated by the Euclidean algorithm as in part b]. (As before, by a sign change we mean that p_{i}(α)p_{i+k}(α) < 0 and p_{i+1}(α) = p_{i+2}(α) = ⋯ = p_{i+k−1}(α) = 0.)

Hint: Without loss of generality, we may assume that p_{m} = 1 (why?), so P_{k} = p_{k} for each k. First, note that v(x) may change magnitude only if p_{i}(x) = 0 for some i. By continuity of the polynomials p_{i}, and by part c], v(x) is constant on any subinterval of [a, b] that does not contain a zero of p_{0}.

We are now reduced to considering the behavior of v(x) at the zeros of p_{0}. Consider the various possibilities for the behavior of p_{1}/p_{0} at the zeros of p_{0} by considering the four possible changes of signs of p_{0} and p_{1} at the zeros of p_{0} and the effect this has on the Cauchy index and the increase and decrease of v(x). (Note that v(x) decreases by 1 if the Cauchy index is 1 while v(x) increases by 1 if the Cauchy index is −1.) ⊓ ⊔

e] Zeros on an Interval. Suppose p is a real polynomial and p(a)p(b) ≠ 0. Then I_{a}^{b}(p′/p) equals the number of distinct zeros of p in [a, b]. This also equals v(a) − v(b), where v(x) is as in d] with p_{0} := p and p_{1} := p′.

Hint: Use

\frac{p′(z)}{p(z)} = ∑ \frac{m_{k}}{z − ζ_{k}}

and note that each distinct pole in (a, b) contributes +1 to the Cauchy index. ⊓ ⊔

f] Show that the number of zeros of a real polynomial on the real line equals

lim_{a→∞}(v(−a) − v(a)) ,

where v(x) is computed as in part e]. (That is, v(x) is the number of sign changes in the sequence (p_{0}, p_{1}, …, p_{m}) generated by the Euclidean algorithm with p_{0} := p and p_{1} := p′.)

g] The Number of Zeros in H := {z ∈ ℂ : Im(z) > 0}. Suppose that the monic polynomial p ∈ P_{n}^{c} has (exactly) k real zeros counting multiplicities. Write

p(x) = r(x) + i s(x) , r, s ∈ P_{n} .

Then the number of zeros of p in H is

\frac{1}{2} [n − k − I_{−∞}^{∞}(s/r)] ,

and this can be computed as in part d] by using

I_{−∞}^{∞}(s/r) = lim_{α→∞} I_{−α}^{α}(s/r) .

Hint: First consider the case where k = 0. Consider p on a counterclockwise semicircular contour with base [−α, α] and radius α. Consider the argument of p as the contour is traversed. On the half-circle, for large α, the argument increases by nπ asymptotically; while on the axis the change is −πI_{−α}^{α}(s/r), from which the result follows. ⊓ ⊔

h] Budan-Fourier Theorem. Let p ∈ P_{n}. Let V (x) be the number of sign changes in the sequence

(p(x), p′(x), …, p^{(n)}(x)) .

Then the number of zeros of p in the interval [a, b], counting multiplicities, is V (a) − V (b) − 2m for some nonnegative integer m.

We have followed Henrici [74] in this discussion.

E.7 Zeros in a Disk.

a] The transform \frac{z + i}{1 + iz} maps the unit circle to the real axis and maps

the open unit disk D to the upper half-plane {z ∈ ℂ : Im(z) > 0}. So the algorithms of the previous exercise apply after transformation. b] Zero Counting by Winding Number. Let ω_{m} := e^{2πi/m}, where m is a positive integer. If p ∈ P_{n}^{c} then

z_{p} := \frac{1}{2πi} ∫_{∂D} \frac{p′(z)}{p(z)} dz ∼ ∑_{k=1}^{m} ω_{m}^{k} \frac{p′(ω_{m}^{k})}{p(ω_{m}^{k})}

counts the number of zeros of p in the open unit disk D (assuming no zeros on the boundary ∂D). More precisely, we have

z_{p} = lim_{m→∞} ∑_{k=1}^{m} ω_{m}^{k} \frac{p′(ω_{m}^{k})}{p(ω_{m}^{k})} .

(Note that this lends itself to rapid evaluation by FFT methods.)

c] Show that

\frac{1}{2πi} ∫_{∂D} \frac{1}{z − α} = ∑_{k=1}^{m} \frac{ω_{m}^{k}}{ω_{m}^{k} − α} + ϵ_{m}(α) ,

where

ϵ_{m}(α) ≤ \frac{|α|^{m−1}}{1 − |α|}  if |α| < 1

and

ϵ_{m}(α) ≤ \frac{|α|^{−m−1}}{1 − |α|^{−1}}  if |α| > 1 .

Hint: Write

\frac{1}{z − α} = −\frac{1}{α} ∑_{k=0}^{∞} \frac{z^{k}}{α^{k}}

and use the fact that

0 = ∫_{∂D} p(z) dz = ∑_{k=1}^{m} ω_{m}^{k}p(ω_{m}^{k})

for every p ∈ P_{m−1}^{c}. ⊓ ⊔

d] If p ∈ P_{n}^{c} has no zeros in the annulus ζ ≤ z ≤ 1/ζ and if m in part b] is greater than

⌊\frac{log \frac{1−ζ}{n}}{log ζ} ⌋ + 1 ,

then the error in estimating z_{p} by the sum is less than 1. So this provides an algorithm.

E.8 Computing General Chebyshev Polynomials. Given a Chebyshev system M := (f_{0}, …, f_{n}) of C^{1} functions on, say, [0, 1] how does one compute the associated Chebyshev polynomial T_{n}? That is, how does one find the unique equioscillating form T := ∑_{i=0}^{n} a_{i}f_{i} of Section 3.3?

a] The Remez Algorithm

Step 1. Choose

x_{0} := (0 =: x_{0}^{(0)} < x_{1}^{(0)} < ⋯ < x_{n}^{(0)} := 1)

and find P_{0} ∈ span M such that

P_{0}(x_{i}^{(0)}) = (−1)^{i} , i = 0, 1, …, n .

Step 2. Inductively set x_{m} := (0 =: x_{0}^{(m)} < x_{1}^{(m)} < ⋯ < x_{n}^{(m)} := 1), where

P′_{m−1}(x_{i}^{(m)}) = 0 , i = 1, 2, …, n − 1 .

(That is, find the extrema of P_{m−1}.)

Step 3. Find P_{m} ∈ span M with

P_{m}(x_{i}^{(m)}) = (−1)^{i} , i = 0, 1, …, n .

Then, provided the initial estimate x_{0} is sufficiently good, P_{m} → T_{n} quadratically (see Veidinger [60]).

This is reasonably easy to code. It involves solving an interpolation problem in Steps 1 and 3. The zero finding at Step 2 can be done quite easily since one can find very good starting values for Newton’s method, namely (x_{i}^{(m−1)})_{i=0}^{n}.

b] This algorithm modifies to solve the best approximation problem

min_{p∈span M} ∥ωp − f∥_{[a,b]}

for ω, f ∈ C[a, b], where ω is positive on [a, b]. One uses the Remez agorithm to find an equioscillating form

P (x) = f(x) − ω(x) ∑_{i=0}^{n} c_{i}f_{i}(x)

at n + 2 points. To do this, one solves the system

f(x_{k}) − ω(x_{k}) ∑_{i=0}^{n} c_{i}f_{i}(x_{k}) + (−1)^{k}λ , k = 0, 1, …, n + 1

for both the c_{i} and λ. (This works reasonably well, provided that at each stage ∥P ∥_{[a,b]} occurs at one of the x_{k}. If not, an extra point must be inserted where the maximum norm occurs and one of the original points must be dropped. This is effected in such a way as to maintain the alternations in sign of the error.) For details see Cody, Fraser, and Hart [68] and Veidinger