Root Radii and Subdivision for Polynomial Root-Finding
RRoot Radii and Subdivision for PolynomialRoot-Finding
R´emi Imbach (cid:63) and Victor Y. Pan Courant Institute of Mathematical SciencesNew York University, USAEmail: [email protected]://cims.nyu.edu/~imbach/ City University of New York, USAEmail: [email protected]://comet.lehman.cuny.edu/vpan/
Abstract.
The recent subdivision algorithms for univariate polynomialComplex Root Clustering (CRC) and Real Root Isolation (RRI) approx-imate all roots in a fixed Region of Interest (RoI) and, like the algorithmof Pan (1995, 2002), achieve near optimal bit complexity for the so calledbenchmark problem. On the other hand, user’s current choice for theapproximation of all complex roots of a polynomial is the package MP-Solve, implementing Ehrlich – Aberth’s iterations. Their very fast empir-ical global convergence (right from the start) is well-known although hasno formal support. Subdivision iterations allow robust implementations,one of which is currently the user’s choice for solving the RRI problem,including the task of the approximation of all real roots. Another imple-mentation is slower than MPSolve (by several orders of magnitude) forfinding all roots but outperforms MPSolve for solving the CRC problemwhere the RoI contains only a small number of roots. We analyze and ex-tend a 2000 variant of Sch¨onhage’s algorithm of (1982), which efficientlycomputes narrow real intervals that contain the moduli of all roots, thusdefining a set of annuli covering all the complex roots. We present animplementable version of this algorithm and by using the computed setsof annuli improve in practice subdivision algorithms for CRC and RRIwhile supporting near optimal bit complexity.
Keywords:
Real root isolation, Complex root clustering, Root radii algorithm,Subdivision iterations
For a polynomial P of degree d in Z [ z ] and a Gaussian integer c ∈ G , let α ( P, c ) , . . . , α d ( P, c ) be the d non-necessarily distinct roots of P | α ( P, c ) − c |≥ | α ( P, c ) − c |≥ . . . ≥ | α d − ( P, c ) − c |≥ | α d ( P, c ) − c | (cid:63) corresponding author; a r X i v : . [ c s . S C ] F e b Imbach-Pan
For all 1 ≤ i ≤ d , write r i ( P, c ) for | α i ( P, c ) − c | ; one has r ( P, c ) ≥ r ( P, c ) ≥ . . . ≥ r d − ( P, c ) ≥ r d ( P, c )In the case where c = 0, write α i ( P ) := α i ( P,
0) and r i ( P ) := r i ( P, Root Radii Covering (RRC) ProblemGiven: a polynomial P ∈ Z [ z ] of degree d , a Gaussian integer c ∈ G , a realnumber δ > Output: d positive numbers ρ c, , . . . , ρ c,d satisfying ∀ s = 1 , . . . , d, ρ c,s δ ≤ r s ( P, c ) ≤ (1 + δ ) ρ c,s (1) ρ c, , . . . , ρ c,d of Eq. (1) for fixed c ∈ G and δ > A c of d c ≤ d disjointconcentric annuli centered at c , covering the roots of P , said to be an annuli cover of the roots of P , and are going to be used in subdivision root-finding iterations. Two Root-Finding Problems.
We count roots with multiplicity and considerdiscs D ( c, r ) := { z s.t. | z − c |≤ r } on the complex plane. For a positive δ let δ ∆and δB denote the concentric δ -dilation of a disc ∆ and a box or real line segment B . We consider two root-finding problems. First, the Complex Root Clustering (CRC) ProblemGiven: a polynomial P ∈ Z [ z ] of degree d , ε > Output: (cid:96) ≤ d couples (∆ , m ) , . . . , (∆ (cid:96) , m (cid:96) ) satisfying: – the ∆ j ’s are pairwise disjoint discs of radii ≤ ε , – ∆ j and 3∆ j contain m j > – each complex root of P is in a ∆ j .Second the Real Root Isolation (RRI) ProblemGiven: a polynomial P ∈ Z [ z ] of degree d Output: (cid:96) ≤ d couples ( B , m ) , . . . , ( B (cid:96) , m (cid:96) ) satisfying: – the B i ’s are disjoint real intervals, – each B i contains a unique real root of multiplicity m j , – each real root of P is in a B i .It is quite common for the RRI problem to pre-process P ∈ Z [ z ] in order tomake it square-free, with m j = 1 for all j , but we do not use this option. We canstate both CRC and RRI problems for P with rational coefficients and readilyreduce them to the above versions with integer coefficients by scaling.Write (cid:107) P (cid:107) := (cid:107) P (cid:107) ∞ and call the value log (cid:107) P (cid:107) the bit-size of P . oot Radii and Subdivision for Polynomial Root-Finding 3 The benchmark problem.
For the bit-complexity of the so called benchmark root-finding problem of the isolation of all roots of a square-free P ∈ Z [ z ] ofdegree d and bit-size τ the record bound of 1995 ([Pan02]) is (cid:101) O (( d + τ ) d ),near optimal for τ > d and based on divide and conquer approach. It wasreached again in 2016 ([BSSY18], [BSS + Our contributions
The complexity of subdivision root-finders is dominated atits bottleneck stages of root-counting and particularly exclusion tests, at whichcostly
Taylor’s shifts , aka the shifts of the variable, are applied. We significantlyaccelerate the root-finders for both RRI and CRC problems by means of decreas-ing the use of exclusion tests and root-counting and hence Taylor’s shifts. Weachieve this by limiting complex root-finding to the intersection of three annulicovers of the roots and real root-finding to the intersection of a single annulicover with the real line. Our progress relies on a rather simple idea: unlike theknown solution of the RRC problem for a polynomial P with any coefficients weassume that P lies near a polynomial with integer coefficients. Based on exten-sive and even tedious but rewarding analysis we proved that the bit-complexityof our computations is within (cid:101) O ( d ( d + log (cid:107) P (cid:107) )) for a disc D ( c, r ), | c |∈ O (1)and r ∈ d O ( − even if we solve the RRC problem for O (1) distinct centers.Our techniques promise to strengthen the known recipes for the initializationof various polynomial root-finding iterations; so far we have only implementedand tested this for the subdivision root-finders, for the CRC and RRI problems,in which case the gain at the initialization is amplified by saving most of invoca-tions of exclusion tests and Taylor’s shifts. For implementation of our algorithmssee the C library Ccluster ; our tests confirm the predicted significant acceler-ation of the known subdivision root-finders. We compare our implementationwith ANewDsc (see [KRS16], implementing [SM16]) and
MPSolve (see [BR14],implementing Ehrlich’s iterations) which are the current user’s choices for solv-ing respectively the RRI and the CRC problems.
Previous works
We departed from the subdivision polynomial root-finding forthe CRC and RRI problems in [BSS +
16] and [SM16], resp., and from the algo-rithms for the RRC problem in [Sch82] (see [Sch82][Cor. 14.3], and [Gou93][Algo. 2])and [Pan00]. We achieved practical progress by complementing these advancedworks with our deciding trick of restricting the input class, which made com-putation of O (1) annuli covers of the roots feasible. We rely on the customaryframework for the analysis of root-finding algorithms and cite throughout therelevant sources of our auxiliary techniques. Organization of the paper.
In Sec. 2 we describe an algorithm for solvingthe RRC problem. In Secs. 3 and 4 we present our algorithms for solving theRRI and CRC problem, respectively. In Subsec. 1.1 we introduce definitions. In https://github.com/rimbach/Ccluster Imbach-Pan
Subsec. 1.2 we briefly describe the subdivision algorithm for the CRC problemof [BSS +
16] and its adaptation to the solution of the RRI problem.
Write P := P d (cid:81) di =1 ( z − α i ( P )) = (cid:80) dj =0 P j z j Root squaring iterations.
For a positive integer (cid:96) write P [ (cid:96) ] := ( P d ) (cid:96) d (cid:89) i =1 ( z − α i ( P ) (cid:96) )and so P [0] = P , P [ (cid:96) ] = ( P [ (cid:96) − ) [1] for (cid:96) ≥
1, and | α ( P [ (cid:96) ] ) |≥ | α ( P [ (cid:96) ] ) |≥ . . . ≥ | α d − ( P [ (cid:96) ] ) |≥ | α d ( P [ (cid:96) ] ) | P [ (cid:96) ] is called the (cid:96) -th root squaring iteration of P , aka the (cid:96) -th Dandelin-Lobachevsky-Gr¨affe (DLG) iteration of P . L -bit approximations, oracle polynomials. For any number c ∈ C , we saythat (cid:101) c ∈ C is an L -bit approximation of c if (cid:107) (cid:101) c − c (cid:107)≤ − L . For a polynomial P ∈ C [ z ], we say that (cid:101) P ∈ C is an L -bit approximation of P if (cid:107) (cid:101) P − P (cid:107)≤ − L ,or equivalently if (cid:107) (cid:102) P j − P j (cid:107)≤ − L for all j .A function O P : N → C [ z ] such that O P ( L ) is an L -bit approximation of P is said to be oracle for P . Boxes, Quadri-section, Line segments, Bi-section [ a − w/ , a + w/
2] + i [ b − w/ , b + w/
2] is the box B of width w centered at c = a + i b . The disc∆( B ) := D ( c, w ) is the minimal cover of B .Partition B into four congruent boxes (children of B ), of width w/ a ± w ) + i ( b ± w ).∆( B ) := D ( c, w/
2) is a minimal disc that covers a real line segment B :=[ c − w/ , c + w/
2] of width w centered at c ∈ R .Partition the segment B into two segments (children of B ) of width w/ c ± w )). A segment B = [ c − w/ , c + w/
2] defines implicitly a boxof width w centered at c , also denoted B . [BSS +
16] describes an algorithm for solving a local version of the CRC problem:for an initial RoI B (a box) it finds clusters of roots with pairwise distance lessthan ε in a small inflation of B . Since our CRC problem is for input polynomialsin Z [ z ], one can define a RoI containing all the roots by using, for instance, theFujiwara bound (see [Fuj16]). oot Radii and Subdivision for Polynomial Root-Finding 5 Subdivision iterations
The algorithm in [BSS +
16] uses subdivision iterationsor Quad-tree algorithms (inherited from [HG69], see also [Pan00]), which con-structs a tree rooted in the RoI B whose nodes are sub-boxes of B . A node B is included only if 2 B contains a root, excluded only if it contains no root.A node is active if it is neither included nor included. At the beginning of asubdivision iteration, each active node B is tested for exclusion. Then for eachconnected component C made up of active boxes, a root-counter is applied to∆( C ) and 3∆( C ), where ∆( C ) is a minimal disc covering C . If ∆( C ) and 3∆( C )contain the same number m > C ) has radius less than ε ,then(∆( C ) , m ) is returned as a solution and the boxes of ∆( C ) are marked asincluded. Each remaining active node is quadrisected into its four active chil-dren, to which a new subdivision iteration is applied. Incorporation of Newton’siterations enables quadratic convergence toward clusters of radii ε . Solving the RRI problem
Using a root separation lower bound (e.g., of[Mig95]), one derive from [BSS +
16] a solution of the RRI problem based onsymmetry of roots of P ∈ Z [ z ] about the real axis. Let disc D ( c, r ) with c ∈ R contain m roots of P . For m = 1 the root in D ( c, r ) is real. If m ≥ r ≤ sep( P ) then D ( c, r ) contains a real root of multiplicity m , where sep( P ) isa root separation lower bound for P . For the RRI problem the RoI B is a linesegment, and the subdivision tree of B is built by means of segment bisection. The T and T ∗ tests In the algorithm of [BSS + T ∗ (∆ , P ) returns an integer k ∈ {− , , . . . , d } such that k ≥ P has k roots in ∆. A result k = − P are close to the boundary of ∆. For a given disc ∆, the exclusion test T (∆ , P ) returns 0 if T ∗ (∆ , P ) returns 0 and returns − T ∗ (∆ , P ) returns anon-zero integer. Based on Pellet’s theorem we obtain the following results. Proposition 1 (see [BSSY18], Lemmas 4 and 5)
Let B be the box centeredin c with width w . The total cost in bit operations for carrying out T ∗ (∆( B ) , P ) or T (∆( B ) , P ) is bounded by (cid:101) O ( d (log (cid:107) P (cid:107) + d log max(1 , c, w ) + L (∆ , P ))) (2) T ∗ (∆( B ) , P ) returns an integer k ∈ {− , , . . . , d } ; if k ≥ then ∆( B ) contains k roots. T (∆( B ) , P ) returns an integer k ∈ {− , } ; if k = 0 then P has no root in ∆( B ) ; if k = − then P has a root in B . Bit complexity
Prop. 1 enables one to bound the Boolean cost of exclusiontests and root-counting as well as the size of subdivision tree and hence the costof solving the benchmark problem in [BSS + + (cid:101) O ( d ( d + τ )) for the benchmark problem. Imbach-Pan
Implementations
A modified versions of [BSS +
16] for the CRC problem havebeen implemented and made public within the C library called Ccluster ; see[IPY18] and [IP20]. We refer to this implementation as Ccluster . An implemen-tation of the modified algorithm of [BSS +
16] solving the RRI problem is alsoavailable within
Ccluster , it will be called
Risolate in the sequel.
We describe and analyse an algorithm for solving the RRC problem for an oracle O P for P .Let c ∈ G and P c ( z ) := P ( c + z ), so that r s ( P, c ) = r s ( P c ) for all 1 ≤ s ≤ d .Hence the RRC problem for a c (cid:54) = 0 reduces to the RRC problem for c = 0 atthe cost of shifting the variable.The next remark reduces the RRC problem for c = 0 and any δ > δ = 4 d by means of DLG iterations: Remark 2
Let g = (cid:100) log log(4 d )log(1 + δ ) (cid:101) , ρ (cid:48) d < r s ( P c [ g ] ) < (4 d ) ρ (cid:48) and ρ = ( ρ (cid:48) ) g and recall that r s ( P c [ g ] ) = r s ( P, c ) g . Then ρ δ < r s ( P, c ) < (1 + δ ) ρ and g is in O (log d ) if δ is in d − O (1) (for instance d − or d − ). Now define the RRC ∗ problem corresponding to the RRC problem for 1+ δ =4 d and c = 0: RRC ∗ problemGiven: a polynomial P ∈ G [ z ] of degree d , satisfying P (0) (cid:54) = 0 Output: d positive real numbers ρ , . . . , ρ (cid:48) d satisfying ∀ s = 1 , . . . , d, ρ (cid:48) s d < r s ( P ) < (4 d ) ρ (cid:48) s In this setting, we assume that 0 is not a root of P and thus P (cid:54) = 0. When 0is a root of multiplicity m , then r d ( P ) = . . . = r d − m +1 ( P ) = 0 and P = . . . = P m − = 0, which is easily detected (since P ∈ G [ z ]) and treated accordingly.The rest of this section is organized as follows. In Subsec. 2.1 we present analgorithm SolveRRC ∗ that takes as an input an M -bit approximation (cid:101) P of P and solves the RRC ∗ problem. Proposition 3
For an input being a -bit approximation of P , Algorithm SolveRRC ∗ solves the RRC ∗ problem by involving O ( d log (cid:107) P (cid:107) ) bit operations. https://github.com/rimbach/Ccluster oot Radii and Subdivision for Polynomial Root-Finding 7 In Subsec. 2.2 we complete formal support for this algorithm.
Theorem 4
An algorithm
SolveRRC of Subsec. 2.3 takes as an input an oraclepolynomial O P and solves the RRC problem for δ = d − at a Boolean cost in (cid:101) O ( d log d (log d + d log( | c | +1) + log (cid:107) P (cid:107) )) This bound turns into (cid:101) O ( d ( d + log (cid:107) P (cid:107) )) for | c |∈ O (1) and into (cid:101) O ( d log d + d log d log (cid:107) P (cid:107) )) for | c | = 0 . We make two remarks concerning Thm. 4. (i) According to our analysis ofthe RRC problem, the cost of its solution for O (1) centers c such that | c |∈ O (1)is dominated by a near optimal bit-complexity of root-finding. (ii) For c = 0,our algorithm has a larger bit complexity than the algorithm of [Sch82] (see[Sch82][Cor. 14.3] and [Gou93][Algo. 2]), which is in (cid:101) O ( d log d ) when log (cid:107) P (cid:107)∈ O ( d ). Our algorithm, however, is easier to implement: Sch¨onhage rescales thevariable at each DLG iteration for each index s . ∗ problem Recall that P = (cid:80) di =0 P i z i and define, for i = 0 , . . . , d , p i = (cid:26) log | P i | if P i (cid:54) = 0 , −∞ otherwise.According to the following result (see Prop. 4.2 and its proof in [Pan00]), theRRC ∗ problem can be solved by computing the upper part of the convex hull ofthe set of points { ( i, p i ) | i = 0 , . . . , d } , provided that a point ( i, −∞ ) does not lieabove any line in the plane. Proposition 5
Given the integer s , let t (cid:48) and h (cid:48) be integers s.t. ( i (cid:48) ) t (cid:48) < d + 1 − s ≤ t (cid:48) + h (cid:48) ≤ d , and ( ii (cid:48) ) ∀ ≤ i ≤ d , the point ( i, p i ) lies below the line (( t (cid:48) , p t (cid:48) ) , ( t (cid:48) + h (cid:48) , p t (cid:48) + h (cid:48) )) Then ρ (cid:48) s = (cid:12)(cid:12)(cid:12)(cid:12) P t (cid:48) P t (cid:48) + h (cid:48) (cid:12)(cid:12)(cid:12)(cid:12) h (cid:48) satisfies: ρ (cid:48) s d < r s ( P ) < (2 d ) ρ (cid:48) s . Call CH the upper part of the convex hull of the points { ( i, p i ) | i = 0 , . . . , d } . CH can be computed exactly (the P i ’s are Gaussian integers) at the price ofexpansive integer arithmetic. Instead we use an M -bit approximation (cid:101) P of P with M ≥
2. For i = 0 , . . . , d , write (cid:101) p i := log | (cid:101) P i | if | (cid:101) P i | > − M , − − M ≤ | (cid:101) P i |≤ − M −∞ otherwise . (3)Let (cid:103) CH be the upper part of the convex hull of { ( i, (cid:101) p i ) | i = 0 , . . . , d } and letpoints ( i, −∞ ) lie below any line in the plane. Given an index s bound the slopeof the edge of CH above d + 1 − s in terms of the slope of the edge of (cid:103) CH above d + 1 − s and M . Imbach-Pan
Algorithm 1
SolveRRC ∗ ( (cid:101) P ) Input: An M -bit approximation (cid:101) P ∈ C [ z ] of P ∈ G [ z ] of degree d , s.t. P (0) (cid:54) = 0. Output: d positive real numbers (cid:101) ρ (cid:48) , . . . , (cid:101) ρ d (cid:48) .1: for i = 0 , . . . , d do
2: Compute (cid:101) p i as defined in Eq. (3).3: (cid:103) CH ← { ( i k , (cid:102) p i k ) | k = 0 , . . . , (cid:96) } , the upper part of the convex hull of { ( i, (cid:101) p i ) | i =0 , . . . , d } for k = 1 , . . . , (cid:96) do for s = d + 1 − i k , . . . , d + 1 − i k − do (cid:101) ρ s (cid:48) ← | (cid:101) P ik − (cid:101) P ik | ik − ik − //double precision floating point return (cid:101) ρ (cid:48) , . . . , (cid:101) ρ d (cid:48) Proposition 6
Let s , t, h, t (cid:48) , and h (cid:48) be integers such that ( i ) t < d + 1 − s ≤ t + h ≤ d , ( ii ) ∀ ≤ i ≤ d , the point ( i, (cid:101) p i ) lies below the line (( t, (cid:101) p t ) , ( t + h, (cid:103) p t + h ))( i (cid:48) ) t (cid:48) < d + 1 − s ≤ t (cid:48) + h (cid:48) ≤ d , and ( ii (cid:48) ) ∀ ≤ i ≤ d , the points ( i, p i ) lie below the line (( t (cid:48) , p t (cid:48) ) , ( t (cid:48) + h (cid:48) , p t (cid:48) + h (cid:48) )) .Then (cid:103) p t + h − (cid:101) p t h − − M +2 ≤ p t (cid:48) + h (cid:48) − p t (cid:48) h (cid:48) ≤ (cid:103) p t + h − (cid:101) p t h + 2 − M +2 . (4)We postpone the proof of Prop. 6. Remark that 2 − L ≤ − L for L ≥ Corollary 7 (of Prop. 6)
Let s, t, h be as in Prop. 6. Define (cid:101) ρ s (cid:48) as (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:101) P t (cid:93) P t + h (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) h .Then (cid:101) ρ s (cid:48) (2 d )(1 + 2 − M +2 ) < r s ( P ) < (2 d )(1 + 2 − M +2 ) (cid:101) ρ s (cid:48) (5)We are ready to describe our Algo. 1, which solves the RRC ∗ problem. In steps1 and 2, approximations (cid:101) p i of p i = log | P i | are computed from 2-bit approxima-tions (cid:101) P i of P i , for i = 0 , . . . , d . This requires O ( d log (cid:107) P (cid:107) ) bit operations. In step 3we compute the convex hull (cid:103) CH of a polygon with d +1 vertices (0 , (cid:101) p ) , . . . , ( d, (cid:101) p d )ordered with respect to their ordinates. Using Graham’s algorithm of [GY83],we only need O ( d ) arithmetic operations (additions) with numbers of magnitude O (log (cid:107) P (cid:107) ). In steps 4,5,6, the (cid:101) ρ s (cid:48) ’s for s = 0 , . . . , d are computed as in Cor. 7.This task is performed with double precision arithmetic with rounding and re-quires O ( d ) bit operations. Finally, if M ≥
2, (1 + 2 − M +2 ) ≤
2; thus when the M -bit approximation of P in input of Algo. 1 satisfies M ≥
2, the (cid:101) ρ s (cid:48) ’s in theoutput satisfy ∀ s = 1 , . . . , d, (cid:102) ρ s (cid:48) d < r s ( P ) < (4 d ) (cid:101) ρ s (cid:48) and Prop. 3 follows. oot Radii and Subdivision for Polynomial Root-Finding 9 Eq. (3) implies that (cid:101) p i is an ( M − p i if (cid:101) p i > −∞ ; write P i = (cid:101) P i + ε with ε ≤ − M . Then (cid:12)(cid:12)(cid:12) log | P i |− log | (cid:101) P i | (cid:12)(cid:12)(cid:12) ≤ (cid:12)(cid:12)(cid:12) log( | (cid:101) P i | + | ε | ) − log | (cid:101) P i | (cid:12)(cid:12)(cid:12) = log (cid:32) | ε || (cid:101) P i | (cid:33) ≤ | ε || (cid:101) P i | ≤ − M +1 since | (cid:101) P i |≥ (cid:103) CH − (cid:103) CHCH tt (cid:48) t (cid:48) + h (cid:48) t + hd + 1 − s (cid:103) CH + slope: (cid:94) pt + h − (cid:102) pth + 2 − M +2 (cid:102) p i − (cid:102) p i slope: (cid:94) pt + h − (cid:102) pth − − M +2 (cid:102) p i + j i j i j i j i ij i j i j Fig. 1.
The convex hulls CH and (cid:103) CH . For i = 0 , . . . , d , define (cid:101) p i + and (cid:101) p i − as (cid:101) p i ∗ = (cid:26) (cid:101) p i ∗ − M +1 if | (cid:101) P i | > −∞−∞ otherwisewhere ∗ stays above for either + or − . Recall that (cid:103) CH is the upper part ofthe convex hull of { ( i, (cid:101) p i ) | i = 0 , . . . , d } . Suppose that it is the poly-line passingthrough the points { ( i k , (cid:102) p i k ) | k = 0 , . . . , (cid:96) } . It defines two poly-lines: – (cid:103) CH + , the poly-line with vertices { ( i k , (cid:102) p i k + | k = 0 , . . . , (cid:96) } , and – (cid:103) CH − , the poly-line with vertices { ( i k , (cid:102) p i k − ) | k = 0 , . . . , (cid:96) } .Recall also that CH is the upper part of the convex hull of { ( i, p i ) | i = 0 , . . . , d } ,and suppose it is the poly-line with vertices { ( j k , p j k ) | k = 0 , . . . , (cid:96) (cid:48) } . For demon-stration see Fig. 1 where d = 8, the (cid:101) p i ’s are drawn with black circles, the intervals[ (cid:101) p i − , (cid:101) p i + ] with bold vertical bars, (cid:103) CH with a bold blue poly-line, (cid:103) CH + and (cid:103) CH − with dashed blue poly-lines, and CH with a bold red line. One has: Proposition 8
The poly-line CH lies below the poly-line (cid:103) CH + and above thepoly-line (cid:103) CH − . Proof of Prop. 8:
In order to show that CH lies below (cid:103) CH + , we show thatif j t , i k , i k (cid:48) is a triple of integers j t , i k , i k (cid:48) such that ( j t , p j t ) is a vertex of CH and [( i k , (cid:102) p i k + ) , ( i k (cid:48) , (cid:103) p i k (cid:48) + )] is an edge of (cid:103) CH + , then ( j t , p j t ) lies on or below theline (( i k , (cid:102) p i k + ) , ( i k (cid:48) , (cid:103) p i k (cid:48) + )). Suppose it is not the case, i.e. ( j t , p j t ) lies strictlyabove the line (( i k , (cid:102) p i k + ) , ( i k (cid:48) , (cid:103) p i k (cid:48) + )). Since p j t ≤ (cid:102) p j t + , (cid:102) p j t + lies strictly above(( i k , (cid:102) p i k + ) , ( i k (cid:48) , (cid:103) p i k (cid:48) + )), thus (cid:102) p j t lies strictly above (( i k , (cid:102) p i k ) , ( i k (cid:48) , (cid:103) p i k (cid:48) )) and (cid:103) CH is not the convex hull of { ( i, (cid:101) p i ) | i = 0 , . . . , d } , which is a contradiction.In order to show that (cid:103) CH − lies below CH , we show that for a given triple ofintegers i t , j k , j k (cid:48) such that ( i t , (cid:102) p i t − ) is a vertex of (cid:103) CH − and [( j k , p j k ) , ( j k (cid:48) , p j k (cid:48) )]is an edge of CH , ( i t , (cid:102) p i t − ) lies on or below the line (( j k , p j k ) , ( j k (cid:48) , p j k (cid:48) )). Supposeit is not the case. Since p i t ≥ (cid:102) p i t − , p i t lies strictly above (( j k , p j k ) , ( j k (cid:48) , p j k (cid:48) )) and CH is not the convex hull of { ( i, p i ) | i = 0 , . . . , d } , which is a contradiction. (cid:117)(cid:116) Proof of Prop. 6:
Given the integer s , let t, h, t (cid:48) , h (cid:48) be integers such thatconditions ( i ) , ( ii ) , ( i (cid:48) ) and ( ii (cid:48) ) hold.From conditions ( i (cid:48) ) and ( ii (cid:48) ), ]( t (cid:48) , p t (cid:48) ) , ( t (cid:48) + h (cid:48) , p t (cid:48) + h (cid:48) )] is the edge of CH whose orthogonal projection onto the abscissa axis contains d + 1 − s , and p t (cid:48) + h (cid:48) − p t (cid:48) h (cid:48) is the slope of that edge.From conditions ( i ) and ( ii ), ]( t, (cid:101) p t ) , ( t + h, (cid:103) p t + h )] is the edge of (cid:103) CH whoseorthogonal projection onto the abscissa axis contains d + 1 − s . Consider the twosegments ]( t, (cid:101) p t − ) , ( t + h, (cid:103) p t + h − )] and ]( t, (cid:101) p t + ) , ( t + h, (cid:103) p t + h + )] that are the edgesof (cid:103) CH − and (cid:103) CH + , respectively, whose orthogonal projections onto the abscissaaxis also contain d + 1 − s .From Prop. 8 CH is a poly-line enclosed by (cid:103) CH − and (cid:103) CH + and since thefirst coordinates of its vertices are integers, its slope p t (cid:48) + h (cid:48) − p t (cid:48) h (cid:48) above d + 1 − s is bounded below by (cid:103) p t + h − (cid:101) p t h − − M +2 and above by (cid:103) p t + h − (cid:101) p t h +2 − M +2 , whichproves Prop. 6. See Fig. 1 for an illustration. (cid:117)(cid:116) According to Rem. 2, solving the RRC problem amounts toStage 1: computing a 2-bit approximation (cid:101) P [ g ] c of P c [ g ] where g is defined as in Rem. 2Stage 2: solving the RRC ∗ problem for polynomial (cid:101) P [ g ] c Stage 3: computing ρ c,s = ( ρ (cid:48) s ) g for s = 1 , . . . , d with double precision arithmeticwith correct rounding, where the ρ (cid:48) s ’s are the output of stage 2.To estimate the cost at Stage 1, recall the following:1. (cid:107) P c [ i ] (cid:107)≤ ( d + 1)( (cid:107) P c [ i − (cid:107) ) ≤ . . . ≤ ( d + 1) i ( (cid:107) P c (cid:107) ) i
2. the loss of precision while computing P c [ i ] from P c [ i − is in O (log d +log (cid:107) P c [ i − (cid:107) )3. (cid:107) P c (cid:107)≤ (cid:107) P (cid:107) ( | c | +1) d oot Radii and Subdivision for Polynomial Root-Finding 11 From (1) and (2) above, the loss of precision while computing P c [ g ] from P c is O (( g − (cid:88) i =0 i )(log d + log (cid:107) P c (cid:107) )) = O (2 g (log d + log (cid:107) P c (cid:107) ))and since solving the RRC ∗ problem requires a 2-bit approximation of P c [ g ] , oneuses (3) to deduce that this requires an N -bit approximation of P c where N ∈ O (2 g (log d + d log( | c | +1) + log (cid:107) P (cid:107) ))Notice also that log (cid:107) P c [ g ] (cid:107) is an O ( N ). According to [KS13][Theorem 12], com-puting an N -bit approximation of P c requires (cid:101) O ( d ( d + log (cid:107) P (cid:107) + d log max(1 , | c | ) + N ))bit operations, that is, (cid:101) O ( d (log d + d log( | c | +1) + log (cid:107) P (cid:107) ))for g ∈ O (log d ). Computing an N -bit approximation of P c requires an L -bitapproximation of P where L ∈ (cid:101) O ( d + log (cid:107) P (cid:107) + d log max(1 , | c | ) + N ) = (cid:101) O ( N ).Performing any of the g DLG iterations with an initial N -bit approximationof P c involves (cid:101) O ( d ( N + log (cid:107) P c [ g ] (cid:107) ) = (cid:101) O ( dN )bit operations. Hence performing g DLG iterations involves (cid:101) O ( gdN ) = (cid:101) O ( dg g (log d + d log( | c | +1) + log (cid:107) P (cid:107) ))= (cid:101) O ( d log d (log d + d log( | c | +1) + log (cid:107) P (cid:107) ))bit operations. This dominates the cost bounds at Stages 1, 2 (which is O ( dN ))and 3 (clearly), and we complete the proof of Thm. 4.We give in Algo. 2 a procedure for solving the RRC problem, and we assumea computational framework where an analysis of the precision loss is performedwith arithmetic operations on approximated numbers. Within this framework,when computing (cid:101) P [ g ] c from O P ( L ) in Step 6, one also obtains an M s.t. (cid:107) (cid:101) P [ g ] c − P c [ g ] (cid:107)≤ − M . This can be implemented with ball arithmetic (see [Joh17]). Thewhile loop in steps 4 to 7 of Algo. 2 terminates in the worst case for an L ∈ (cid:101) O ( N ).In step 1, we use a heuristic to predict the required precision. In this section, we propose to use root distances to 0 in order to improve inpractice subdivision approaches to real root isolation, and in particular the onedescribed in Subsec. 1.2. Let us define the notion of annuli cover of the roots of P ∈ Z [ z ]. Algorithm 2
SolveRRC ( O P , c, δ ) Input: O P an oracle for P ∈ Z [ z ] of degree d , a center c ∈ G and a relative precision δ > Output: d positive real numbers ρ c,s , . . . , ρ c,d solving task S.1: L ← × (cid:100) log( d/ (cid:101)− //a heuristic for predicting the precision M ← g ← (cid:100) log log(4 d )log(1 + δ ) (cid:101) while M < do L ← L ; (cid:101) P ← O P ( L )6: Compute (cid:101) P [ g ] c with internal precision L
7: Let M be s.t. (cid:107) (cid:101) P [ g ] c − P c [ g ] (cid:107)≤ − M ρ (cid:48) , . . . , ρ (cid:48) d ← SolveRRC ∗ ( (cid:101) P [ g ] c )9: for s = 0 , . . . , d do ρ c,s ← ( ρ (cid:48) s ) / g return ρ c, , . . . , ρ c,d Definition 1.
A set A c of disjoint concentric annuli centered in c is an annulicover of the roots of P of degree d if1. ∀ A ∈ A c , there are indexes t ( A ) and h ( A ) such that α t ( A ) ( P, c ) , α t ( A )+1 ( P, c ) , . . . , α t ( A )+ h ( A ) ( P, c ) ∈ A ∀ i ∈ { , . . . , d } , there is an A ∈ A c so that α i ( P, c ) ∈ A .For an annulus A ∈ A c , r ( A ) and r ( A ) are the interior and exterior radiiof A , respectively. Write s + ( A ) := P ( c + r ( A )) × P ( c + r ( A )) and s − ( A ) := P ( c − r ( A )) × P ( c − r ( A )) . An annuli cover A centered at 0 enables us to skip many calls for exclusion androot- counting based on the following: Remark 9
Let A ∈ A such that r ( A ) > .1. if h ( A ) = 0 and s + ( A ) > (resp. s − ( A ) > ), A ∩ R + (resp. A ∩ R − )contains no real root of P .2. if h ( A ) = 0 and s + ( A ) < (resp. s − ( A ) < ), A ∩ R + (resp. A ∩ R − )contains one real root of P .3. if h ( A ) > and s + ( A ) < (resp. s − ( A ) < ), A ∩ R + (resp. A ∩ R − )contains at least one real root of P . In Subsections 3.1 and 3.2 we describe our exclusion test and root counterbased on Rem. 9. In Subsec. 3.3 we describe our algorithm solving the RRIproblem. In Subsec. 3.4 we give numerical results. oot Radii and Subdivision for Polynomial Root-Finding 13
Algorithm 3 C R ( B, P, A ) Input:
A polynomial P ∈ Z [ z ] of degree d , an interval B of R , an annuli cover A ofthe roots of P centered in 0. Assume 0 / ∈ B . Output: an integer in {− , } ; if 0 then P has no real root in B ; if − B .1: Compute n ( A , B ), n ( A , B ) and n ≥ ( A , B )2: if n ( A , B ) = n ( A , B ) then return if n ≥ ( A , B ) > = 1 then return − return T (∆( B ) , P ) Let B be a real interval that does not contain 0, let A be an annuli cover centeredin 0 and A ∈ A . We define: – s B , the sign of B , as − B <
B > – R s B ( A ) as A ∩ R − if s B < A ∩ R + otherwise – ss B ( A ) as s − ( A ) if s B < s + ( A ) otherwise – n ( A , B ): the number of annuli A ∈ A s.t. A ∩ B (cid:54) = ∅ , – n ( A , B ): the number of annuli A ∈ A s.t.( A ∩ B (cid:54) = ∅ ) ∧ ( h ( A ) = 0) ∧ ss B ( A ) > – n ≥ ( A , B ): the number of annuli A ∈ A s.t.( A ∩ B (cid:54) = ∅ ) ∧ ( h ( A ) ≥ ∧ ss B ( A ) < ∧ R s B ( A ) ⊆ B By virtue of Rem. 9, if n ( A , B ) = n ( A , B ), then all the annuli intersecting B contain no root, thus B contains no root. If n ≥ ( A , B ) ≥ B contains atleast one real root.Our exclusion test C R is described in Algo. 3. For computing n ( A , B ), n ( A , B )and n ≥ ( A , B ) in Step 1, we use double precision interval arithmetic, so thatStep 1 involves O ( d ) bit operations. As a consequence, one has Proposition 10
Let B be an interval with / ∈ B , and let A be an annuli coverof the roots of P centered in . The cost for carrying out C R ( B, P, A ) is boundedby the cost for carrying out T (∆( B ) , P ) . C R ( B, P, A ) returns an integer k in {− , } . If k = 0 , then P has no realroots in B . If k = − , then P has a root in B . In order to describe our root counter, we define:
Algorithm 4 C ∗ R ( B, P, A ) Input:
A polynomial P ∈ Z [ z ] of degree d , an interval B of R , an annuli scover A ofthe roots of P centered in 0. Assume 0 / ∈ B . Output: an integer k in {− , , , . . . , d } ; if k ≥ P has k roots in B ; if − B \ (1 / B .1: Compute n ( A , B ), n ( A , B ), n ( A , B ) and n (cid:48)≥ ( A , B )2: if n ( A , B ) = n ( A , B ) + n ( A , B ) then return n ( A , B )4: if n (cid:48)≥ ( A , B ) ≥ then return − return T ∗ (∆( B ) , P ) – n ( A , B ): the number of annuli A ∈ A s.t. :( A ∩ B (cid:54) = ∅ ) ∧ ( h ( A ) = 0) ∧ ( ss B ( A ) < ∧ ( R s B ( A ) ⊆ B ) – n (cid:48)≥ ( A , B ): the number of annuli A ∈ A s.t.:( A ∩ B (cid:54) = ∅ ) ∧ ( h ( A ) ≥ ∧ ss B ( A ) < ∧ ( R s B ( A ) ⊆ B \ (1 / B )By virtue of of Rem. 9, if n ( A , B ) = n ( A , B ) + n ( A , B ), B contains exactly n ( A , B ) real roots. If n (cid:48)≥ ( A , B ) ≥ P has at least one real root in2 B \ (1 / B .Our root counter is described in Algo. 4. For computing n ( A , B ), n ( A , B ), n ( A , B ) and n (cid:48)≥ ( A , B ) in Step 1, we use double precision interval arithmetic,thus step 1 involves O ( d ) bit operations. Proposition 11
Let B be an interval with / ∈ B and let A be an annuli coverof the roots of P centered in . The cost of carrying out C ∗ R ( B, P, A ) is boundedby the cost of carrying out T ∗ (∆( B ) , P ) . C ∗ R ( B, P, A ) returns an integer k in {− , , . . . , d } . If k ≥ , then P has k roots in ∆( B ) . If k = − , then P has a root in B \ (1 / B . Consider the procedure:Stage 1: Compute A by calling SolveRRC ( P, , d − )Stage 2: Apply the subdivision procedure of Subsec. 1.2 while using C R ( B, P, A )(resp. C ∗ R ( B, P, A )) as the exclusion test (resp. root counter) for interval B of the subdivision treeBy virtue of Thm. 4 and Propositions 10 and 11, this procedure solvesthe RRI problem and its bit-complexity is bounded by the bit-complexity ofthe algorithm described in [BSS + oot Radii and Subdivision for Polynomial Root-Finding 15 The procedure given in Subsec. 3.3 has been implemented within the library
Ccluster ; we call this implementation
RisolateR . Next we present our exper-imental results. Comparison of
RisolateR with
Risolate reveals practical im-provement due to using our root radii algorithms in subdivision process. We alsocompare
RisolateR with the subdivision algorithm of [SM16] whose implemen-tation
ANewDsc is described in [KRS16] and is currently the user’s choice for realroot isolation.
Test polynomials
We consider the following polynomials.The Bernoulli polynomial of degree d is B d ( z ) = (cid:80) dk =0 (cid:0) dk (cid:1) b d − k z k where the b i ’s are the Bernoulli numbers.The Wilkinson polynomial of degree d is W d ( z ) = (cid:81) di =1 ( z − i ).For integers n >
0, we define polynomials with (2 n + 1) × (2 n + 1) roots onthe nodes of a regular grid centered at 0 as P (2 n +1) × (2 n +1) ( z ) = (cid:89) − n ≤ a,b ≤ n ( z − a + i b )We also consider dense polynomials of degree d with coefficients randomlychosen within [ − τ − , τ − ] (under the uniform distribution). Results
In our test polynomials with several degrees/bit-sizes, we used
Risolate , RisolateR and
ANewDsc to solve the RRI problem. Our non-random exampleshave only simple roots and for those examples
ANewDsc is called with option -S1 to avoid testing input polynomial for square-freeness.Times are sequential times in seconds on a
Intel(R) Core(TM) i7-8700CPU @ 3.20GHz machine with Linux. We report in Tab. 1:- d , τ and d R , that is, the degree, the bit-size and the number of real roots,respectively- t (resp. t ), the running time of Risolate (resp.
RisolateR ),- n (resp. n ), the number of T -test in Risolate (resp.
RisolateR ),- n (cid:48) (resp. n (cid:48) ), the number of T ∗ -test in Risolate (resp.
RisolateR ),- t , the time required to compute the annuli-cover in RisolateR ,- t , the running time in second of ANewDsc .For random polynomials, we give averages over 20 examples of those values.Compare columns n , n (cid:48) and n , n (cid:48) in Tab. 1: using the annuli cover both inexclusion tests and root counter reduces dramatically the number of Pellet’s testperformed in the subdivision process, and significantly decreases the runningtime (see column t /t ). In the cases where the ratio σ/d is low most of the timein RisolateR is spent on solving the RRC problem (see column t /t ). Finally, ANewDsc remains faster than
RisolateR for polynomials having a few real rootsor a low bit-size, whereas this trend seems to reverse when the ratios of thenumber of real roots and/or bit-size over the degree increase (see columns t and t ). Risolate RisolateR ANewDsc d τ d R t n , n (cid:48) t n , n (cid:48) t /t (%) t /t (%) t
20 monic random dense polynomials per degree/bit-size256 16 5.00 .255 57.1,26.5 .047 0.0,.850 82.9 18.4 .015256 8192 6.00 3.28 130.,98.3 .454 3.55,24.6 24.2 13.8 .507256 16384 6.70 5.32 157.,130. .516 3.65,23.6 17.5 9.70 2.02256 32768 7.10 7.87 153.,132. .498 5.50,26.5 16.8 6.33 1.28391 16 4.90 .591 61.0,29.9 .132 .400,2.15 63.4 22.3 .036391 8192 7.50 9.31 158.,124. 1.06 4.40,24.7 16.1 11.4 3.08391 16384 8.80 11.2 178.,134. 1.46 4.55,28.9 15.8 12.9 5.81391 32768 8.80 24.6 200.,169. 1.77 4.80,33.2 12.0 7.17 3.77512 16 5.10 1.46 65.1,30.9 .280 .400,1.90 74.9 19.1 .077512 8192 7.10 37.5 163.,122. 2.98 5.85,26.3 13.3 7.94 2.33512 16384 7.20 48.6 156.,123. 4.72 5.45,28.3 5.79 9.72 3.68512 32768 7.00 62.3 168.,136. 2.17 4.65,25.5 16.5 3.49 5.42Bernoulli polynomials256 1057 64 1.23 292, 82 0.10 12, 3 57.0 8.35 0.21391 1810 95 2.93 460, 145 0.44 12, 2 77.5 15.0 1.14512 2591 124 6.63 528, 144 0.40 14, 3 66.0 6.06 1.62791 4435 187 16.3 892, 264 2.70 20, 1 86.5 16.5 10.61024 6139 244 55.9 1048, 283 2.46 12, 3 77.4 4.40 15.0Wilkinson polynomials256 1691 256 3.66 1030, 281 0.18 0, 10 45.7 4.93 1.62391 2816 391 17.8 1802, 539 0.70 0, 10 53.9 3.96 5.75512 3883 512 26.6 2058, 531 1.14 0, 11 49.1 4.28 27.6791 6489 791 163. 3698, 1108 6.43 0, 11 57.6 3.92 166.1024 8778 1024 260. 4114, 1047 8.02 0, 12 53.8 3.08 320.Polynomials with roots on a regular grid289 742 17 0.38 86, 30 0.13 0, 16 82.9 36.1 0.09441 1265 21 0.85 106, 36 0.20 0, 20 77.8 24.5 0.41625 1949 25 1.56 118, 39 0.91 0, 24 89.3 58.5 0.80841 2801 29 3.19 154, 51 1.71 0, 28 88.1 53.7 2.591089 3829 33 7.89 166, 55 2.09 0, 32 78.9 26.5 4.64
Table 1.
Runs of
Risolate , RisolateR and
ANewDsc on our test polynomials.
In this section, we propose to use root distances to three centers, namely 0, 1and i to improve practical performance of subdivision algorithms for complexroot clustering.Using three annuli covers A , A and A i of the roots of P , one can computea set D of O ( d ) complex disks containing all the roots of P , and use this setin the exclusion test both to filter the subdivision tree and to save expensivecomputations in Pellet’s theorem based exclusion.In Subsec. 4.1 we describe an exclusion test using the set D of discs containingthe roots of P , and in Subsec. 4.2 we present a procedure solving the CRC withnear optimal bit complexity. In Subsec. 4.3 we show experimental results. Let D be a set of O ( d ) complex discs covering all the roots of P , i.e. any rootof P is in at least one disc in D . A box B such that B ∩ D = ∅ cannot containa root of P .We define an exclusion test based on the above consideration, called C C -testand described in Algo. 5. For a box B having a nonempty intersection with thereal line, the number n ≥ ( A , B ) of annuli intersecting B and containing at leastone real root in B ∩ R is used to save some T -tests. oot Radii and Subdivision for Polynomial Root-Finding 17 Algorithm 5 C C ( B, P, D , A ) Input:
A polynomial P ∈ Z [ z ] of degree d , a box B of C , a set D of O ( d ) complexdiscs covering all the roots of P , an annuli cover A centered in 0 Output: an integer in {− , } ; if 0 then P has no real root in B ; if − B .1: Compute the number n of discs in D having nonempty intersection with B .2: if n = 0 then return if B ∩ R (cid:54) = ∅ then
5: Compute n ≥ ( A , B )6: if n ≥ ( A , B ) > = 1 then return − return T (∆( B ) , P ) Proposition 12
Let D contain O ( d ) discs covering the roots of P and let A be an annuli cover of the roots of P centered in . The cost of performing C C ( B, P, D , A ) is bounded by the cost of performing T (∆( B ) , P ) . C C ( B, P, D , A ) returns an integer k in {− , } . If k = 0 , then P has no rootin B . If k = − , then P has a root in B . Consider the procedure:Stage 1: For c = 0 , , i , compute A c by calling SolveRRC ( P, c, d − )Stage 2: Use A , A and A i to compute a set D of at most 2 d disks covering theroots of P .Stage 3: Apply the Complex Root Clustering Algorithm described in Subsec. 1.2 whileusing C C ( B, P, D , A ) instead of T (∆( B ) , P ) as the exclusion test for boxes B of the subdivision treeAccording to Thm. 4, Stage 1 involves (cid:101) O ( d ( d + log (cid:107) P (cid:107) ∞ )) bit operations.In the stage 2, D is computed as follows: using double precision floatingpoint arithmetic with correct rounding, first compute complex discs containingall possible intersections of an annulus in A with an annulus in A , and obtain aset D of at most 2 d complex discs containing all roots of P . Then, for each disc∆ in D check if ∆ and its complex conjugate ∆ have a nonempty intersectionwith at least one annulus of A i , and remove ∆ from D if it does not. This stephas cost in O ( d ).By virtue of Prop. 12, the cost of performing stage 3 is bounded by the costof performing the algorithm described in Subsec. 1.2. This procedure solves theCRC problem and supports near optimal complexity for the benchmark problem. The procedure of Subsec. 4.2 is implemented within
Ccluster ; below we callthis implementation
CclusterR and present experimental results that highlight practical improvement due to using our root radii algorithms in subdivisionprocess.We used
Ccluster and
CclusterR with input value ε = 2 − to find clustersof size at most ε . We also used MPSolve-3.2.1 , with options -as -Ga -o16 -j1 to find approximations with 16 correct digits of the roots.For our test polynomials (see 3.4) we report in Tab. 2:- d and τ respectively the degree and the bit-size,- t (resp. t ), the running time of Ccluster (resp.
CclusterR ),- n (resp. n ), the number of T -test in Ccluster (resp.
CclusterR ),- t , the time to compute the three annuli-cover in CclusterR ,- t , the running time in second of MPSolve .For random polynomials, we give averages over 20 examples of those values.As for the real root isolator presented in Sec. 3, using root radii allows to decreasesignificantly the number of exclusion tests based on Pellet’s theorem (comparecolumns n and n ) in the subdivision process, and yields a speed-up factor about3 (see column t /t ) for our examples. This speed-up increases as the numberof real roots increases (see, e.g., Wilkinson polynomials) because exclusion testsfor boxes B containing the real line are avoided by using the number n ≥ ( A , B )in the C C test. The time spent to compute the three annuli covers remains lowcompared to the running time of CclusterR (see column t /t ). MPSolve remainsthe user’s choice for approximating all complex roots.
We revisited and advanced the variant of 2000 of Sch¨onhage’s efficient approx-imation of polynomial root radii of 1982. We supported our simple but novelidea with extensive and even tedious analysis. This rewarded us with significantupgrading of the initialization of the known iterative polynomial root-finders.So far we implemented our progress only in the case of subdivison root-finding,for which the initialization gain is amplified at the subsequent stages. Our testsshow significant progress vs. the known algorithms.
References
BR14. Dario A Bini and Leonardo Robol. Solving secular and polynomial equa-tions: A multiprecision algorithm.
Journal of Computational and AppliedMathematics , 272:276–292, 2014.BSS +
16. Ruben Becker, Michael Sagraloff, Vikram Sharma, Juan Xu, and Chee Yap.Complexity analysis of root clustering for a complex polynomial. In
Pro-ceedings of the ACM on International Symposium on Symbolic and AlgebraicComputation , ISSAC ’16, pages 71–78, 2016.BSSY18. Ruben Becker, Michael Sagraloff, Vikram Sharma, and Chee Yap. A near-optimal subdivision algorithm for complex root isolation based on Pellet testand Newton iteration.
Journal of Symbolic Computation , 86:51–96, May-June2018.oot Radii and Subdivision for Polynomial Root-Finding 19
Ccluster CclusterR MPSolve d τ t n t n t /t (%) t /t (%) t
20 monic random dense polynomials per degree128 128 4.24 2655. 1.35 467. 8.11 31.9 .030191 191 11.9 3778. 3.96 690. 4.33 33.2 .061256 256 24.4 4990. 8.16 914. 7.21 33.4 .103391 391 74.0 7628. 23.3 1451. 3.67 31.5 .227512 512 149. 9919. 46.2 1852. 7.02 30.9 .394Bernoulli polynomials128 411 3.86 2954 1.25 548 7.48 32.3 0.07191 690 12.2 4026 4.51 942 8.07 36.8 0.16256 1057 24.7 5950 10.1 1253 6.57 41.1 0.39391 1810 75.1 8322 27.4 1907 16.2 36.5 0.97512 2591 133. 11738 49.9 2645 12.7 37.5 2.32Wilkinson polynomials128 722 8.43 3786 1.09 14 14.4 12.9 0.17191 1184 25.4 5916 2.99 18 27.9 11.7 0.51256 1691 50.7 7500 6.34 18 21.7 12.4 1.17391 2816 201. 12780 23.1 22 36.2 11.4 4.30512 3883 379. 14994 51.3 22 35.6 13.5 9.33Polynomials with roots on a regular grid169 370 7.37 3072 1.99 592 4.03 27.1 0.05289 742 27.1 5864 10.2 1573 3.18 37.9 0.13441 1265 81.4 9976 24.4 1713 4.28 29.9 0.56625 1949 228. 15560 70.2 2508 15.0 30.7 1.16841 2801 493. 19664 169. 4294 5.75 34.2 3.84
Table 2.
Runs of
Ccluster , CclusterR and
MPSolve on our test polynomials.Fuj16. Matsusaburˆo Fujiwara. ¨Uber die obere schranke des absoluten betrages derwurzeln einer algebraischen gleichung.
Tohoku Mathematical Journal, FirstSeries , 10:167–171, 1916.Gou93. Xavier Gourdon. Algorithmique du theoreme fondamental de l’algebre. Re-search Report RR-1852, INRIA, 1993.GY83. Ronald L Graham and F Frances Yao. Finding the convex hull of a simplepolygon.
Journal of Algorithms , 4(4):324–331, 1983.HG69. Peter Henrici and Irene Gargantini. Uniformly convergent algorithms forthe simultaneous approximation of all zeros of a polynomial. In
Construc-tive Aspects of the Fundamental Theorem of Algebra , pages 77–113. Wiley-Interscience New York, 1969.IP20. R´emi Imbach and Victor Y. Pan. New practical advances in polynomial rootclustering. In Daniel Slamanig, Elias Tsigaridas, and Zafeirakis Zafeirakopou-los, editors,
Mathematical Aspects of Computer and Information Sciences ,pages 122–137, Cham, 2020. Springer International Publishing.IPY18. R´emi Imbach, Victor Y. Pan, and Chee Yap. Implementation of a near-optimal complex root clustering algorithm. In
Mathematical Software – ICMS2018 , pages 235–244, 2018.Joh17. F. Johansson. Arb: efficient arbitrary-precision midpoint-radius intervalarithmetic.
IEEE Transactions on Computers , 66:1281–1292, 2017.KRS16. Alexander Kobel, Fabrice Rouillier, and Michael Sagraloff. Computing realroots of real polynomials ... and now for real! In
Proceedings of the ACM onInternational Symposium on Symbolic and Algebraic Computation , ISSAC’16, pages 303–310, 2016.KS13. Alexander Kobel and Michael Sagraloff. Fast approximate polynomial mul-tipoint evaluation and applications. arXiv preprint arXiv:1304.8069 , 2013.0 Imbach-PanMig95. Maurice Mignotte. On the distance between the roots of a polynomial.
Appli-cable Algebra in Engineering, Communication and Computing , 6(6):327–332,Nov 1995.Pan00. Victor Y Pan. Approximating complex polynomial zeros: modified weyl’squadtree construction and improved newton’s iteration.
J. of Complexity ,16(1):213–264, 2000.Pan02. Victor Y Pan. Univariate polynomials: nearly optimal algorithms for nu-merical factorization and root-finding.
J. of Symb. Comp. , 33(5):701–733,2002.Sch82. Arnold Sch¨onhage. The fundamental theorem of algebra in terms of compu-tational complexity.
Manuscript. Univ. of T¨ubingen, Germany , 1982.SM16. Michael Sagraloff and Kurt Mehlhorn. Computing real roots of real polyno-mials.