AN EFFICIENT RELIABLE ALGORITHM FOR THE

AN EFFICIENT RELIABLE ALGORITHM FOR THE APPROXIMATION OF ALL POLYNOMIAL ROOTS BASED ON THE METHOD OF D. BERNOULLI ¨ HERBERT MOLLER Abstract. The present versions of Bernoulli’s method use quotients of successive sums of powers of the roots to extract dominant roots. These algorithms can be ameliorated considerably through two new basic results, namely, a recursion formula for the quotients just mentioned and a representation of the minimal absolute value of the roots only using the quotients calculated earlier. In the case of (nearly) equimodular roots the second result is combined with the search for local minima of the absolute value of the polynomial over the circle which has the approximated minimal modulus as radius. Shifting to these minimum points turns the handicap of Bernoulli’s method to an optimal fitting. Error bounds and stopping criteria are obtained with the help of “Laguerre disks”. The reliability is ´ n’s method. A supplementary stage guaranteed by an adaptation of Tura is developed to separate all roots in disjoint Laguerre disks and to ensure at least quadratic convergence of Newton-Raphson iteration. Further details, programs and visualizations can be found on the webpage [7] .

1. Introduction The design of the complete algorithm is based on twelve theorems. Since most of them share the same conditions and abbreviations, we first state these specifications. Then we summarize the contents of this paper and put them in order with the efficient reliable polynomial rootfinding methods already known. Specifications. Let Jm : = {1, . . . , m} for m ∈ N and let z1 , . . . , zk be complex numbers with zg 6= zh for g, h ∈ Jk , g 6= h, and 0 < |z1 | ≤ . . . ≤ |zk | . A P i polynomial p(z) is given by p(z) = m i=0 ai z with ai ∈ C (i ∈ Jm ) , am 6= 0 and a0 = −1 . It denotes either the original polynomial or a polynomial obtained Q by deflation. In the product representation p(z) = am kj=1 (z − zj )mj with mj ∈ N the roots z1 , . . . , zk are complex numbers with properties as stated above. The index j always ranges in Jk , the index n is in N0 . Throughout this paper the following abbreviations are used: sn : =

k X

mj zj−n ,

j=1

Date: October 22, 2000. 1991 Mathematics Subject Classification. Primary: 65; Secondary: 68. Key words and phrases. Polynomial rootfinder, polynomial zerofinder. 1

2

¨ HERBERT MOLLER

k(n) : = min {k ∈ N | k > n and sk 6= 0} , qn : =

sn sk(n)

,

Z : = {z1 , . . . , zk } , Z1 : = {z ∈ C | p 0 (z) = 0 } , µj : = min {d ∈ R | There is a g ∈ Jk \ {j} with d = |zg − zj | } , µ : = min {µ1 , . . . , µk } and nmo p (u) ht (u) : = u − t 0 for t ∈ Jm ∪ and u ∈ C \ Z1 . p (u) 2 In Section 2 we start with Newton’s formulas connecting the sums of powers sn with the coefficients of p(z). Since the sums of powers are fundamental for the mathematical derivation given here, the method is called “Sums of Powers Algorithm” (SPA). The essential theoretical results were published in [6]. But many of the facts which let the corresponding rootfinding program fulfill the ¨ vel who planned criteria of [4] have been developed and tested by D. Gunstho to write a dissertation in this field. For example, in Theorem 1 he has found the explicit a priori lower bound for the number of steps needed to get the desired z1 accuracy depending on a fixed upper bound of . z2 The basic recursion formula for qn is stated in Theorem 2. This is one of three results which are also contained in [6] and for which the proofs are not given here. The english translation of these proofs can be found on the webpage [7]. In Section 3 the error bounds and two stopping criteria are developed. They are connected with Laguerre disks which, by definition, are the closed disks having any u ∈ C \ (Z ∪ Z1 ) and hm (u) as endpoints of a diameter. It is known from function theory that every Laguerre disk contains at least one root of p(z) . Moreover, the radius is bounded from above and below by constant multiples of the distance between {u} and Z if the distance is sufficiently small. In Section 4 the “second phase” of the SPA is established which solves the z1 convergence problem of the Bernoulli method if is close or equal to 1 . z2 At first, the minimal absolute value of the roots is approximated by a number r only using the quotients qn calculated in the “first phase”. Afterwards, small local minima of t → | p(r cos t + ir sin t)|2 , t ∈ [0, 2π] are determined with the help of the Fourier expansion of this function. Then, after shifting to the corresponding points in C the SPA gives quotient sequences which converge sufficiently fast to roots with absolute values close to r. Section z1 5 contains an amelioration of this method if is quite close to 1 . z2 ´ n are used to handle the In Section 6 inequalities of Buchholtz and Tura “ultimate case”. In this way, the SPA becomes absolutely reliable with good

APPROXIMATION OF POLYNOMIAL ROOTS

3

worst case complexity. Moreover, the interplay of quotient sequences and Laguerre disks enables us in Section 7 to separate all roots in disjoint disks and to overcome the deflation error. The last theoretical part is devoted to Newton-Raphson iteration. The same condition which has been used for the first stopping criterion in Section 3 now ensures at least quadratic convergence of the sequence (un )n generated by un+1 : = h1 (un ) , n ∈ N 0 , if the single root contained in the Laguerre disk of u0 has multiplicity 1. Section 8 contains further hints for the implementation of the SPA both in a compiler version and in a CAS. Several survey articles and bibliographies on polynomial rootfinding have been published during the last fifteen years (e.g. [1]). Therefore, only some supplementary remarks shall help to compare the SPA with present versions of efficient reliable methods already known. The properties of the SPA are similar to those of the Jenkins-Traub algorithm [3] which indeed starts with a reformulation of Bernoulli iteration. There are two major differences, namely, the SPA doesn’t directly lead to quadratic convergence, but it makes up for this by not depending on any unsure decisions. In fact it can be used to determine the critical number s in Stage Two of the Jenkins-Traub algorithm. As well, it is able to produce satisfactory starting values for the Dochev-Durand-KernerWeierstrass algorithm [2] which simultaneously determines approximations for all roots. Since the SPA immediately works with numbers which have the same order as the roots, its efficiency is quite good compared with other methods which use disks or circles, namely, the Cohn-Lehmer-Schur algorithm [5] which encloses ¨ nhage-Weyl the roots by two-dimensional bisection and the Henrici-Scho algorithm [8] which factors polynomials with the help of contour integrals. 2. Properties of sums of powers If p(z) is a polynomial as specified above, then the sums of powers sn for n ∈ N can be determined recursively with Newton’s formulas:  n−1 P   aj sn−j + nan when n ≤ m ,  j=1 (1) sn = m P   aj sn−j when n > m .  j=1

Theorem 1. Let z1 , . . . , zk be complex numbers as specified above with |z1 | < |z2 | and |z1 | ≤ b , where b may be chosen suitably. If ε and κ are real numbers with 0 < ε ≤ 2b and zz12 ≤ κ < 1, then |qn − z1 | < ε holds for all n ∈ N0 with " # 1 1 −1 (2) n> log log 4mb . κ ε Proof. If k = 1, then we have sn = mz1−1 and therefore qn = z1 for every n ∈ N0 . In the case k > 1 we write sn = z1−n (m1 + m2 ζ2n + . . . + mk ζkn ) with

¨ HERBERT MOLLER

4

ζj : = zz1j for j = 2, . . . , k . Since | m2 ζ2n + · · · + mk ζkn | ≤ mκn , it follows that log(2m) | m1 + m2 ζ2n + · · · + mk ζkn | ≥ 12 for all n ≥ = : n0 . In particular, we − log κ get sn 6= 0 and therefore k(n) = n + 1 for n ≥ n0 . For these n, we have ! !−1 k k X X qn − z 1 = z 1 mj ζjn (1 − ζj ) m1 + mj ζjn+1 j=2

j=2

from which |qn − z1 | ≤ |z1 | 4mκn follows. Hence for every n ∈ N0 with n > log 4mbε−1 ≥ n0 , we get |qn − z1 | < ε . − log κ 1 Based on extensive tests, the value log 0.71...

−1

= 3.0 has been chosen for

the implementation of the SPA. Since (sn )n or s−1 may be unbounded, even if (qn )n is bounded, a new k(n) n

recursion formula for qn is decisive for the efficiency of the SPA. For instance, in the second case of (1) with sk 6= 0 for all k ∈ N it follows that −1 −1 qn−1 = sn+1 s−1 n = a1 + a2 sn−1 sn + · · · + am sn+1−m sn

= a1 + a2 qn−1 + · · · + am qn−1 qn−2 · · · qn+1−m = a1 + qn−1 a2 + qn−2 (· · · ) . Theorem 2. Let p(z) be a polynomial as specified above. Then k(0) = −1 min{k ∈ Jm | ak 6= 0} and q0 = m k(0) ak(0) . Assume k(s) and qs to be known for s = max {0 , n − m} , . . . , n − 1 . Then k(n) and qn can be calculated in the following way: For s = max {0 , n − m} , . . . , n − 1 and for  when  qs (s+t) 1 when qs := −1  t k(0) ak(0) when at when qs 6= 0 , (s+t) at := 0 when qs = 0 .

t ∈ Jm , introduce s > 0 and s > 0 and s = 0,

qs = 6 0, qs = 0 , and

In the case k(n − 1) = n define min{k,m}

(3)

pn,k : = ak−n +

X t=k−n+1

(k) at

n−1 Y

! qs(k)

s=k−t

for k = n + 1, . . . , n + m . Then k(n) = n + min {j ∈ Jm | pn,n+j 6= 0} and qn = p−1 n,k(n) hold. For k(n − 1) > n, it follows that k(n) = k(n − 1) and qn = 0 .


5

3. Error bounds and stopping criteria The error bounds and stopping criteria are based on the following result of E. Laguerre (Œuvres, Vol. 1, Paris 1898): Theorem 3. If u ∈ C \ (Z ∪ Z1 ) , then every circle passing through u and hm (u) has the property that either there are elements of Z inside and outside the circle or that all numbers of Z lie on the circle. Proof. We start from the well-known formula k

(4)

p 0 (u) X mj = p(u) u − zj

for u ∈ C \ Z .

j=1

P a a With v : = hm (u), we get 0 = kj=1 mj u−z − u−v for every a ∈ C . If C is j a circle through u and v, then there exists a unique number a ∈ C with |a| = 1, a a − u−v , w ∈ C \ {u} maps such that the injective function f : = w → u−w v to 0, C \ {u} to the real axis and the interior of C to the upper half-plane P {z ∈ C | Im z > 0} . Moreover we have 0 = kj=1 mj f (zj ) . Therefore either there is at least one element of f (Z) in the upper half-plane and at least one in the lower half-plane or f (Z) contains only real numbers. Then the mapping properties of f yield the corresponding result with respect to the circle C . Especially, the closed disk n o Lu = Lu (p) : = w ∈ C w − h m2 (u) ≤ u − h m2 (u) , which we call Laguerre disk of u, contains at least one root of p . The following theorem proves that the radius of Lu is less than m |u − zi | , if u is sufficiently close to zi ∈ Z. From this we obtain the stopping criterion in the first two stages. Theorem 4. Let p(z) be a polynomial as specified above with k ≥ 2 . Then p(u) ≤ 2m − 1 |u − zj | (5) |u − zj | < m 0 p (u) 2mj − 1 and Lu ∩ Z = {zj } hold for every j ∈ Jk and for all u ∈ C with 0 < |u − zj | ≤ 1 2m µj . Proof. From (4) for u ∈ C \ Z1 , we get p(u) m = m 0 |u − zj | with p (u) |mj + Sj |

Sj : =

X g∈Jk \{j}

mg

u − zj . u − zg

The assumption on u leads to 2m |u − zj | ≤ |zg − zj | ≤ |u − zg | + |u − zj | . u−zj 1 This gives u−zg ≤ 2m−1 for every g ∈ Jk \ {j} . Now, with (6)

|Sj | ≤

m − mj , 2m − 1

¨ HERBERT MOLLER

6

we have

p(u) m 2m − 1 ≤ m 0 |u − zj | ≤ |u − zj | . p (u) mj − |Sj | 2mj − 1

The proofs of the first inequality and of the intersection property may be found in [6] . They are carried out by contradiction assuming that the interior of Lu contains more than one root or none. Since the a priori number of steps determined by (2) covers a wide range of convergence quotients zz21 , it is desirable to have a simple pre-check instead of p(qn ) ≤ ε from (5). This can be directly testing the stopping condition m 0 p (qn ) arranged with the help of the following theorem. Theorem 5. Let p(z) be a polynomial as specified above with k ≥ 3 . If |z1 | < |z2 | < |z3 | , then the sequences p 0 (qn ) qn+1 − qn+2 and (qn+1 − qn ) qn − qn+1 n p (qn ) n are convergent with limits zz12 and m1 zz12 − m1 respectively. From the second result it follows that p (qn ) (7) |qn+1 − qn | ≤ 2 m 0 p (qn ) holds for sufficiently large n . g P for g ∈ N , we have Proof. Using the abbreviation τg : = m2 + kj=3 mj zz2j τ − ζ2 τn+1 from the proof of Theorem 1 and we get qn − z1 = z1 ζ2n n n+1 m1 + ζ2 τn+1 n+1 2 2 τn τn+2 − τn+1 n m1 τn − 2ζ2 τn+1 + ζ2 τn+2 + ζ2 qn − qn+1 = z1 ζ2 . m1 + ζ2n+1 τn+1 m1 + ζ2n+2 τn+2 q − qn+2 Since lim τn = m2 and lim ζ2n = 0, it follows that n+1 qn − qn+1 n is convern→∞ n→∞ z gent with limit ζ2 = z12 . From (4), we obtain k X mj p 0 (qn ) qn+1 − qn = m1 + (qn+1 − qn ) . p (qn ) qn − z 1 qn − z j j=2 qn+1 − qn As above, we find that we have qn − z1 n converges to ζ2 − 1 . Therefore 0 limn→∞ (qn+1 − qn ) pp(q(qnn)) = m1 zz21 − m1 . With m1 ≤ m − 2 and zz12 < 1, it p (qn ) for follows that there exists a number n1 such that |qn+1 − qn | ≤ 2 m 0 p (qn ) all n ∈ N with n ≥ n1 .

(qn+1 − qn )


7

4. The minimum method Let " (8)

nε : =

1 log κ

−1

1 log 4mb ε

#

be the bound of (2) in Theorem 1. Even if z1 is not approximated sufficiently by the first nε terms of the sequence (qn )n , they can be used to approximate |z1 | in the following way. Theorem 6. Let z1 , . . . , zk be complex numbers as specified above. If h(n) is defined recursively by h(0) : = k(0) and h(j + 1) : = k (h(j)) for every j ∈ N 0 1 h(n) and if ms−1 is abbreviated by cn for n ∈ N , then |z1 | is the limit of the h(n) sequence (σn )n with n o σn : = min a ∈ R There exists a j ∈ N with h(j) ≤ n and a = cj for n ∈ N which consists of the “successive minima” of (cn )n . Using the abbreviation λn : = h(n) log cn , the terms log cn can be calculated recursively by λ0 = log |q0 | and λn+1 = λn + log qh(n) for every n ∈ N0 . An elementary proof is given in [6] . Since the log-function is monotone λn increasing, the successive minima of (cn )n correspond with those of h(n) . n

Every root of p(z) generates a spherical point of the “analytical landscape” (x, y, w) ∈ R3 | w = | p(x + iy)|2 which has no other hollows. Therefore we use the following method to obtain origin shifts of the current polynomial to good approximations of roots close to the “minimal distance circle” with radius r : = σnε given by Theorem 6. If we look at the intersection of the analytical landscape and the cylinder over the circle with radius r and with the origin as midpoint, we see that every sufficiently small local minimum of the function πr : = t → | p(r cos t + ir sin t)|2 , t ∈ [0, 2π] belongs to the root in the hollow which produces the minimum. Now it suggests itself to proceed by approximating these minimum points. First we determine the coefficients fk : = rk

m−k X j=0

Re (aj+k aj ) r2j and gk : = rk

m−k X j=0

Im (aj+k aj ) r2j ,

¨ HERBERT MOLLER

8

k = 1 , . . . , m, of the Fourier expansion πr (t) − f0 =

m X 3π fk 2 cos(kt) + gk 2 cos + kt . 2 k=1

Then we calculate the values of πr (t) − f0 at equidistant arguments t using a π with n = 0 , . . . , 4N −1, which can be generated recursively table of 2 cos n 2N by 2 cos ((n + 2) ∆) = (2 cos ∆) (2 cos((n + 1) ∆)) − 2 cos (n∆) π . To enable easy doubling of this table using for n = 0 , . . . , N − 3 with ∆ : = 2N the same recursion, we take powers of 2 for N weakly depending on the degree of p . Comparing every three successive values of πr (n∆) − f0 for n = 0 , . . . , 4N − 1, we find approximations of the minima. Let x1 , . . . , xe be these numbers arranged so that πr (x1 ) ≤ . . . ≤ πr (xe ) holds. Then we choose r cos xj + 3π ir cos 2 + xj for j = 1 , . . . , e as origin shifts for the current polynomial as long as the first nε terms of the corresponding sequence (qn )n lead to a sufficient approximation of a root. Since at least those roots which have caused the bad behavior of the quotient sequence in the first phase are quite close to −1 the minimal distance circle, we diminish the “stepfactor” log κ−1 of nε in (2) to 2.0. For random polynomials with normal distribution of coefficients, it is known that most of the roots are uniformly distributed around the unit circle (cf. [1], p.173). Furthermore, deflation with a large zero may cause instability. Therefore we introduce a method which we call “dual deflation”. At first, we deflate p(z) with the approximations lying in the central disk with radius 1 − √0.1 . m This gives a polynomial p1 (z) and a set A1 of approximations with associated multiplicities. Afterwards, we apply the same steps with the polynomial p2 (z) : = z m p z −1 which has the roots z1−1 , . . . , zk−1 . Let A2 be the set of reciprocals of the approximations found with p2 (z) and p3 (z) be the polynomial obtained by deflating p1 (z) with the elements of A2 . For the remaining roots, we use the minimum method to determine the local minima of |p3 (z)|2 on the unit circle. Then we proceed by shifting the current polynomial p1 (z) to these minimum points. Continuing as above, it may happen that a root belonging to an element of A2 is approximated once again. In this rare case, the correct multiplicity is determined by local deflation of the shifted original polynomial. If the sum of the multiplicities of the approximations in A1 ∪ A2 is less than m − 2 , then p3 (z) is deflated with the new approximations and the minimum method is repeated until all roots are approximated with the exception of at most two.


9

5. The meteorite method If |z2 | is equal or very close to |z1 | , then the minimum method described above represents an optimal balancing complement of the pure quotient method. In this case the numbers qn visualized as points of the complex plane don’t show any regularity. If p is a polynomial of small degree (say up to 15), we often observe that the points qn lie on a spiral which closes very slowly. Trying to use this extra information, we find that the centers of the circles through every three successive points qn lie in a small parallel strip containing z1 and z2 . To prove this observation we suppose that z1 , . . . , zk are complex numbers as specified above with k ≤ 3 and |z1 | < |z2 | < |z3 | . Then qn can be approximated by (9)

z1

m1 + m2 ζ2n = z1 + (z2 − z1 ) m1 + m2 ζ2n+1 1+

1 ζ2−n−1 ,

m1 m2

for sufficiently large n ∈ N . Therefore we only have to consider sequences −1 with c ∈ R\{0} . Writing ζ2−1 = : ρ exp (2πiα) with α, ρ ∈ R , 1 + c ζ2−n ρ > 1, 0 < α
1 , and let xj,k + i yj,k with xj,k , yj,k ∈ R be the center of the circle through u(j +kq), u(j +kq +1) and u(j +kq +2) for j ∈ {0 , . . . , q −1} and k ∈ N 0 .

p p sin 2π(j + 1) cos 2π(j + 1) q q If Sj : = ρ2 − 1 , Cj : = ρ2 − 1 and p p 2 sin 2π 2 sin 2π q q −1 Tj : = Sj2 + ρ2 , then (2 xj,k − 1)2 (2 yj,k + Cj Sj Tj )2 − =1 ρ2 Tj ρ2 Cj2 Tj2

(10)

holds for every k ∈ N 0 . The proof is long but straightforward. We use the representation of the center coordinates as quotients of three-rowed determinants. With the abbreviation Rj,k : = c ρ qk+j , we get xj,k =

1 − Sj Rj,k Cj Rj,k and yj,k = 2 2 2 1 − 2 Sj Rj,k − ρ Rj,k 1 − 2 Sj Rj,k − ρ2 Rj,k

¨ HERBERT MOLLER

10

which may be considered as a parametric representation of points on a curve, the parameter being Rj,k . Therefore (10) is obtained by simply eliminating Rj,k . Equation (10) means that all centers lie on hyperbolas. The absolute values of the slopes of the asymptotes and of the ordinates of the central points 2 ρ2 − 1 ρ2 − 1 and respectively. Since the are bounded by 2 ρ sin 2π pq 8 ρ2 sin2 2π pq sine function is continuous, we can replace

p q

by α . Therefore the centers of

the circles through three successive terms of (qn )n are quite near to the line {z1 + (z2 − z1 ) t | t ∈ R} which passes through z1 and z2 . We use this result with polynomials of degree from 3 to 14 to determine approximations of z1 and z2 by calculating the intersection points of the minimal distance circle with a simplified adjustment line through suitable subsets of the centers. We call this procedure meteorite method , because the visualization of the centers and of the minimal distance circle looks like a meteorite approaching a planet. ´n 6. Inequalities of Buchholtz and Tura In Theorem 6, we proved that |z1 | is the limit of a sequence which can be calculated with the help of the quotient sequence (qn )n . Since we stop at the relatively small number nε of terms approximating the radius of the minimal distance circle, we can’t be sure to find a root in this way, despite the good numerical experience. The following theorem contains the adaptation of ´ n to our terminology. They fill the gap inequalities of Buchholtz and Tura in the first stage of the SPA, namely, we use the inequalities of the theorem to implement an “ultimate case”. Theorem 8. Let z1 , . . . , zk be complex numbers as specified above and let σn for n ∈ N be defined as in Theorem 6 . Then (11)

1

5− N σmN < |z1 | ≤ σmN

holds for every N ∈ N . The case N = 1 is due to Buchholtz. The general case was derived by ´ n from Buchholtz’s inequality. The proofs can be found in [9]. Tura With respect to (11), we call every set o n 1 TN (ξ) : = z ∈ C 5− N σmN < |z − ξ| ≤ σmN with N ∈ N ´ n annulus , if σmN is calculated with the quotient sequence belonging to Tura the current polynomial shifted with ξ . Then the ultimate case is based on the ´ n who developed a whole method of rootfinding with it. following idea of Tura


11

´ n annulus T16 (0) and let ξi with i ∈ J12 Let r be the outer radius of a Tura be equidistant points on the outer boundary of T16 (0) . Then at least one of the ´ n annuli T16 (ξi ) has an outer radius less than 53 r . Dual deflation twelve Tura √ 0.1 √ < 1.11 1 + . enables us to start with r ≤ 16 5 1 + √0.1 m m ´ n annulus contains at least one root and since the minimal Since every Tura ´n distance of the roots is µ > 0 , this procedure, which selects Tura a suitable 1 annulus at every turn, approximates a root after at most O log µ repetitions. In every round we need O m2 operations for the quotients and shifts. The calculation of the logarithms with accuracy O m−2 using the mantissa and the exponent of |qn | can be done with O(m log m) operations. Therefore we 1 3 get the worst case complexity O m log µ for the approximation of all roots in Stage One where µ may be replaced by the desired accuracy ε . 7. Separation of the roots At the end of Stage One we obtain two lists containing the approximations of all roots and the distances of the corresponding Laguerre disks calculated with the original polynomial. If all Laguerre disks are disjoint and if the diameters are less then the desired accuracy, we are ready. Sometimes, this is not true, since by default we start with medium accuracy to get an overview. Due to the deflation error it may happen that we miss roots located outside the Laguerre disks while others are associated with more than one disk. Therefore, the main goal of Stage Two is the amelioration of the approximations so that all roots will be contained in Laguerre disks. Principally, these results are obtained by stepwise increasing the precision of the basic operations and the accuracy of the approximations which can be achieved dynamically with most CAS. Since every Laguerre disk contains at least one root, the ameliorations are always carried out with the SPA after an origin shift of the given polynomial to the center of the Laguerre disk. Moreover, to avoid repeated calculation of the same root we look out for all other Laguerre disks which contain this root. Such a root is stored with the disk data and will be deflated before the SPA is used to search for further roots in that disk. This process is repeated until all zeros in each Laguerre disk are determined. Since the number of roots contained in a disk can be calculated with the help of the Cohn-Schur transformation (cf. [5]), we have a precise stopping condition. If having passed all Laguerre disks we miss some roots due to the deflation error in Stage One, we go on with deflating the given polynomial with the known roots one after the other. Then the resulting polynomial is treated as in Stage One. The following theorem shows that by increasing the accuracy, the deflation error can be made as small as it is necessary to locate all roots of the given polynomial in disjoint Laguerre disks.

¨ HERBERT MOLLER

12

Theorem 9. Let f (z) be a polynomial, w be a root of f (z) with multiplicity g , η be a nonzero complex number and δ = δ(w, η) : = min {t ∈ R | There is a v ∈ C with f (v) = η and |v − w| = t }. If |η| is sufficiently small, then δ≤

(12)

2g! |η| |f (g) (w)|

1 g

.

Proof. Let u ∈ C be a solution of f (u) = η with |u − w| = δ. By elementary methods it can be shown that 1 1 f (u) − f (w) − f 0 (w)(u − w) − · · · − f (g) (w)(u − w)g = q(u, w)(u − w)g+1 1! g! with an explicitly known polynomial q(u, w). Here we have f (u) = η, f (w) = f 0 (w) = · · · = f (g−1) (w) = 0 and f (g) (w) 6= 0. If m is the degree of f (z), then Vieta’s root theorem leads to (13)

1

|u − w| = δ ≤ |η| m .

There are well known upper bounds B of |w| depending on the coefficients of 1 f (z). Restricting |η| by an upper bound M > 0 we can take B + M m as an upper bound for |u| by (13). Therefore, for every M > 0, we can find a uniform upper bound K > 0 of |q(u, w)| such that |η − 1 f (g) (w)(u − w)g | ≤ Kδ g+1 for g! all η with |η| ≤ M. It follows that 1 (g) g (14) δ f (w) − Kδ ≤ |η|. g! (g) |f (w)|m |f (g) (w)| If |η| ≤ min . Now, from (14) m , M , then (13) gives δ ≤ 2g!K (2g!K) we obtain δ g

|f (g) (w)| ≤ |η| which is equivalent with (12). 2g!

The SPA uses η in the following way. If w0 is another zero of f (z) and if w0∗ is an approximation of w0 determined with the stopping condition f (w0∗ ) (15) m 0 ∗ ≤ ε , f (w0 ) where ε is the desired upper bound of |w0 − w0∗ | , then deflation with w0∗ leads to a polynomial f1∗ (z) satisfying f (z) = (z − w0∗ ) f1∗ (z) + f (w0∗ ) . With η : = f (w0∗ ) Theorem 9 gives an upper bound for the deflation error |u − w| = δ , where u is the zero of f (z) − η = (z − w0∗ ) f1∗ (z) corresponding to w.


13

ε . Therefore, by repeated application of By (15) we have |η| ≤ |f 0 (w0∗ )| m Theorem 9, the total deflation error δd (ε) after d deflations fulfills |δd (ε)| ≤

(16)

d X

1

ci ε gi

i=1

with nonnegative constants ci only depending on f and with multiplicities gi which are defined successively by Theorem 9 for the sequence of roots in the order of the corresponding deflations. Now, let w be any zero of p , w∗ be the corresponding zero of a deflation polynomial after d deflations and u be an approximation of w∗ . By stepwise decreasing ε, the inequality (17)

max {ε , |δd (ε)|} ≤

can be satisfied which leads to |u − w∗ | ≤ 1 2m

1 4m

1 µ 4m

µ by (15) and to |w − w∗ | ≤

1 4m

µ

by (16) . Therefore we get |u − w| ≤ µ and Theorem 4 gives Lu ∩ Z = {w}. As mentioned above this enables us to determine disjoint Laguerre disks which contain all roots. If w has multiplicity g > 1, then there may be g different approximations of w which generate up to g Laguerre disks. All these disks contain w but no other roots. Moreover, w is an interior point of the intersection. Therefore, the approximations of w can be ameliorated locally so that they finally lie in only one of the disks which are not adapted to the new approximations. The practical use of condition (17) raises two problems. The deflation error δd (ε) is not known explicitly and the minimal distance µ of the roots can only be approximated in a complicated manner. In [6] a simple lower bound µ0 of µ is given. But, up to now, this result can’t be used efficiently because the method is either very time-consuming or extremely sensitive to roundoff errors. Moreover, the approximation of all roots in Stages One and Two often runs more than twice as fast as the calculation of the rational number (µ0 )2 even using integer arithmetic. Therefore, we never estimate the minimal distance µ of the roots in advance. Instead of that, the precision of the operations and the accuracy of the results are ameliorated stepwise with ε = 10−[ 2+log m] j , j = 1, . . . , J , where the stopping number J is defined through successful separation of the roots which is guaranteed by the above considerations. From (16) and γ : = max {g ∈ N | There is an i ∈ Jd with g = gi } (17) with 1 we get J = O γ log µ , where γ does not really depend on m because the multiplicities gi are determined by successive deflation polynomials. Since we use quotients of power sums, shifts as in Stage One, the worst case and deflations complexity in Stage Two is O m3 log µ1 . Moreover, it should be mentioned that two shifts are as costly as m deflations with decreasing degree.

14

¨ HERBERT MOLLER

8. Newton-Raphson iteration This optional Stage Three has two aims. First, it provides at least quadratic convergence even in the case of multiple roots thus generalizing the method of [6]. Second, it helps to find the correct multiplicities of roots lying in the same Laguerre disk. If a Laguerre disk L contains only one root zj with any multiplicity mj and if L has no common elements with the other Laguerre disks determined in Stage Two, then the minimal distance between L and the other Laguerre disks is a lower bound for µj in the following theorem which states the “normal case”. Theorem 10. Let p(z) be a polynomial as specified above with k ≥ 2. Then hm (u) − zj ≤ m − mj |u − zj | (18) j 2mmj − m 1 µj . holds for every j ∈ Jk and for all u ∈ C with 0 < |u − zj | ≤ 2m The sequence (un )n with u0 : = u and un+1 : = hmj (un ) for n ∈ N0 converges to zj and the convergence is at least quadratic which means that there exists a constant K > 0 such that

(19)

|un+1 − zj | ≤ K |un − zj |2 for all n ∈ N0 .

The error can be estimated a posteriori by m − mj (20) |un+1 − zj | ≤ |un − un+1 | for all n ∈ N0 . 2m(mj − 1) + mj Proof. Similar to the proof of Theorem 4 we get m − mj |Sj | hm (u) − zj = 1 − mj |u−zj | ≤ |u−zj | ≤ |u−zj | j mj + Sj mj − |Sj | 2mmj − m m = : q |u − zj | . With q = 2m1j −1 1 − mj < 1 complete induction leads to 1 µj q n for all n ∈ N0 2m which means that zj is the limit of (un )n . 1 µj , the function hmj (u) Since p 0 (u) 6= 0 for all u ∈ C with |u − zj | ≤ 2m is infinitely differentiable at zj . Therefore, we can find a constant K > 0 such that hmj (u) − hmj (zj ) − h0mj (zj )(u − zj ) ≤ K(u − zj )2 . With hmj (zj ) = zj and h0mj (zj ) = 0 it follows that hmj (u) − zj ≤ K|u − zj |2 . From (18), we have |un+1 − zj | ≤ q |un − zj | ≤ q |un − un+1 | + q |un+1 − zj | and therefore q |un+1 − zj | ≤ 1−q |un − un+1 | which is equivalent with (20). |un − zj | ≤ |u0 − zj | q n ≤

At the end of Stage Two all roots are enclosed in small disjoint Laguerre disks. But sometimes, we do not know the number of different roots in such a disk. This will be achieved by using the revealing behavior of the sequences 1 (un )n defined by 0 < |u0 − zj | ≤ 2m µj and un+1 : = hf (un ), n ∈ N0 , with f ∈ Jmj .


15

Theorem 11. Let p(z) be a polynomial as specified above with k ≥ 2, F > 1 be the sum of the multiplicities of all roots of p(z) contained in a Laguerre disk 1 L, zj be a root in L and u0 be a complex number with 0 < |u0 − zj | ≤ 2m µj . Then the terms of the sequence (un )n defined by un+1 : = hf (un ) for n ∈ N0 fulfill un+2 − un+1 (21) un+1 − un = O (|un − zj |) if f = mj and un+2 − un+1 1 (22) un+1 − un ≥ 2F (2F − 1)(2m − 1) if f ∈ JF ∩ J2mj −1 \ {mj } for all n ∈ N0 with un 6= un+1 . Proof. i) In the case f = mj , we start with |un+2 − un+1 | ≤ |un+2 − zj |+ |un+1 − zj |. From (19), (20) and (18) we have |un+2 − zj | ≤ K |un+1 − zj |2 2 m − mj ≤ K |un − un+1 |2 2m(mj − 1) + mj ≤

K(m − mj )2 (2m − 1)mj (2m(mj − 1) + mj )2 (2mj − 1)m

|un − zj ||un − un+1 |

≤ 2K(m − 1)2 |un − zj ||un − un+1 |. Splitting the right-hand side of (19) we get |un+1 − zj | ≤ K|un − zj |2 ≤ K|un − zj | (|un − un+1 | + |un+1 − zj |) and therefore |un+1 − zj | ≤

K|un − zj | |un+1 − un | ≤ 2K|un − zj ||un+1 − un | 1 − K|un − zj |

for all n ∈ N0 with K|un − zj | ≤ 12 . Now, summation and division gives un+2 − un+1 2 un+1 − un ≤ 2Km |un − zj | for all n ∈ N0 with un 6= un+1 and K|un − zj | ≤ 12 . ii) If mj > 1, we look at the remaining cases f ∈ J F ∩ J2mj −1 \ {mj }. Here, f we write |un+1 − zj | = |hf (un ) − zj | = 1 − m + S |un − zj | with j n X |un − zj | Sn : = mg . |un − zg | g∈Jk \{j} |mj − f | + |Sn | f Using the abbreviation ϕn : = 1 − m + S we get ϕn ≤ ≤ mj − |Sn | j n m−m mj − 1 + |Sn | . Complete induction starting with |S0 | ≤ 2m − 1j from (6) leads mj − |Sn | to m − mj m − mj −1 1 ϕn ≤ mj − 1 + mj − =1− 2m − 1 2m − 1 m

16

¨ HERBERT MOLLER

m−m and |Sn | ≤ 2m − 1j for all n ∈ N0 . Now, from |un+2 − un+1 | ≥ ϕ−1 n − 1 |un+2 − zj | and |un+1 − un | ≤ |un+1 − zj | + |un − zj | it follows that !−1 un+2 − un+1 un+1 − zj −1 un+2 − zj ≥ ϕ−1 1+ n −1 un+1 − un un+1 − zj un − zj −1 −1 −1 = ϕ−1 = (1 − ϕn ) 1 + ϕ−1 n − 1 ϕn 1 + ϕn n −1 |mj − f | + |Sn | mj + |Sn | ≥ 1− 1+ mj − |Sn | |mj − f | − |Sn | (mj − |mj − f | − 2|Sn |) (|mj − f | − |Sn |) = (|mj − f | + mj ) (mj − |Sn |) m−1 m−1 1 ≥ 1−2 1− 2m − 1 2m − 1 (2F − 1)F 1 ≥ . 2F (2F − 1)(2m − 1) Since side of (21) tends to zero, there exists an n0 ∈ N0 such the right-hand un+2 −un+1 1 that un+1 −un < 2F (2F −1)(2m−1) for all n ∈ N0 with n ≥ n0 if f = mj . 9. Remarks on the implementation of the algorithm References [1] D. V. Chudnovsky and G. V. Chudnovsky, Computer Algebra in the Service of Mathematical Physics and Number Theory (D. V. Chudnovsky, R. D. Jenks eds.), Computers in Mathematics, Marcel Dekker, New York-Basel, 1987, 109–232. [2] T. E. Hull and R. Mathon, The Mathematical Basis and a Prototype Implementation of a New Polynomial Rootfinder with Quadratic Convergence, ACM Transactions on Mathematical Software, 22 (1996), No. 3, 261–280. [3] M. A. Jenkins and J. F. Traub, A Three-Stage Variable-Shift Iteration for Polynomial Zeros and Its Relation to Generalized Rayleigh Iteration, Numer. Math., 14 (1970), 252–263. , Principles for Testing Polynomial Zerofinding Programs, ACM Transactions on [4] Mathematical Software, 1 (1975), No. 1, 26–34. [5] D. Loewenthal Improvements on the Lehmer-Schur Root Detection Method, Journal of Computational Physics, 109 (1993), 164–168. [6] H. M¨ oller, Algorithmische Lineare Algebra, Verlag Vieweg, Braunschweig-Wiesbaden, 1997 (German). [7] , http://wwwmath1.uni-muenster.de/u/mollerh. [8] A. Sch¨ onhage, Equation Solving in Terms of Computational Complexity, Proceedings of the International Congress of Mathematicians, Berkeley, Calif., USA, 1986, 131–153. [9] P. Tur´ an, On a New Method of Analysis and Its Applications, John Wiley & Sons, New York e.a., 1978. ¨ nster, Einsteinstr. 62, D-48149 Mu ¨ nster, Math. Institute, University of Mu Germany E-mail address: [email protected]