Lifting/descending processes for polynomial zeros

Bernard Mourrain

SAGA, INRIA, 06902, Sophia Antipolis, Cedex, France internet: mourrain@sophia.inria.fr

Victor Y. Pan^{∗} Mathematics and Computer Science Department Lehman College, CUNY, Bronx, NY 10468 internet: vpan@alpha.lehman.cuny.edu

Abstract

The recently proposed Chebyshev-like lifting map for the zeros of a univariate polynomial was motivated by its applications to splitting a univariate polynomial p(x) numerically into factors, which is a major step of some most effective algorithms for approximating polynomial zeros. We complement the Chebyshev-like lifting process by a descending process, decrease the estimated computational cost of performing the algorithm, demonstrate its correlation to Graeffe’s lifting/descending process and generalize lifting from Graeffe’s and Chebyshev-like maps to any fixed rational map of the zeros of the input polynomial.

Key words: polynomial zeros, numerical splitting a polynomial into factors, lifting polynomial zeros.

1991 Math. Subject Classification: 26C10, 30C15, 65H05, 68Q40.

1 Introduction and background.

Graeffe’s map squares the zeros of a fixed univariate polynomial

p(z) = ∑_{i=0}^{n} p_{i}z^{i} = p_{n} ∏_{j=1}^{n}(z − z_{j}), p_{n} ≠ 0,  (1.1)

∗Supported by NSF Grants CCR 9625344 and 9732206 and PSC-CUNY Award 669363

that is , transforms p(x) into the polynomial p̃(z) = p̃_{n} ∏_{j=1}(x − z_{j}^{2}). The transformation is simple (it essentially amounts to polynomial multiplication and costs O(n log n) operations) and is used in various algorithms for polynomial rootfinding. Hereafter, we will write ”ops” as abbreviation for ”arithmetic operation”.

In [BP96] this map was extended to the transformation of p(x) into the polynomial p̂_{n} ∏_{j}(x−(z_{j} +z_{j}^{−1})/2) with the zeros (z_{j} +z_{j}^{−1})/2, and some applications were shown to splitting p(x) numerically into factors, which is the basic step of some of the most efficient known polynomial rootfinders [P95], [P96]. In [BP96], this transformation was called Chebyshev-like lifting and the estimate O(n log^{2} n) ops was shown for performing it.

In this paper we observe show that the actual cost is O(n log n) (cf. sections 2 and 3), complement Chebyshev-like lifting map by the descending map (section 2), demonstrate correlation between Chebyshev-like and Graeffe’s maps (section 3), and generalize Graeffe’s and Chebychev-like maps to computing a polynomial p(x) = p_{n} ∏_{j=1}^{n}(x − r(z_{j})) for any fixed rational function r(z) (section 4).

2 Chebyshev-like lifting/descending steps.

Chebyshev-like lifting algorithm was introduced in [BP96] as a means of the transition from the polynomial p(z) of (1.1) to the polynomial

p̂(z) = p_{n}p_{0} ∏_{j=0}^{n}(z − \frac{z_{j} + z_{j}^{−1}}{2} )  (2.1)

with the zeros \frac{z_{j}+z_{j}^{−1}}{2} , provided that p_{n}p_{0} ≠ 0. For the sake of completeness of our presentation, we will next recall this algorithm, based on the equations

p̂(x) = p(z)p(z^{−1}), x = \frac{z + z^{−1}}{2} ,  (2.2)

and we will refine its computational cost estimate, versus the pessimistic one claimed in [BP96].

Algorithm 2.1, Chebyshev-like lifting.

Input. Degree n and the coefficients p_{0}, p_{1}, …, p_{n} of the polynomial p(z) of (1.1) having zeros z_{1}, …, z_{n}, where p_{0}p_{n} ≠ 0.

Output. The coefficients p̂_{0}, p̂_{1}, …, p̂_{n} of the polynomial p̂(x) whose zeros are \frac{z_{i}+z_{i}^{−1}}{2} , i = 1, …, n.

Computation.

1. Compute the values α_{i} = p(ζ^{i}), i = 0, …, 2n + 1, where ζ = exp(\frac{2π\sqrt{−1}}{2n+2}) is a primitive (2n + 2)-nd root of 1.

2. Compute the values β_{i} = α_{i}α_{2n+2−i}, i = 0, …, n, which the function p(z)p(z^{−1}) takes on at the (2n+2)-nd roots of 1. These values coincide with the values that the polynomial p̂(x) takes on at the (Chebyshev-like) points x_{i} = (ζ^{i} + ζ^{−i})/2 = cos(πi/(n + 1)), i = 0, …, n.

3. Recover the coefficients of p̂(x) by interpolating to the n-th degree polynomial having the values of β_{i} at the points cos(πi/(n + 1)), i = 0, …, n.

Correctness of this algorithm immediately follows from the observations that p̂(x) is a polynomial of a degree at most n and that the values of p(z)p(z^{−1}) at the (2n + 2)-nd roots of 1 coincide with the values of the polynomial p̂(x) at the points x_{i} = (ζ^{i} + ζ^{−1})/2 = cos(πi/(n + 1)), i = 0, …, n. Algorithm 2.1 belongs to the class of the evaluation-interpolation algorithms (cf. [BP94], ch. 1). It requires O(n log n), n + 1, and O(n log n) ops at its stages 1, 2, and 3, respectively.

Stage 1 is handled by FFT, and the known algorithms of [CHQZ87], [For95], and [P98] for Chebyshev-like interpolation in O(n log n) ops can be used at stage 3; [BP96] relied on the pessimistric cost bound, O(n log^{2} n), and extended it to the Chebyshev-like lifting algorithm. Now, with the refinement of the Chebyshev interpolation cost, the estimated overall cost of Chebyshev-like lifting immediatley goes down to O(n log n), versus the one claimed in [BP96]:

Proposition 2.1 O(n log n) ops suffice to perform Chebyshev-like lifting algorithm 2.1.

Recursive application of algorithm 2.1 enables us to isolate from each other two groups of the zeros of a polynomial. In the next section, we will describe an alternative algorithm supporting proposition 2.1. Namely, such an application moves towards 1 the zeros of p(x) having positive real parts and moves towards –1 the zeros of p(x) having negative real parts. When the isolation is achieved, we may apply the known splitting algorithms of [C96] or [Sc82], in order to split the resulting polynomial into two factors whose zeros are separated by the imaginary coordinate line L.

Our next goal is the descending from such splitting to splitting the original polynomial p(z) over the line L. To simplify the notation, let us assume that already the polynomial p̂(x) has been split into two factors,

p̂(x) = F̂(x)Ĝ(x),  (2.3)

where all the zeros of F̂(x) [respectively, Ĝ(x)] have positive (respectively, negative) real parts. Then we will split p(z) into two factors as follows:

Algorithm 2.2, Chebyshev-like descending.

Input. Polynomials p(z) of (1.1), F̂(x) and Ĝ(x)of (2.3), of degrees n, k, and n − k, respectively.

Output. Two polynomials, F (z) of degree k and G(z) of degree n − k, satisfying p(z) = F (z)G(z) where the real parts of all the zeros of F (z) are positive and of all the zeros of G(z) are negative.

Computation. Compute and output F (z) being the greatest common divisor (gcd) of F̃(z) = z^{k}F̂(\frac{z+z^{−1}}{2} ) and p(z). Then compute and output G(z) = p(z)/F (z), so that G(z) = gcd(p(z), G̃(z)) and G̃(z) = z^{n−k}Ĝ(\frac{z+z^{−1}}{2} ).

To show correctness of the algorithm, recall that F̂(\frac{z+z^{−1}}{2} )Ĝ(\frac{z+z^{−1}}{2} ) = F̂(x)Ĝ(x) = p̂(x) = p(z)p(z^{−1}); therefore, p(z)p_{r}(z) = p(z)p(z^{−1})z^{n} = z^{k}F̂(\frac{z+z^{−1}}{2} )(z^{n−k}Ĝ(\frac{z+z^{−1}}{2} )) = F̃(z)G̃(z). It remains to observe that the real parts of all the zeros of the polynomials F̂(x) and F̃(z) are positive, whereas the real parts of all the zeros of Ĝ(x) and G̃(z) are negative. □

We may obtain the polynomial F (z) by computing the (k, n−k) entry (F (z), G̃(z)/G(z)) of the Padé approximation table for the formal power series defined by the ratio p(z)/G̃(z).

Computation of F (z) both as the gcd and from the entry of the Padé table can be done in O(n log^{2} n) ops, by means of the Euclidean algorithm [BGY80]. The computation of F (z) from the Padé table can be also reduced to solving a Toeplitz linear system of k or n − k equations; the solution of such a system costs O(k log^{2} k) or O((n − k) log^{2}(n − k)) ops, respectively [BGY80]. The techniques of [P96] enable us to express the condition of this system via the isolation ratio of the maximum zero-free annulus on the complex plane \{w\} that separates the zeros of F (\frac{1−w}{1+w}) from ones of G(\frac{1−w}{1+w}). (The isolation ratio is defined as the ratio of the radii of the 2 boundary circles of the annulus (cf. [P95], [P96]).)

3 Correlation between Chebyshev-like and Graeffe’s processes.

Chebyshev-like lifting is quite similar to Graeffe’s process:

p_{0}(x) = p(x)/p_{n}, p_{i+1}(x) = (−1)^{n}p_{i}(\sqrt{x})p_{i}(−\sqrt{x}), i = 0, 1, ….  (3.1)

(Actually, ”Graeffe’s process” was first discovered by Dandelin and soon thereafter independently by Lobachevsky.) Let us formalize this similarity. Recall that the i-th step of (3.1) squares the zeros of p_{i}(x), that is,

p_{i}(x) = p_{n,i} ∏_{j=0}^{n}(x − z_{j}^{2^{i}}), i = 0, 1, …,  (3.2)

whereas (2.2) transforms z_{j} into (z_{j} + z_{j}^{−1})/2, j = 1, …, n. Let us map z into a variable w according to the next map:

z = \frac{1 − w}{1 + w}, w = \frac{1 − z}{1 + z} ,  (3.3)

and let us similarly map x into y:

x = \frac{1 + y}{1 − y} , y = \frac{x − 1}{x + 1}.  (3.4)

Then the map of z_{j} into z_{j}^{2} can be obtained as a composition of the three maps, (3.3), (3.4), and x = \frac{w+w^{−1}}{2} . Namely, we may write

x_{j} = \frac{1 + y_{j}}{1 − y_{j}} ,

w_{j} = \frac{1 − z_{j}}{1 + z_{j}} ,

x_{j} = \frac{w_{j} + w_{j}^{−1}}{2} = (\frac{1 − z_{j}}{1 + z_{j}} + \frac{1 + z_{j}}{1 − z_{j}})/2 = \frac{1 − z_{j}^{2}}{1 + z_{j}^{2}} ,

which implies that y_{j} = z_{j}^{2}.

Likewise, the map of z_{j} into (z_{j} + z_{j}^{−1})/2 can be obtained as a composition of the three maps, (3.3), (3.4), and y = w^{2}. Indeed, we write

x_{j} = \frac{1 + y_{j}}{1 − y_{j}} , y_{j} = w_{j}^{2}, z_{j} = \frac{1 − w_{j}}{1 + w_{j}}

and similarly deduce that x_{j} = \frac{z_{j}+z_{j}^{−1}}{2} .

Let us rewrite the maps (3.3) and (3.4) as follows:

z = \frac{1 − w}{1 + w} = \frac{2}{1 + w} − 1, w = \frac{1 − z}{1 + z} = \frac{2}{1 + z} − 1,

x = \frac{1 + y}{1 − y} = \frac{2}{1 − y} − 1, z = \frac{x − 1}{x + 1} = 1 − \frac{2}{x + 1},

reducing them to scaling, shifts and inversion of the variable. For the respective transformations of a fixed polynomial p(z), (implied by scaling, shift and inversion of the variable), we need O(n), O(n log n) (cf. e. g. [BP94], page 15), and 0 ops, respectively (the inversion of the variable amounts to reversion of the order of the polynomial coefficients). Squaring of the polynomial zeros is supported by (3.1), at the cost of O(n log n) ops. This enables us to perform a single step (2.2) in O(n log n) ops, which is an alternative derivation of proposition 2.1.

4 Generalization of lifting process.

Graeffe’s and Chebyshev-like processes can be generalized to define a polynomial p̃(x) = p̃_{r}(x) of degree at most n with zeros x_{j} = r(z_{j}), j = 1, …, n, for any fixed rational function x = r(z) = ν(z)/δ(z), ν(z) and δ(z) being two given relatively prime polynomials, δ(z) is monic, so that r(z) is a polynomial if and only if δ(z) = 1. We will specify two expressions for p̃(x) = p̃_{r}(x) via p(x) and r(x):

p̃(x) = R(x) = R(p(z), δ(z)x − ν(z)),  (4.1)

R(u(z), v(z)) denoting the resultant of two polynomials u(z) and v(z) (cf. e.g. [BP94], page 149, for the definition of resultants), and

p̃(x) = det(xI − r(F_{p})),  (4.2)

F_{p} denoting the Frobenius (companion) matrix of the monic polynomial p(z)/p_{n},

F_{p} = ( \begin{matrix}0 & 1 & & 0 \\ 0 & & ... & 0 \\ 0 & & & 1 \\ −p′_{0} & & … & −p′_{n−1} \end{matrix} ) , p′_{i} = p_{i}/p_{n}, i = 0, 1, …, n − 1.

It is immediatlely verifed that in both cases (4.1) and (4.2) we have p̃(x) = ∏_{j=1}^{n}(x − r(z_{j})).

We may evaluate the resultant R(x) of (4.1) for x = aω^{i}, i = 0, 1, …, n, ω denoting a primitive 2^{h}-th root of 1, for h = ⌈log_{2}(n + 1)⌉, and then compute the coefficients of p̃(x) by applying the inverse of FFT. The cost of this computation is dominated by the cost O((d + n)n log^{2}(d + n)) of the 2^{h} evaluations of the resultant R(x), where

d = max{deg(ν(z)), deg(δ(z))}.  (4.3)

An alternative algorithm relies on (4.2). We first compute the matrix r(F_{p}) (which belongs to the matrix algebra A_{p} generated by the matrix F_{p}) and then compute its characteristic polynomial, p̃(x) of (4.2). An addition, a subtraction and a multiplication in A_{p} require O(n log n) ops, the inversion of a matrix in A_{p} takes O(n log^{2} n) ops [C96]. Therefore, the computation of the matrix r(F_{p}) involves O((d+μ log n)n log n) ops, for d of (4.3), where μ = 0 if δ(x) = 1, μ(x) = 1 otherwise.

Since every matrix of A_{p} has a displacement rank at most 2 [C96], the evaluation of the characteristic polynomial of such a matrix takes O(n^{2} log n) ops (cf. [P92] or [BP94], pages 189–190). In particular, this cost estimate applies to the matrix r(F_{p}) ∈ A_{p}, so that O((n + d)n log n) ops suffice for computing the coefficients of p̃(x) based on (4.2).

If the next objective is the splitting of the polynomial p̃(x) into the product of two factors by means of algorithm 4.1 of [BP96] or an algorithm of [C96], then it is sufficient to use the implicit representation of p̃(x) by the matrix r(F_{p}), which can be computed at the smaller cost of O((d + μ log n)n log n).

The descending process, from the factors of p̃(x) to the factors of p(z), is more involved in the case of the general rational function r(x). If, however, d is small relative to n and if all the n zeros x_{j} = r(z_{j}) of p̃(x) have been approximated, we may recover the approximation to z_{j} for a fixed j among the d candidates z_{j}^{(i)}, i = 1, …, d.

To achieve this, we may solve in z_{j} the polynomial equation δ(z_{j})x_{j} = ν(z_{j}) for each j = 1, …, n, then compute p(z_{j}^{(i)}) for all i and j, and finally select the n zeros of p(z). The computation of the values p(z_{j}^{(i)}) requires a high precision in the case where some other zeros of p(z) lie near z_{j}^{(i)}.

References

[BGY80] R. P. Brent, F. G. Gustavson, D. Y. Y. Yun, Fast Solution of Toeplitz Systems of Equations and Computation of Padé Approximations, Journal of Algorithms, 1, pp. 259–295, 1980.

[BP94] D. Bini, V. Y. Pan, Polynomial and Matrices Computations, vol. 1: Fundamental Algorithms, Birkhaeuser, Boston, 1994.

[BP96] D. Bini, V. Y. Pan, Graeffe’s, Chebyshev-like, and Cardinal’s Processes for Splitting a Polynomial into Factors, Journal of Complexity, 12, pp. 492–511, 1996.

[C96] J. P. Cardinal, On Two Iterative Methods for Approximating the Roots of a Polynomial, Proc. of Workshop, Math. of Numerical Analysis: Real Number Algorithms, (Park City, Utah, 1995), J. Renegar, M. Shub, and S. Smale editors, vol. 32 of Lecture Notes in Applied Math., pp. 165–188, Amer. Math. Society Press, Providence, R.I. 1996.

[CHQZ87] C. Canuto, M. Y. Hussaini, A. Quateroni, T. A. Zang, Spectral Methods for Fluid Dynamics, Springer, 1987.

[For95] B. Fornberg. A Practical Guide to Pseudospectral Methods. Cambridge University Press, 1995.

[P89] V. Y. Pan, Fast Evaluation and Interpolation at the Chebyshev Set of Points, Applied Math. Letters, 2, 3, pp. 255–258, 1989.

[P92] V. Y. Pan, Parametrization of Newton’s Iteration for Computations with Structured Matrices and Applications, Computers & Mathematics (with Applications), 24, 3, pp. 61–75, 1992.

[P95] V. Y. Pan, Optimal (up to Polylog Factor) Sequential and Parallel Algorithms for Approximating Polynomial Zeros, in Proc. 27th Annual ACM Symposium on the Theory of Computing, pp. 741–750, ACM Press, New York, 1995.

[P96] V. Y. Pan, Optimal and Nearly Optimal Algorithms for Approximating Polynomial Zeros, Computers and Math. (with Applications), 31, 12, pp. 97–138, 1996.

[P98] V. Y. Pan, New Fast Algorithms for Polynomial Interpolation and Evaluation on the Chebyshev Node Set, Computers & Math. (with Applications), 35, 3, 125–129, 1998.

[Sc82] A. Schönhage, The Fundamental Theorem of Algebra in Terms of Computational Complexity, Manuscript, Dept. of Math., Univ. of Tübingen, Tübingen, Germany, 1982.