An adaptive subdivision method for root finding of univariate polynomials

An adaptive subdivision method for root finding of univariate polynomials

Journal of Computational and Applied Mathematics 352 (2019) 146–164 Contents lists available at ScienceDirect Journal of Computational and Applied M...

759KB Sizes 0 Downloads 34 Views

Journal of Computational and Applied Mathematics 352 (2019) 146–164

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

An adaptive subdivision method for root finding of univariate polynomials ∗

Juan Luis García-Zapata a , , Juan Carlos Díaz Martín b , Álvaro Cortés Fácila b a

Departamento de Matemáticas, Escuela Politécnica, Universidad de Extremadura, Avda. de la Universidad s/n, 10071, Cáceres, Spain Departamento de Tecnología de los Computadores, Escuela Politénica, Universidad de Extremadura, Avda. de la Universidad s/n, 10071 Cáceres, Spain

b

article

info

Article history: Received 28 June 2016 Received in revised form 17 May 2018 Keywords: Root finding Subdivision Winding number Binary search tree Contour integration

a b s t r a c t The most common methods for root finding of polynomials are iterative. This type of methods shows some drawbacks, like the impossibility of parallelization, the unpredictable nature of the iterations, that prevents to search roots inside a pre-specified region of complex plane, and the moderated degree of the polynomials that they can handle. This work describes a recursive root finding method. It is based on the winding number of plane curves, that is related to the number of zeros of a polynomial in a plane region. By our previous work (García Zapata and Díaz Martín, 2014) we find the winding number at reasonable computational cost, so we can approximate the roots by recursive subdivision of the search region. This subdivision approach is parallelizable, with a known bound of the computational cost, and the search can be restrained, in contrast to the iterative methods. We use error returns to identify and avoid singular curves, adapting the subdivision to the roots. This allows us to handle high degree polynomials without resorting to the expensive multiprecision arithmetics of other recursive methods. A significant contribution is that we formally prove its correctness. © 2018 Elsevier B.V. All rights reserved.

1. Introduction: Recursive root finding We focus on the problem of finding roots (real or complex) of univariate polynomials (with coefficients real or complex). The most common methods to solve this problem are iterative. These methods are based on a sequence of error estimation and correction that, in well-defined conditions, leads to a complex point as close to a root as needed. The methods of this type are fast (better than linear convergence, in general) and its analysis (that is, the proof of its correction, and the determination of the resources needed) relies on relatively standard numerical techniques [1]. They are the common approach to the polynomial root finding problem, named after Newton, Müller, Legendre, etc. For example, the procedure built in mathematical packages (as LAPACK [2] or MATLAB [3]) is an iterative method which is based on the eigenvectors of the companion matrix of the polynomial (see [4]). An iterative method usually converges rapidly for most polynomials of moderate degree. However, even for these degrees, the iterative methods are not well suited for certain classes of polynomials (as those with clusters or with multiple roots, depending on the method), or for specific polynomials that show ill-conditioning (like the Wilkinson example [5]). To cope with these issues, it is necessary to apply heuristic rules in case of difficulties (to choose another initial approximation, or change the method applied, or use high precision arithmetic [6]). In practice, most applications of root finding rely on iterative methods, with a proper set of heuristics. Among recent developments along this line we can highlight [7] or [8]. ∗ Corresponding author. E-mail addresses: [email protected] (J.L. García-Zapata), [email protected] (J.C. Díaz Martín), [email protected] (Á. Cortés Fácila). https://doi.org/10.1016/j.cam.2018.11.023 0377-0427/© 2018 Elsevier B.V. All rights reserved.

J.L. García-Zapata, J.C. Díaz Martín and Á. Cortés Fácila / Journal of Computational and Applied Mathematics 352 (2019) 146–164

147

Fig. 1. On the left, rectangular curves generated in the recursive procedure RDP for a polynomial with degree 322. On the right, the 322 roots of the polynomial.

The iterative methods, inherently sequential, make difficult to take advantage of parallel High Performance Computing platforms. And there is a demand for root finding tools, because of the growing order of the models used in mathematical modeling. In contrast, others methods of root finding are not based on point estimations, but instead on the distribution of roots in the complex plane, and if they are contained inside a given circle. Traditionally, roots separation is applied in a first stage of root finding, as ‘‘rules of thumb’’ to separate intervals on the real line (or regions on the complex plane) each containing a desired root. On a second stage, an iterative method approximates the root in every interval, supposedly small enough for its convergence. The stage of separation of roots can be repeatedly applied, as for example in the bisection method for zero finding of real continuous functions. It is based on the Bolzano’s theorem. This theorem says that if a continuous function f : [a, b] → R in an interval [a, b] of the real straight line takes opposed signs at its ends (that is, f (a) · f (b) < 0), then it has at least a root in the interval. The bisection method consists of using this fact recursively to build a succession of nested intervals (and hence progressively smaller) containing a real root. These are the prototype components of a subdivision method: an inclusion test, to decide whether there is any root in a piece (either of the line or of the plane); and a recursive procedure to get progressively smaller localizations of desired roots. The subdivision methods are more difficult to implement than iterative methods, requiring data types for geometric objects and tree-search or backtracking procedures for control flow. Notwithstanding, unlike iterative methods, those based on subdivision are valid for all polynomials equally. This uniformity allows an analysis of the complexity of such methods. The theoretical studies of complexity have been a driving force in the development of subdivision methods [6]. A recent effort to incorporate these theoretical advances in implementations of practical use is [9]. In this work we expose a subdivision method based on an inclusion test, Insertion Procedure with control of Singularity for polynomial Roots (IPSR), and on a recursive procedure giving a binary search tree, Recursive Decomposition Procedure (RDP). The inclusion test uses the winding number (that is, the number of twists of a closed plane curve around the origin). We have shown [10] that IPSR is an efficient procedure to compute the winding number of a nonsingular curve (i.e., that does not cross the origin, hence having a definite winding number). Besides, we gave a bound of its computational cost. Knowing this bound we can early detect if the border of a region that arises by subdivision is singular or near to singular, and interrupt the computation. In this case, the procedure RDP modifies the subdivision, up to achieving subregions with nonsingular borders. The principles under RDP have not just theoretical interest. RDP has been implemented as a C library called Contour, with good performance even on high-degree ill-conditioned cases. Contour can be accessed at the website https://gim. unex.es/contour/. An example of use is Fig. 1. It shows the actual resultant subregions, of rectangular border, for an instance polynomial with degree 322 submitted to RDP. About the structure of the paper, Section 2 summarizes the results of [10] about the winding number computation of the curves that appear in the root finding method. Section 3 describes the search tree procedure, RDP. The recursive partition is done by avoiding passing near a root, because this causes a high cost of IPSR. Section 4 proves that RDP is properly defined, that is, that the procedure effectively divides each region into two smaller subregions. Finally, Section 5 contains a proof of the correctness of the RDP root finding procedure. 2. Winding number as inclusion test This section exposes a procedure to compute the winding number of a closed curve ∆ (of the form ∆ = f (Γ ) for a polynomial f and a closed curve Γ ). We have described it in [11], and we use it to find the roots of a polynomial inside a

148

J.L. García-Zapata, J.C. Díaz Martín and Á. Cortés Fácila / Journal of Computational and Applied Mathematics 352 (2019) 146–164

Fig. 2. The number of roots of the polynomial f (z) = z 3 + 1 inside Γ1 and Γ2 equals the winding numbers of ∆1 = f (Γ1 ) and ∆2 = f (Γ2 ).

closed curve (Insertion Procedure with control of Singularity for polynomial Roots, IPSR). The main properties of IPSR are summarized in Theorem 1 given later in this section. Let us consider closed curves defined parametrically, that is, as mappings of an interval to the complex plane ∆ : [a, b] → C with ∆(a) = ∆(b). The winding number, or index, Ind(∆) of a curve ∆ : [a, b] → C is the number of complete rotations, in counterclockwise sense, of the curve around the point (0, 0). As a particular case of the Cauchy formula of Complex Analysis, the winding number is equal to this line integral: Ind(∆) =

1 2π i



1 ∆

w

dw

The principle of the argument [12] states that the number of zeros (counted with multiplicity) of an analytic function f : C → C, w = f (z), inside a region with border defined by the curve Γ , is equal to the winding number of the curve ∆ = f (Γ ). The principle of the argument can be viewed as a bi-dimensional analogy of Bolzano’s theorem, and it is in the base of several recursive methods to find the zeros of holomorphic functions [13]. We apply the principle of the argument to analytic functions defined by polynomials. The number of zeros of a polynomial complex function f (z) inside a region bounded by a curve Γ is the winding number of the curve ∆ = f (Γ ) (Fig. 2). This can be used as an inclusion test to decide if a given region of the plane of border Γ has any root, that is, Ind(f (Γ )) > 0. The input parameters of the inclusion test IPSR are the curve Γ : [a, b] −→ C, the polynomial f and the parameter θ . We denote a call with IPSR(Γ , θ ), assuming that the other parameter f is clear from the context: the polynomial subject to root-finding. The accuracy and cost of IPSR are described by the following Theorem 1 [10]. Theorem 1. If Γ : [a, b] → C is a closed curve uniformly parameterized, θ a positive real number, and f a polynomial of degree n, then the procedure IPSR verifies the following: a + 1 iterations. (a) IPSR(Γ , θ ) returns in less than b− θ (b) If it exits normally, then it gives the number of roots of f inside Γ . √ (c) If it exits with error, then it gives a point γ of Γ such that there is a root z of f with d(z , γ ) ≤ E(θ ) = (4 + 2 2)nθ .





To summarize, IPSR prevents an excessive number of iterations, controlled by input parameter θ . It effectively computes the winding number of f (Γ ) on normal return. If it returns with error, ⌊ itagives⌋a point on the curve near a root and does not compute the winding number. In any case, IPSR returns in less than b− + 1 iterations. θ The bound of the cost of part (a) is necessary to estimate the resources required when we use IPSR as an inclusion test. Besides, part (c) is crucial, as it detects if we are handling a subregion whose border requires a expensive computation. Hence IPSR can identify and avoid such high cost curves, that are precisely the near to singular ones. The relationship between distance from roots and computational cost is reported in the literature ([14] and references therein), but the ability of IPRS to avoid expensive curves is not present in the inclusion test of other subdivision methods. According to Theorem 1(c), we say that a root z of f is near the curve Γ if d(z , Γ ) ≤ E(θ ). Thus, IPSR returns with error if there is (at least) a root of f near Γ . In that case, we say that the root (or roots) have caused the error. 3. Recursive decomposition of the search region In this section we detail the subdivision algorithm described in the introduction, which is based on the IPSR inclusion test to decide if there is one or more roots in region and on the recursive procedure RDP to subdivide the region and locate

J.L. García-Zapata, J.C. Díaz Martín and Á. Cortés Fácila / Journal of Computational and Applied Mathematics 352 (2019) 146–164

149

Fig. 3. Supporting lines and diameters for a plane region P.

the roots up to a given accuracy. With the result of Theorem 1, IPSR computes the number of roots inside a region P with boundary Γ as the winding number of ∆ = f (Γ ). This method will be used as inclusion test. For the recursive subdivision, in the next Section 3.1 we define precisely how the subregions are obtained. They are created by a horizontal or vertical straight cut done to the starting region. In the second Section 3.2 we describe how to achieve subregions in which IPSR does not return with error. By Theorem 1(c), it is sufficient that the straight cut is sufficiently separated from the roots. We will see that this can be achieved by shifting the cut appropriately. Finally, in the third Section 3.3 we define the recursive partition procedure, RPD, integrating the above. Before this, a comment about the type of regions P that we consider. For simplicity, the initial region PI , as well as those that arise in its subdivisions, must be arc-connected. With this assumption, the border is a simple closed curve Γ (Jordan curve [15]). A non-arc-connected region is, for instance, the interior of two disjoint circles. Non-arc-connected regions, whose border is made of several closed curves, can be treated in a way similar to the arc-connected ones, because the number of roots inside a non-arc-connected region is the sum of the winding numbers of its arc-connected components. However, the assumption of arc-connectedness simplifies the reasonings. The subregions will be produced by cuts along straight lines. To prevent such cuts to cause a non-arc-connected subregion, we impose that the initial region PI be convex (i. e., that if a straight segment has its endpoints in PI , all the segment is inside PI ). So the subregions will be also convex (hence arc-connected) because they are produced by cuts along straight lines. From the practical point of view, when finding the roots in a non-convex region, RDP should be applied to its convex hull, or to its decomposition into convex parts. 3.1. Subregions by straight cuts Now we define the horizontal or vertical lines along which we will make the cuts. We consider each plane region as a closed set, that is, including its border. For a plane region P, its top horizontal supporting line lT is the upper horizontal straight line that has some point in common with P. Likewise, the bottom horizontal supporting line lB of P is the lower horizontal straight line that has some point in common with P. The vertical diameter dmV (P) is the distance between these lines d(lT , lB ). That is dmV (P) = minx∈lT , y∈lB d(x, y). In a similar manner, the left vertical supporting line lL is the leftmost vertical straight line that has some point in common with P, and the right vertical supporting line lR is the rightmost straight line that has some point in common with P. The horizontal diameter dmH (P) is the distance between these lines d(lL , lR ). That is dmH (P) = minx∈lL , y∈lR d(x, y). Fig. 3 depicts these concepts, √ and besides the classical diameter of P, defined by dm(P) = maxx,y∈P d(x, y), and the rectangular diameter dmR (P) = dmH (P)2 + dmV (P)2 .

150

J.L. García-Zapata, J.C. Díaz Martín and Á. Cortés Fácila / Journal of Computational and Applied Mathematics 352 (2019) 146–164

Fig. 4. The action of the operators Tλ , Bλ , Lλ′ and Rλ′ .

Fig. 5. Iterated Shifted Cuts (ISC): being E(θ ) given in Theorem 1(c), the lines mH (P), mH (P) + 2E(θ ), mH (P) − 2E(θ ), mH (P) + 4E(θ ) successively cause error by proximity of the roots z1 , z2 , z3 and z4 respectively. Finally mH (P) − 4E(θ ) does not have a root at distance lesser than E(θ ).

In the case of a region P of border Γ and a real value λ, we define the operator Tλ in the following way (see Fig. 4). We denote by mH (P) the horizontal line equally distant from lT and lB . The horizontal line at distance λ above mH (P) is denoted mH +λ (P). Tλ (P) is the intersection of P with the upper half-plane defined by mH +λ (P). If λ is a negative value, the line mH +λ (P) should be understood below mH (P). Analogously the operator Bλ (P) is the intersection of P with the lower half-plane defined by the same horizontal line, mH +λ (P). In this way Tλ (P) ∪ Bλ (P) = P and Tλ (P) ∩ Bλ (P) is its common border, a segment of mH +λ (P). Likewise, the vertical line equally distant from lL and lR is mV (P), the vertical line at distance λ to the right of mV (P) is denoted mV +λ (P), and Rλ (P) is the intersection of P with the right half-plane defined by mV +λ (P), and Lλ (P) the intersection with the left half-plane defined by the same line mV +λ (P). We say that Tλ and Bλ are horizontal cut operators, and that Lλ and Rλ are vertical cut operators. We say that P has horizontal minor axis if dmV (P) ≥ dmH (P), and that it has vertical minor axis if dmV (P) < dmH (P). 3.2. Iterated shifted cuts In this subsection we define the procedure Iterated Shifted Cuts (ISC) to divide a region P in two subregions on which IPSR returns successfully. To get a decreasing in the diameter (see Proposition 3(b)), we will divide P by a straight line along its minor axis, either horizontal (if dmV (P) > dmH (P)) or vertical (if dmV (P) < dmH (P)). However, the cut will not be, in general, along mH (P) (respectively, mV (P)), but along a parallel line mH +λ (P) (respectively, mV +λ (P)), for a value λ obtained by an iterative computation, with cuts progressively shifted as long as the algorithm produces subregions with error in IPSR. See Fig. 5 for several steps of the algorithm.

J.L. García-Zapata, J.C. Díaz Martín and Á. Cortés Fácila / Journal of Computational and Applied Mathematics 352 (2019) 146–164

151

We now describe the loop of shifting cuts. The loop begins at iteration/step i = 0 with a central cut along mH (P), for horizontal cuts, or mV (P), for vertical cuts, giving rise to two subregions P0 , with border Γ0 , and P1 , with border Γ1 . Then we consider the calls IPSR(Γ0 , θ ) and IPSR(Γ1 , θ ), with a parameter θ defined in the forthcoming Proposition 1. If any of the two calls return with an error, we discard the cut and continue for the next step. In general, if the step i − 1 produces a subregion such that the IPSR return with an error, we consider a cut along a line shifted by a distance hi such that hi = (i + 1)E(θ ), if i is odd, or hi = −iE(θ ), if i is even, for E(θ ) given in Theorem 1(c). This procedure is repeated until both calls of IPSR return without error. We denote by ISC (P , i) the pair of subregions produced from P with a cut along a line shifted a distance hi from the central cut, specifying horizontal or vertical cut if necessary. We name shifting loop the sequence of these ISC (P , i) calls with increasing values of the step i. The value θ used as argument in IPSR, by Theorem 1(a), is inversely proportional to the computational cost. Hence we are interested in a value of θ as high as possible in the calls IPSR(Γ , θ ). We need also that θ (and hence E and hi ) be small enough so the last cut made by ISC in P divides it in two non empty subregions. By the forthcoming Proposition 1, an upper min(dmH (P),dmV (P)) √ bound for θ satisfying this property is given by ZΓ = . 4(4+2 2)(n+2)n

3.3. Final definition of RDP We now properly define the recursive procedure RDP to subdivide the region P and locate the roots up to a given accuracy A. This procedure is shown in the pseudocode given in Fig. 6. To better discuss it, we cite the nested calls of RDP saying that the first call has recursion depth 0, and that the calls to RDP made inside a call of recursion depth v have recursion depth v + 1. The following four points define the strategy of RDP for error management: first, when RDP(P , A) at depth v receives an error from IPSR(Γ , ZΓ ) (line 2), identifies Γ as singular and returns with error to the calling instance, at depth v − 1. RDP returns at line 6, and the following line manages the errors which may occur. Being Υ the cut segment and γ the point returned on error by RDP(P0 , A) (or RDP(P1 , A)), in line 2, we say that this is a cut error if d(γ , Υ ) ≤ E(θ ), with E(θ ) defined in Theorem 1(c), and that it is a deferred error in contrary case. Second, if the error is of cut type (line 8), this RDP instance at depth v − 1 considers the child P0 as singular (or P1 , according to the subregion that causes error), and generates a new cut in P. Third, if the error is of deferred type, the RDP instance considers its parent region P as singular and raises the error to its caller at depth v − 2, which in turn classifies it as cut or deferred, and proceeds accordingly. Fourth, the initial call RDP(PI , A) classifies all errors as cut if there are no roots near the border1 ΓI . We introduce the fourth point because if there is a root near the initial border ΓI , then it can produce an error from which we cannot recover. So RDP is a method that finds the roots in PI with the precondition that its border ΓI is separated from roots. In contrary case, RDP returns a root near ΓI , as the claim on the error output of Fig. 6 says. Fig. 7 shows a trace of the procedure. Beginning with RDP(Γa , A), the rectangular region of border Γa in step 1 is divided by the cut Υ along the minor axis in the subregions P0 and P1 , of borders Γ0 and Γ1 . In step 2 and 3 (calls RDP(Γ0 , A) and RDP(Γ1 , A), respectively), IPSR computes without error the number of roots inside them. Γ0 does not have any root and Γ1 has one, r1 . Although r1 is near the border of Γ1 , in this example it does not cause an IPSR error. Note that Theorem 1(c) asserts that if there is an error, then there is a close root, but the reciprocal is not true in general. In step 3, Γ1 is divided by Υ2 in Γ10 and Γ11 (this last subregion is not shown). Inside RDP(Γ10 , A) (step 4), IPSR(Γ10 , ZΓ10 ) returns with error the point γ = Γ10 (u). RDP(Γ10 , A) returns (step 5) with error this point to its caller (that is RDP(Γ1 , A)), which classifies it as deferred (step 6) because γ is not near to Υ2 . The error is hence raised to RDP(Γa , A), which classifies it as cut error (step 7) because γ is near to Υ . As a result, RDP(Γa , A) shifts the cut Υ (step 8), making the new subregions of borders Γ0′ and Γ1′ , that are later processed by RDP (steps 9 and 10). 4. Properties of the recursive procedure In this section - Proposition 3, Section 4.3 -, we prove that ISC does not reach a shift so large to produce an empty subregion (a subregion without error in IPSR) and that the cut along minor axis decreases the diameter of the subregion. Before that we consider two preliminary results: Proposition 1, in Section 4.1, that deals with the iterations of the shifting loop (lines 5–9 of the procedure RDP on Fig. 6) for partitioning a region P, and Proposition 2, in Section 4.2, that deals with cuts made in several instances of RDP. 4.1. Effectiveness of ISC In the shifting loop of Fig. 5, each root zi is the cause of a new iteration. This allows us to bound the number of iterations of ISC by the number of roots inside P. What is the maximum distance reachable by this bounded number of shifts? In the next Proposition 1 we give the bound of this maximum distance. The bound is a function of θ , the parameter of the calls to IPSR. dmV (P) √ Being nP the number of roots inside P, and n the degree of the polynomial, we define ΘH = and

ΘV =

(8+4 2)(nP +1)n

dmH (P)



(8+4 2)(nP +1)n

. The values hi are the shift distances previously defined (hi = (i + 1)E(θ ) if i is odd, else hi = −iE(θ ),



with E(θ ) = (4 + 2 2)nθ ). 1 This fourth point is consider to avoid the RDP algorithm to produce an error that we cannot recover if there is a root near the initial border.

152

J.L. García-Zapata, J.C. Díaz Martín and Á. Cortés Fácila / Journal of Computational and Applied Mathematics 352 (2019) 146–164

Fig. 6. Recursive Division Procedure with error handling.

Fig. 7. Successive divisions of Γa . The numeric labeling step 1, step 2, . . . , step 10 shows the temporal sequence of actions, described in the text.

J.L. García-Zapata, J.C. Díaz Martín and Á. Cortés Fácila / Journal of Computational and Applied Mathematics 352 (2019) 146–164

Proposition 1. Let P be a region containing nP roots of a polynomial of degree n, ΘH =

dmV (P)



(8+4 2)(nP +1)n

and ΘV =

153 dmH (P)



(8+4 2)(nP +1)n



,

and hi the shift distances previously defined (hi = (i + 1)E in i is odd, hi = −iE, in i is even, with E(θ ) = (4 + 2 2)nθ ). Let us suppose that there are no roots at a distance less or equal than E(θ ) from the border of P. If θ ≤ ΘH (respectively θ < ΘV ) there is an iteration imax of the shifting loop with horizontal (respectively, vertical) cuts such that the call ISC (P , imax ) returns two nonempty subregions with borders ΓT , ΓB (respectively, ΓL , ΓR ), IPSR(ΓT , θ ) and IPSR(ΓB , θ ) (respectively, IPSR(ΓL , θ ) and IPSR(ΓR , θ )) return successfully and that |himax | < (nP + 1)E(θ ). Proof. Let us consider the first iteration ISC(P , 0), the step 0 of a horizontal cut (see Fig. 5). The proof for vertical cuts is similar and will be omitted. It produces two subregions T (P) and B(P) of respective borders ΓT and ΓB , composed from the border Γ of P, and the horizontal cut Υ following mH (P). If IPSR(ΓT , θ ) returns with error, as ΓT is Υ concatenated with Γ ∪ T (P) and by hypothesis there are no roots near Γ (neither near its portion Γ ∪ T (P)) then by Theorem 1c there is a point in Υ at distance E(θ ) or less from a root. In that case, with the segment Υ of mH (P) near a root, say z1 , in the next iteration the procedure ISC(P , 1) proceeds by cutting at 2E(θ ) above mH (P). The new cut is at distance greater than E(θ ) from z1 , hence z1 cannot cause error. Perhaps there is another root z2 near this first shifted line causing again an error in IPSR. In such case ISC(P , 2) performs a horizontal cut by a second shifted line, below mH (P), at distance −2E(θ ). That is, step 1 in ISC is triggered by z1 , step 2 by z2 , and so on. Therefore ISC will perform as much up to step nP (being nP the number of roots inside P). Hence the last horizontal cutting line is the imax -th, for some value verifying imax ≤ nP . This last cutting line is at a distance himax from mH (P), that is himax = (imax + 1)E(θ ) √ if imax is odd, else himax = −imax E(θ ). In both cases |himax | ≤ (imax + 1)E(θ ), that is lesser or equal than (nP + 1)E(θ ) = (4 + 2 2)(nP + 1)nθ , as we wanted to show. It remains to show that the subregions are both non empty, in other words that the last shifted line actually cuts P. This line dmV (P) dmV (P) √ is at a distance himax from the central line, so to actually cut P means |himax | < . By hypothesis, θ < ΘH = , 2 2(4+2 2)(nP +1)n



hence |himax | ≤ (nP + 1)E(θ ) <

(4+2 2)(nP +1)n dmV (P)



2(4+2 2)(nP +1)n

=

dmV (P) , 2

and we conclude the proof. □

4.2. Decreasing diameters In this section we will show that the successive subdivisions of recursive RDP calls end up in subregions of diameter less than a given accuracy A. To this end, we introduce the notions of cut tolerance and aspect ratio of a region. We consider the cut operators Tλ , Bλ , Lλ and Rλ , generically denoted C . We will consider the subregion that arises after applying several cuts operators (or simply cuts) to a region P, denoting it as a chain Cm · · · C2 C1 (P) of cuts along minor axis. For such a chain of m cuts, let us suppose that it has a number tH of operators of horizontal type and a number tV of operators of vertical type, in such a way that m = tH + tV . We define t0 = min(tH , tV ), that is, the minimum of the numbers of operators of each type. The reduction in rectangular diameter is related to the length m of the cut chain (that is, the recursion depth m reached by RDP) through the minimum number t0 . log

dmR (P)

From Lemma 6 (in Appendix A) we deduce that, given a tolerance A, after a chain with t0 ≥ log2 (1/AK ) operators, the 2 rectangular diameter of the resulting region is less or equal than A. We will now introduce the notion of tolerance that is related with the impact of shifting on the diameter reduction. Let us denote by Hλ a generic horizontal operator, that is, either Tλ or Bλ . Similarly Vλ denotes either Lλ or Rλ . For a horizontal 2|λ| operator applied to a region P, that is Hλ (P), we define its tolerance as dm (P) . It depends on λ and P, in such a way that V operators with zero tolerance make cuts through the middle line mH (P), with tolerance 1/2 cut halfway between the middle line and one of the two horizontal supporting lines, and with tolerance 1 or greater the cut is degenerated (that is, one of 2|λ| the subregions that arise has null area or is empty). Likewise, we say that a vertical operator Vλ (P) has a tolerance of dm (P) . H ′ In Fig. 4, both cuts have tolerance 1/5, although λ > λ . We have chosen this term from the field of mechanical engineering, referring to the permissible range of variation in a dimension of an object without harming its function. We consider that a cut along the middle line is at the reference point, and a cut with some shifting moves away from this reference. If Ci , for i = 1, . . . m, denotes an operator either Hλ or Vλ′ , the composition Cm Cm−1 · · · C2 C1 is a chain of operators. The concept of tolerance is extended from operators to chains of operators. For an integer m ≥ 1, a real µ ≥ 0 and a region P, we say that a chain of m operators applied to P (that is, Cm · · · C1 (P)) has tolerance up to µ if, for each i = 1, . . . m, the ith operator Ci has tolerance lesser or equal than µ. Note that Ci is applied to the region Ci−1 · · · C1 (P), and hence the tolerance of Ci depends on the diameter (horizontal or vertical according to the type of Ci ) of Ci−1 · · · C1 (P), and not on the diameter of P. We are interested in the reduction in )diameter as function of the number of cuts m. Lemma 6 in Appendix A gives us ( 1+µ t0 , as function of the minimum number t0 of operators of each type in the chain the factor of reduction in diameter, 2 Cm Cm−1 · · · C1 (P). To relate the length m of a chain with the minimum t0 , we need to discuss the elongation, or aspect ratio of P. The fact that the aspect ratio is connected with the relationship between t0 and m can be seen considering an initial region with a highly elongated shape across the width (Fig. 8). The chain of cuts along the minor axis starts with a first subsequence of vertical cuts. The third cut produces a subregion roughly as high as wide. After these vertical cuts, it comes a second subsequence of cuts, also along minor axis, that happens to be alternating between horizontal and vertical types. That is, the chain C4 C3 C2 C1 (P) has m = 4 and t0 = 1. For a region with higher elongation, the number m of cuts up to the first type change is greater, while t0 = 1.

154

J.L. García-Zapata, J.C. Díaz Martín and Á. Cortés Fácila / Journal of Computational and Applied Mathematics 352 (2019) 146–164

Fig. 8. Six cuts along minor axis made in a elongated region P. They start with C1 , C2 and C3 , each with some shifting from the central dashed line, produced by ISC. After these, C4 , C5 and C6 alternate horizontal and vertical type.

Fig. 9. A pair of horizontal and vertical cuts along middle line does not alter aspect. However, after pairs of horizontal and vertical cuts with biased shifts, the resulting aspect grows indefinitely.

In general, a chain of cuts along minor axes does not have a first subsequence of cuts of the same type and, after that, a second and last subsequence of cuts of alternating type. By virtue of the shifting of ISC, the subregion that arises can have high elongation, like in Fig. 8 after C6 , requiring in turn a subsequence of cuts along minor axis all of the same type, vertical in this case. To precise the relationship between the length m of a chain with the minimum t0 of operators of some type, we introduce the aspect ratio. dm (P) dm (P) The horizontal aspect ratio of a region P is aspH (P) = dmH (P) , and the vertical aspect ratio is its reciprocal aspV (P) = dmV (P) . V

H

The aspect ratio (for short, aspect) is asp(P) = max(aspH (P), aspV (P)). Fig. 9(a) shows the subregions that arise after the application of cuts along the middle line, alternatively T and R. The initial region of the figure, P0 , is the unit square (of aspect 1), and its subregion P1 = T (P0 ) has aspH (P1 ) = 2 and aspV (P1 ) = 1/2 . The subregion P2 = R(P1 ) of P1 has aspect 1. Calling Pi = T (Pi−1 ) if i is odd and Pi = R(Pi−1 ) if i is even, the aspect of the regions Pi alternates between 1 and 2. In general, after a horizontal and a vertical cut (by the middle line), the resulting region has equal aspect than the initial. With shifted cuts a phenomenon appears: the aspect can increase with the number of cuts. Furthermore, the increasing can be an exponential function of this number. This situation arises if the shifts of these cuts are repeatedly biased towards a determinate direction. Fig. 9(b) shows the subregions obtained after alternative application of Tλ and Rλ′ , being λ and λ′ corresponding with tolerance µ = 1/2 . The initial region P0′ is the unit square, like in the previous example, and Pi′ = Tλ (Pi′−1 ) if i is odd and Pi′ = Rλ′ (Pi′−1 ) if i is even. The aspects of P1′ , P2′ , P3′ , . . . are respectively 4, 3, 12, 9, 36, 27, . . . It can be seen ′ ′ i−1 that asp(P2i ) = 3i and asp(P2i . +1 ) = 4 · 3 Using Lemma 10 (in Appendix B) we prove the following result for the decreasing in diameter up to the accuracy A: Proposition 2. If P is a convex region and Cm Cm−1 · · · C2 C1 (P) is a chain of operators of tolerance up to µ < 1 such that Ci is a dm (P)

cut along minor axis of Ci−1 · · · C2 C1 (P), and m ≥

R log2 A L 1−log2 (1+µ) 1

+ L0 , with L0 = 2 +

log2 (asp(P)) 1−log2 (1+µ)

dmR (Cm Cm−1 · · · C2 C1 (P)) ≤ A. dm (P)

Proof. The hypothesis log2

dmR (P) A



R log2 A L 1−log2 (1+µ) 1

m − L0 L1

+ L0 ≤ m is equivalent to the following inequality:

(1 − log2 (1 + µ))

1−log (1−µ)

and L1 = 1 + 2 1−log2 (1+µ) , then 2

J.L. García-Zapata, J.C. Díaz Martín and Á. Cortés Fácila / Journal of Computational and Applied Mathematics 352 (2019) 146–164

Calling L2 = log2

m−L0 L1

this implies:

dmR (P)

≤ L2 (1 − log2 (1 + µ))

A

L2 log2 (1 + µ) + log2

( log2

dmR (P)

≤ L2 )A dmR (P) ≤ L2 (1 + µ)L2 A

(1 + µ)

(

155

L2 dmR (P)

1+µ

A

≤ 2L2

)L2 dmR (P) ≤ A

2

The inequality of Lemma 10 (Appendix B) in these terms is

( dmR (Cm · · · C1 (P)) ≤

1+µ

) m−L L0

2

1

( dmR (P) =

1+µ

)L2

2

dmR (P),

Chaining with the above, we have dmR (Cm · · · C1 (P)) ≤ A. □ 4.3. The shifting loop meets the requirements We will prove that the division loop using ISC along minor axis in RDP algorithm (Fig. 6) meets the requirements: decreasing diameter and IPSR without error. Proposition 1 bounds the maximum shifting |himax | as function of the parameter θ of IPSR. We will show in the following Proposition 3 that giving to θ the value ZΓ =

min(dmH (P), dmV (P))



4(4 + 2 2)(n + 2)n

.

(1)

we achieve the two goals: first, each final cut of ISC produces two non-empty regions where IPSR returns successfully: second, this cut has a tolerance lesser than 1/2 . Therefore, with this value of the parameter θ the chains of cuts produced by RDP (Fig. 6) have tolerance µ < 1/2 . Proposition 3. Let P be a region containing nP roots of a polynomial of degree n such that none of these roots is at a distance less √ or equal to E(ZΓ ) = (4 + 2 2)nZΓ from its border Γ , with ZΓ given by (1). (a) The shifting loop ends giving us two regions, P0 and P1 , such that IPSR(P0 , ZΓ ) and IPSR(P1 , ZΓ ) return successfully, with a n +1 cut of tolerance lesser or equal than µ = 2(nP +2) < 21 . (b) The subregions that arise from P after applying m recursive calls to RDP have rectangular diameter lesser or equal than a given tolerance A, if m verifies the hypothesis of Proposition 2 for the above value of µ. Proof. For (a), we discuss the case of a shifting loop of horizontal ISC, being similar the vertical case. Since n ≥ nP , we have dmV (P) √ < ΘH , where ΘH is given in Proposition 1. Therefore, by Proposition 1 if we consider θ = ZΓ , ISC(P , i) ZΓ = 4(4+2 2)(n+2)n gives us two subregions that do not return with error in IPSR, with separation from the central line of

√ |hi | ≤ (nP + 1)E(ZΓ ) = (4 + 2 2)(nP + 1)n =

dmV (P) nP + 1 2

2(n + 2)

dmV (P)



4(4 + 2 2)(n + 2)n

.

By definition, any horizontal cut with a separation of λ has a tolerance of given by ISC(P , i) is

2 dmV (P)

|hi | ≤

dmV (P) nP +1 2 dmV (P) 2 2(n+2)

=

nP +1 2(n+2)

<

1 . 2

2λ . dmV (P)

Hence the tolerance

2|hi | dmV (P)

of the cut

For (b), note that the subregions produced by m cuts, horizontal or vertical, are the result of a chain of operators n +1 Cm Cm−1 · · · C2 C1 (P), verifying the hypothesis of Proposition 2 for the value µ = 2(nP +2) . By this proposition we conclude that the rectangular diameters of the subregions are lesser or equal than A. □ According to the previous result, the bound for the number m of recursive calls required to reach a region with diameter lesser or equal than a given tolerance A starting from a region P is greater or equal than a quantity KL0 + L1 , where L0 depends dm (P) on the logarithm of the aspect of P and K on the logarithm of the ratio AR . These logarithmic dependences were expected, because the term L0 counts the number of cuts required to decrease, by division, the aspect of P, if it is very elongated. The

156

J.L. García-Zapata, J.C. Díaz Martín and Á. Cortés Fácila / Journal of Computational and Applied Mathematics 352 (2019) 146–164

resultant region after L0 cuts, of low aspect, is then reduced by division down to the required accuracy A, with a number of cuts accounted by K . We can give another lower bound, with simpler expression (that is, without terms L0 or L1 ) for the number of cuts required: Corollary 1. The regions obtained after m recursive calls of RDP have rectangular diameter lesser or equal than A, if m≥

6 − log2 3 (2 − log2 3)2

log2

dmR (P)

+

A

Proof. The bound of above proposition on the tolerance µ =

nP +1 . 2(n+2)

dmR (P) A

1 − log2 (1 + 1/2)

=

log2

dmR (P) A

2 − log2 3

with L′ = 1 + log2

L′ + 2 +

L′ + 2 +

1−log (1−1/2) 2 1−log2 (1+1/2) 2 dmR (P) A

1 − log2 (1 + µ)

2 − log2 3 log2

dmR (P) A

L 1−log2 (1+µ) 1

+ 2. log2 (1−µ) + L0 , with L1 = 1 + 2 11− and L0 = 2 + −log (1+µ)

The terms on this expression,

expression is also increasing. As µ = log2

log2 (asp(P))

=

nP +1 2(n+2)

<

dmR (P) log2 A

1−log2 (1+µ)

2

log2 (asp(P)) , 1−log2 (1+µ)

depends

, L1 and L0 , are increasing functions of µ, hence the

the bound is lesser than the value taken by the expression for 1/2 , that is

1 , 2

log2 (asp(P)) 1 − log2 (1 + 1/2)

log2 (asp(P)) 2 − log2 3

6−log2 3 . 2−log2 3

log2

Therefore

dmR (P) A

6 − log2 3 log2 (asp(P)) +2+ 2 − log2 3 2 − log2 3 2 − log2 3 6 − log2 3 dmR (P) log2 (asp(P)) = log2 + +2 □ (2 − log2 3)2 A 2 − log2 3

L1 + L0 ≤

5. Correctness of the recursive procedure In this section we show that the output of the RDP procedure is an array of approximations to the roots of f inside P, with their respective multiplicities. The following theorem ensures that the recursive procedure RDP (see Fig. 6) ends and that its outputs√ verify the claimed properties if the initial region PI does not have any root of f a distance less or equal than E(ZΓ ) = (4 + 2 2)nZΓ from its border ΓI , for ZΓ given by (1). Theorem 2. For a polynomial f of degree n, if PI is a plane region with dmR (PI ) > A without roots of f at a distance less or equal √ than E(ZΓ ) = (4 + 2 2)nZΓ from its border ΓI , for ZΓ given by (1), the procedure RDP applied to PI with accuracy A verifies the following properties. (a) It ends, and the recursive calls reach a maximum recursion depth of

⌈ lmax =

6 − log2 3 (2 − log2 3)2

log2

dmR (PI ) A

+

log2 (asp(PI )) 2 − log2 3



+2 .

(b) At the end, the plane regions of array Π = (P1 , P2 , . . . , Pk ) are approximations of all the roots of f in PI , each containing the number of roots given by N = (n1 , n2 , . . . , nk ). Proof. For (a), suppose that we reach a call of recursion depth m, that is RDP(Cm Cm−1 · · · C2 C1 (PI ), A). The parameter region arises after performing m cuts Ci . Each cut Ci is made by the last iteration of the shifting loop in a RDP call of depth i. All these 6−log 3 dm (P ) cuts are made by minor axis, with tolerance up to 12 by Proposition 3(a). Hence by Corollary 1, if m ≥ (2−log 23)2 log2 RA I + 2

log2 (asp(PI )) 2−log2 3

+ 2 then the rectangular diameter of the subregion Cm · · · C1 (PI ) is lesser or equal than A, returning therefore non ⌈ ⌉ 6−log 3 log (asp(P )) dm (P ) recursively by exit 2. Then, we can reach recursion depth m with m = (2−log 23)2 log2 RA I + 22−log 3I + 2 , but not 2

2

higher. For the finiteness of each RDP instance, in general a recursive procedure that makes a bounded number of recursive calls in each call, and that has a bounded recursion depth, ends in a finite number of steps. And this is the case of RDP, because as there can be only nP cut errors, the loop of division is done at most nP + 1 times. There are two recursive calls inside the loop, hence each call RDP(P , A) makes⌈2(nP + 1) or less recursive calls. Besides ⌉ we have shown that the maximum recursion depth reached by RDP in PI is lmax = steps.

6−log2 3 (2−log2 3)2

log2

dmR (PI ) A

+

log2 (asp(PI )) 2−log2 3

+ 2 . Hence RDP(PI , A) ends in a finite number of

J.L. García-Zapata, J.C. Díaz Martín and Á. Cortés Fácila / Journal of Computational and Applied Mathematics 352 (2019) 146–164

157

Fig. 10. Call tree of RDP, starting from RDP(PI , A). Each node is labeled with the region name and its A-height.

Now we prove part (b). We say that a plane region P requires an A-height v if the procedure RDP(P , A) reaches up a recursion depth v . Fig. 10 shows an example of the call tree of RDP. That is RDP(PI , A) recursively calls RDP(P0 , A) and RDP(P1 , A), that in turn call RDP(P00 , A) and RDP(P01 , A), and RDP(P10 , A) and RDP(P10 , A), respectively, and so on. The region of the call at depth 0, PI , has A-height 3, and P0 and P1 , though they both have depth 1, show A-height 1 and ⌈ 2, respectively. Every plane region P without roots at E(ZΓ ) or less from border Γ has a finite A-height, lesser or equal than

dmR (P) A

+

log2 (asp(P)) 2−log2 3



6−log2 3 (2−log2 3)2

log2

+ 2 as we saw in (a), because the reasoning made for PI is also valid for a general region P. We will prove

part (b) by induction on the A-height v . That is, we first prove (b) for regions P of A-height 0 (the base case), and then for regions of higher A-height, assuming it true for lower heights. So, the hypothesis√ of induction, that we denote HI(v ), is ‘‘At the end of the procedure RDP(P , A) (for any plane region P without roots at (4 + 2 2)nZΓ or less from border Γ , of A-height v ) the plane regions of Π have rectangular diameter lesser than A, contain the number of roots given by N, and these are all the roots of f contained in P’’. We will prove first the base case HI(0), and then the general case: if HI(v ′ ) for each v ′ < v , then HI(v ). Note that this induction is complete, meaning that the v th case depends on all the previous cases, not only the immediately preceding v − 1. It is a structural induction [16,17], a method to prove properties of recursive procedures. The base case (v = 0) applies to the regions that do not require recursive calls, those without roots (that return by exit 1) or with one or more roots and with rectangular diameter lesser than A, (that return by exit 2). On returning by exit 1 the assertion HI(0) is trivially verified because P does not have roots and Π (and N) remains empty. On returning by exit 2 the assertion is also verified because Π contains only the region P , and IPSR computes without error the winding number (by hypothesis HI(0), there are no roots near the border of P). Therefore N contains the number of roots inside P . For the general case (v ≥ 1), we consider the two recursive calls RDP(Pi , A), being P0 and P1 disjoint regions that cover P, produced by some iteration of ISC, that does not return with error in IPSR. P has A-height v , that is, RDP(P , A) reaches up a recursion depth of v . As RDP(P0 , A) and RDP(P1 , A) are one recursion level under RDP(P , A), the A-heights of P0 and P1 , say v0 and v1 , are not necessarily the same, but both are equal or lesser than v − 1. Let we call Πi , for i = 0 or 1, the array of regions added to Π by RDP(Pi , A). Similarly Ni is the array of numbers added to N. By induction hypothesis HI(vi ) (that we suppose true because vi < v ) at the end of the recursive call RDP(Pi , A) the regions pertaining to Πi are all of rectangular diameter lesser than A, contain the number of roots given by Ni , and these are all the roots in Pi . Then the concatenation of Π1 and Π2 (that is Π ) and also that of N1 and N2 (that is N) verify that its regions have diameters and number of roots as was respectively specified. That these are all the roots inside P comes from the fact that the Pi are disjoints and cover P. This is what we want to show, HI(v ). □ 6. Conclusions We have developed a method for root finding of polynomials. It allows to restrict the search to an area of interest in the complex plane, unlike the more common iterative algorithms. It is a combination of an inclusion test, IPSR (to see if a plane region contains some roots) based on the winding number, with a recursive method of partition in subregions, RDP, like other known algorithms [6,18–21]. The contribution of our method is to use a robust inclusion test, that never fails. This provides two properties. The first is that RDP always ends with an approximation to the solution within its accuracy specification, which has made possible to develop a practical implementation. The second is that it allows a formal proof of correctness, provided in this article, and also of computational cost. RDP is the first geometric method with these properties. The computational cost has been calculated for the worst case using techniques of analysis of algorithms. For a given degree, as function of the required accuracy A, is of order ( standard ) O log A1 . Although our interest in this development is mainly practical, this result can be compared with other theoretical bounds, being better than [22] and [21]. This is only indicative and certainly more work need to be done to compare the performance of the proposed method with others existing in the literature. The order with respect to the degree, which is O(n2 nI ), does not look promising, being O(n log n) that of [21]. However it should be noted that the factor nI is the number of roots contained in the search region, which may be much less than the total number of roots.

158

J.L. García-Zapata, J.C. Díaz Martín and Á. Cortés Fácila / Journal of Computational and Applied Mathematics 352 (2019) 146–164

As a practical result of this research, we have produced a computer implementation, Contour, that efficiently solve illconditioned polynomials up to degree 500 with the standard floating-point resolution of current commodity computers. It can be tested in the web address http://gim.unex.es/contour. As future work, we are developing a parallel implementation, taking advantage of the task parallelism implicit in geometrical methods. Acknowledgments This work has been funded by the ‘‘IV Plan Regional de I+D+i’’, ‘‘Consejeria de Empleo, Empresa e Innovacion’’ of the Extremadura Local Government. Appendix A We prove here Lemma 6 about diameter decreasing, after some preliminary results. Lemma 1. If two plane regions P1 , P2 verify P1 ⊂ P2 , then dm(P1 ) ≤ dm(P2 ), dmH (P1 ) ≤ dmH (P2 ), and dmV (P1 ) ≤ dmV (P2 ). Proof. For the classical diameter is clear because the maximum, in the definition of dm(P2 ), is taken over a greater set than in dm(P1 ). For the horizontal diameter, note that the two vertical supporting lines of P1 are between the two vertical supporting lines of P2 , perhaps coinciding with some of them. Hence its distance is lesser or equal. Likewise for the vertical diameter. □ Lemma 2. For a plane region P, it is verified: max(dmH (P), dmV (P)) ≤ dm(P) ≤ dmR (P) Proof. For the first inequality, note that each supporting line of P has as least one point in common with P. Being x0 one of these points for the upper horizontal supporting line, that is x0 ∈ P ∩ lT , and y0 for the lower one, y0 ∈ P ∩ lB , using the definitions of dmV (P) and dm(P) we have dmV (P) = minx∈lT , y∈lB d(x, y) ≤ d(x0 , y0 ) ≤ maxx, y∈P d(x, y) = dm(P). Likewise dmH (P) ≤ dm(P), hence max(dmH (P), dmV (P)) ≤ dm(P). For the last inequality, we define the HV-envelope of P, EnvHV (P), as the bounded rectangle delimited by the horizontal and vertical supporting lines. As P ⊂ EnvHV (P), by Lemma 1 we have that dm(P) ≤ dm(EnvHV (P)). Furthermore, the diameter of the rectangle EnvHV (P), of base dmH (P) and height dmV (P), is the distance between opposed vertexes, then dm(EnvHV (P)) =



dmH (P)2 + dmV (P)2 = dmR (P). Chaining with the previous inequality, we conclude. □

Lemma 3. If λ ≥ 0, then B(P) ⊂ Bλ (P) and R(P) ⊂ Rλ (P). If λ ≤ 0, then B(P) ⊃ Bλ (P) and R(P) ⊃ Rλ (P). For the operators Tλ and Lλ , in the same hypothesis the inclusions are reversed. Proof. It is immediate considering the position of mH +λ (P) and mH (P) (or of mV +λ (P) and mV (P)).



Hλ denotes a generic horizontal operator, that is, either Tλ or Bλ . Similarly Vλ denotes either Lλ or Rλ . Note that for λ (or −λ) large enough in absolute value, the result of an operator can be the empty set, whose diameter is conventionally zero. Lemma 4. For a horizontal operator Hλ , dmH (Hλ (P)) ≤ dmH (P) and For a vertical operator Vλ , dmV (Vλ (P)) ≤ dmV (P) and

dmH (P) 2

dmV (P) 2

− |λ| ≤ dmV (Hλ (P)) ≤

− |λ| ≤ dmH (Vλ (P)) ≤

dmH (P) 2

dmV (P) 2

+ |λ|.

+ |λ|.

Proof. We discuss the horizontal operators, being similar the vertical ones. The first inequality is a consequence of Lemmas 1 dmV (P) and 3. For the two chained inequalities, note that if |λ| > , then Hλ (P) is the empty set and the inequalities are trivially 2 dmV (P) true. In the alternative case, |λ| ≤ , the defining line of Hλ (P), mH +λ (P), will be placed above mH (P), if either λ ≥ 0 2 and Hλ is Tλ , or λ ≤ 0 and Hλ is Bλ , and below in other cases. Anyway it is at a distance of |λ| from mH (P). This horizontal dm (P) line mH (P) is at a distance of 2V from lT if Hλ is Tλ (alternatively lB if Hλ is Tλ ). As mH +λ (P) and lT (or lB ) are the supporting lines of dmV (Hλ (P)), hence the inequalities are verified. □ Lemma 5. If Hλ (P) has a tolerance of µ, then 1+µ dmV (P) ≤ dmV (Hλ (P)) ≤ dmV (P). 2 2 Likewise, if Vλ (P) has a tolerance of µ, then 1−µ

1−µ 2

dmH (P) ≤ dmH (Vλ (P)) ≤

1+µ 2

dmH (P).

J.L. García-Zapata, J.C. Díaz Martín and Á. Cortés Fácila / Journal of Computational and Applied Mathematics 352 (2019) 146–164

Proof. By the definition of tolerance |λ| = µ operator: dmV (Hλ (P)) ≥

dmV (P)

− |λ| =

2

dmV (P) 2

dmV (P) 2

159

. We substitute this into the inequalities of Lemma 4 for a horizontal

−µ

dmV (P)



dmV (P)

2

= (1 − µ)

dmV (P)

= (1 + µ)

dmV (P)

2

.

and dmV (Hλ (P)) ≤

dmV (P)

+ |λ| =

2

dmV (P) 2

2

2

.

The case of a vertical operator is similar. □ If a horizontal or vertical cut has tolerance greater or equal than 1, the above lemma is trivially true, and insubstantial, because one of the subregions coming out of the cut has zero diameter, and the other subregion has the same diameter than P. Lemma 6. For m ≥ 1, let us suppose that the chain Cm Cm−1 · · · C1 has a number tH of operators of type horizontal and tV of type vertical, so that tH + tV = m. If Cm · · · C1 (P) has tolerance up to µ with µ < 1, then

( dmR (Cm Cm−1 · · · C1 (P)) ≤

1+µ

) t0 dmR (P)

2

being t0 = min(tH , tV ).

( 1+µ )tV

( 1+µ )tH

dmH (P) and dmV (Cm Cm−1 · · · C1 (P)) ≤ dmV (P), Proof. We will first show that dmH (Cm Cm−1 · · · C1 (P)) ≤ 2 2 and then the claim about rectangular diameter. With respect to the horizontal diameter dmH (Cm Cm−1 · · · C1 (P)), we consider in general Ci Ci−1 · · · C1 (P) for i = m, m − 1, . . . , 1. If Ci is horizontal, by the inequality dmH (Hλ (P)) ≤ dmH (P) of Lemma 4, dmH (Ci Ci−1 · · · C1 (P)) ≤ dmH (Ci−1 · · · C1 (P)). 1+µ 2

If Ci is vertical, by the inequality dmH (Vλ (P)) ≤ dmH (Ci Ci−1 · · · C1 (P)) ≤

1+µ 2

dmH (P) of Lemma 5,

dmH (Ci−1 · · · C1 (P)).

Starting from dmH (Cm Cm−1 · · · C1 (P)), we will apply m times either one of the above inequalities, the first one if Ci is horizontal and the second one if Ci is vertical, in a iterative way. So, by the first inequality, if Cm is horizontal, or else by the second, if Cm is vertical, we have

( dmH (Cm · · · C1 (P)) ≤

1+µ

)βm dmH (Cm−1 · · · C1 (P))

2

being βm = 0 in the first case and βm = 1 in the second. Then we apply for second time one of the inequalities, to the factor dmH (Cm−1 · · · C1 (P)) that appears in the second member above. We apply the first inequality if Cm−1 is horizontal, the second one in the other case, giving

(

1+µ

)βm

( dmH (Cm−1 · · · C1 (P)) ≤

2

1+µ

)βm (

2

1+µ

)βm−1 dmH (Cm−2 · · · C1 (P))

2

being βm−1 = 0 if Cm−1 is horizontal and βm−1 = 1 if it is vertical. We apply for third time either one of the inequalities to the second member above, giving

(

1+µ

)βm (

2

( ≤

1+µ

)βm−1 dmH (Cm−2 · · · C1 (P))

2 1+µ

)βm (

1+µ

2

)βm−1 (

1+µ

2

)βm−2 dmH (Cm−3 · · · C1 (P))

2

where βm−2 is 0 if Cm−2 is horizontal and 1 in contrary case. Doing this for i = m, m − 1, . . . , 1, we have a chain of inequalities

( dmH (Cm · · · C1 (P)) ≤

1+µ

)βm dmH (Cm−1 · · · C1 (P))

2

( ≤

1+µ

)βm (

1+µ

2

( ≤

)βm−1

2

1+µ 2

)βm (

1+µ 2

dmH (Cm−2 · · · C1 (P)) ≤ · · ·

)βm−1

( ···

1+µ 2

)β1 dmH (P)

160

J.L. García-Zapata, J.C. Díaz Martín and Á. Cortés Fácila / Journal of Computational and Applied Mathematics 352 (2019) 146–164

where βi = 0 if Ci is horizontal and βi = 1 if Ci is vertical. Hence the sum of the exponents

( 1+µ )tV

vertical operators in Cm , Cm−2 , . . . , C1 . Therefore dmH (Cm · · · C1 (P)) ≤ Likewise dmV (Cm · · · C1 (P)) ≤ Hence dmR (Cm · · · C1 (P)) =

√ √

( 1+µ )tH

βi is tV , the number of

dmV (P), using the analogous inequalities of Lemmas 4 and 5 for vertical diameter.

2

1+µ

)2tV

1+µ

( dmH

2

=

i=1

dmH (P).

dmH (Cm · · · C1 (P))2 + dmV (Cm · · · C1 (P))2

(



(

2

∑m

) t0



(

(P)2

1+µ

2

+

1+µ

)2tH dmV (P)2

2

)2(tV −t0 )

( dmH

2

(P)2

+

1+µ 2

)2(tH −t0 ) dmV (P)2

One of the exponents, 2(tV − t0 ) or 2(tH − t0 ), is equal suppose that it is the first (equivalently t0 = tV ), being the ( to 0. Let )2(tus H −t0 ) 1 + µ 1+µ ≤ 1. Therefore other case similar. Besides 2 ≤ 1, which implies 2

(

1+µ

) t0



(

1+µ

2

=

1+µ

) t0



1+µ

(

dmH (P)2 +

2

(

( dmH (P)2 +

2

(



)2(tV −t0 )

)t0 √

2

1+µ

1+µ

)2(tH −t0 ) dmV (P)2

2

)2(tH −t0 ) dmV (P)2

2

( dmH (P)2 + dmV (P)2 =

1+µ

)t0

2

dmR (P). □

Appendix B The aim of this appendix is to prove Lemma 10. Lemma 7. If P is a convex region, and Ck Ck−1 · · · C2 C1 (P) is a chain of at least two operators, of tolerance up to µ, such that Ci is log2 (asp(P)) , then the operators Ci cannot be all of the same type, horizontal a cut along minor axis of Ci−1 · · · C1 (P), and k > 1 + 1−log 2 (1+µ) or vertical. Proof. We will show that the operators cannot be all horizontal, being similar the proof of the other impossibility. If Ci is horizontal for i = 1, 2 . . . , k, we will reach a contradiction. As Ck is horizontal then the minor axis of Ck−1 · · · C1 (P) is horizontal, that is dmH (Ck−1 · · · C1 (P)) ≤ dmV (Ck−1 · · · C1 (P)), or 1≤

dmV (Ck−1 · · · C1 (P)) dmH (Ck−1 · · · C1 (P))

.

Besides, by the inequalities of Lemma 4 (in the denominator) and Lemma 5 (in the numerator) for horizontal operators Ci , beginning with i = k − 1 and going down to i = 1: dmV (Ck−1 · · · C1 (P)) dmH (Ck−1 · · · C1 (P))

( 1+µ )2

Hence 1 ≤

dmH (Ck−2 · · · C1 (P))

dmH (Ck−3 · · · C1 (P))

( 1+µ )k−1 2

( 1≤

(

dmV (Ck−2 · · · C1 (P))

dmV (Ck−3 · · · C1 (P))

2



1+µ 2



1+µ

dmV (P) . dmH (P)

Note that asp(P) =

)k−1

2 2

( 1+µ )k−1 ≤ ··· ≤

asp(P)

)k−1

≤ asp(P) ( ) 2 (k − 1) log2 ≤ log2 (asp(P)) 1+µ 1+µ

(k − 1) ≤

log2 (asp(P)) log2

2 1+µ

.

2

dmV (P)

dmH (P) dmV (P) , dmH (P)

because C1 is horizontal. Then solving for k:

J.L. García-Zapata, J.C. Díaz Martín and Á. Cortés Fácila / Journal of Computational and Applied Mathematics 352 (2019) 146–164

161

Fig. B.11. A convex region P (thick line), the convex hull EC of its cuts with the supporting lines (thin), and the triangle DT that joins any point of IT with the opposite line (dotted).

as

log2 (asp(P))

k>

log2 (asp(P)) , 2 1−log2 (1+µ) log2 1+µ log2 (asp(P)) 1 . □ 1−log2 (1+µ)

=

this implies (k − 1) ≤

log2 (asp(P)) . 1−log2 (1+µ)

That k ≤ 1 +

log2 (asp(P)) 1−log2 (1+µ)

is in contradiction with the hypothesis

+

Lemma 8. Let us suppose that P is convex and µ < 1. If Hλ (P) has a tolerance of µ, then 1−µ 2

dmH (P) ≤ dmH (Hλ (P)).

If Vλ (P) has a tolerance of µ, then 1−µ 2

dmV (P) ≤ dmV (Vλ (P)).

Proof. The intersection of P with the upper supporting line lT is, by convexity, a segment IT perhaps reduced to a point. Likewise IB , IL and IR are the intersections of P with the lines lB , lL and lR , respectively (Fig. B.11). Let us call EC the convex hull [23] of the segments IT , IR , IB , IL . As P is convex, it contains this convex hull EC . We first prove the inequality for horizontal operators in the case of Tλ . We consider the triangle DT that joins any point of IT with the base of the rectangle of supporting lines that encloses P, shown in Fig. B.11. Consider the two points that are at maximum height in IL and IR respectively, and let us take the lowest one. For example in Fig. B.11 this is the unique point that pertains to IR . An horizontal line M over this point divides the height of P in two intervals. We will consider two cases: first, the line mH +λ (P) that produces Tλ (P) is above or coincides with M; second, mH +λ (P) is below M. In the first case, Tλ (P) contains the part of the convex hull EC above mH +λ (P), which in turn contains the part of the triangle DT above this line. That is Tλ (P) ⊃ EC ∩ Tλ (P) ⊃ DT ∩ Tλ (P). This chain of inclusions, by Lemma 1, implies 1−µ dmH (Tλ (P)) ≥ dmH (EC ∩ Tλ (P)) ≥ dmH (DT ∩ Tλ (P)). We will show that dmH (DT ∩ Tλ (P)) ≥ 2 dmH (P). The segment dmV (P) A = mH +λ (P) ∩ DT gives us the horizontal diameter of DT ∩ Tλ (P). A is at height + λ in the triangle DT , which has base 2 ) ( dmH (P) and height dmV (P). Therefore by a rule of three it has a length of

dmV (P) 2

−λ

dmH (P) . dmV (P)

Moreover, as

2|λ| dmV (P)

= µ, then

162

J.L. García-Zapata, J.C. Díaz Martín and Á. Cortés Fácila / Journal of Computational and Applied Mathematics 352 (2019) 146–164

2λ dmV (P)

≤ µ and then λ ≤ µ dm2V (P) , hence −λ ≥ − µ dm2V (P) . So ) ( ) ( dmH (P) 1 µ dmV (P) −λ ≥ − dmH (P) length(A) = 2

dmV (P)

2

2

1−µ 2

1−µ

and dmH (DT ∩ Tλ (P)) = length(A) ≥ dmH (P). By the above chain of inequalities, dmH (Tλ (P)) ≥ 2 dmH (P). 1−µ 1+µ 1+µ In the second case, as µ is positive, −µ ≤ µ hence 2 < 2 , and as µ < 1, 2 < 1. Then chaining these inequalities 1−µ 1−µ < 1, therefore 2 dmH (P) ≤ dmH (P). This was true also in the first case, but in this second case besides Tλ (P) contains 2 at least one point of IL and one of IR (the points at maximum height of IL and of IR are placed in Tλ (P)). Hence the width of 1−µ Tλ (P) coincides with the width of P, that is dmH (P) = dmH (Tλ (P)). Hence 2 dmH (P) ≤ dmH (Tλ (P)) and we conclude. For the other operators Lλ , Rλ and Bλ , the above proof can be rewritten with the proper changes. □ Lemma 9. If P is a convex region, C2 C1 (P) is a chain of operators of tolerance up to µ, one of horizontal type and the other vertical, such that C2 is a cut along minor axis of C1 (P) and C1 along minor axis of P, then asp(C2 C1 (P)) ≤

2(1 + µ) (1 − µ)2

Proof. Let us suppose that C2 is a vertical operator Vλ′ and C1 horizontal, Hλ , being the contrary case similar. dm (V ′ H (P)) dm (V ′ H (P)) By definition asp(Vλ′ Hλ (P)) is the maximum of aspV (Vλ′ Hλ (P)) = dmV (Vλ′ Hλ (P)) and aspH (Vλ′ Hλ (P)) = dmH (Vλ′ Hλ (P)) . By the H λ λ V λ λ inequalities of Lemmas 4 and 5, we have dmV (Vλ′ Hλ (P)) dmH (Vλ′ Hλ (P))



dmV (Hλ (P)) 1−µ 2

=

dmH (Hλ (P))

2

dmV (Hλ (P))

1 − µ dmH (Hλ (P))



2 1−µ

The last inequality comes from the fact that the cut along minor axis of Hλ (P) is vertical, then dmV (Hλ (P)) ≤ dmH (Hλ (P)), dm (H (P)) that is dmV (Hλ (P)) ≤ 1. H λ By the inequalities of Lemma 5 (numerator) and Lemma 8 (denominator), dmH (Vλ′ Hλ (P)) dmV (Vλ′ Hλ (P))

≤ =

1+µ 2 1−µ 2

1+µ

dmH (Hλ (P))

dmH (P)

≤ ( 2 )2 1−µ

dmV (Hλ (P))

2(1 + µ) dmH (P) (1 − µ)2 dmV (P)

2



dmV (P)

2(1 + µ) (1 − µ)2 dm (P)

The last inequality comes from that the cut along minor axis of P is horizontal, then dmH (P) ≤ dmV (P), that is dmH (P) ≤ 1. V As asp(Vλ′ Hλ (P)) = max (aspV (Vλ′ Hλ (P)), aspH (Vλ′ Hλ (P))), and as(the maximum ) of two values is bounded by the maximum of the bounds of each value, we have asp(Vλ′ Hλ (P)) ≤ max

and hence multiplying by

2 , 2 1−µ 1−µ



2(1+µ) . (1−µ)2

2 1−µ

2(1+µ) , (1 , and we conclude, because 1 ≤ −µ)2

1+µ 1−µ



Finally, we apply the above lemmas to get a reduction in diameter: Lemma 10. If P is a convex region and Cm Cm−1 · · · C2 C1 (P) is a chain of operators of tolerance up to µ, with µ < 1, and such that Ci is a cut along minor axis of Ci−1 · · · C1 (P), then

( dmR (Cm · · · C1 (P)) ≤ where L0 = 2 +

log2 (asp(P)) 1−log2 (1+µ)

1+µ 2

) m−L L0 1

dmR (P) 1−log (1−µ)

and L1 = 1 + 2 1−log2 (1+µ) . 2

( 1+µ ) m−L0

m−L

L1 Proof. In case of a short chain (that means m < L0 , equivalently m − L0 < 0, hence L 0 < 0) the factor 2 is greater 1 that 1, and the claim of the Lemma is trivially true. In other case, we will divide the chain Cm Cm−1 · · · C1 in t + 1 subchains, determined by subindexes k1 , k2 , . . . , kt , (for a particular value t). That is, the first subchain is Ck1 · · · C1 , the second Ck2 Ck2 −1 · · · Ck1 +2 Ck1 +1 , etcetera. So we will have Cki+1 Cki+1 −1 · · · Cki +2 Cki +1 for i = 1, . . . , t − 1, and a final subchain Cm Cm−1 · · · Ckt +1 . We will define the ki in such a way that every subchain (except the last one) contains at least one operator of each type. We next bound the number of operators in every subchain, that is the lengths k1 of the first one, ki+1 − ki (for i = 1, . . . , t − 1) and m − kt for the others. We define k1 as the index of the first type change in Cm Cm−1 · · · C1 , that is, Ck1 −1 Ck1 −2 · · · C2 C1 are of the same log2 (asp(P)) , then taking integer part type, horizontal or vertical, different from that of Ck1 . As real number L0 = 2 + 1−log (1+µ)

⌊L0 ⌋ > 1 +

log2 (asp(P)) . We apply Lemma 7 to the chain C⌊L0 ⌋ C⌊L0 ⌋−1 1−log2 (1+µ)

2

· · · C1 (P), that hence has at least two operators of different type. Hence k1 ≤ ⌊L0 ⌋. This implies that the length of the first subchain, k1 , meets k1 ≤ L0 , fact that will be used below.

J.L. García-Zapata, J.C. Díaz Martín and Á. Cortés Fácila / Journal of Computational and Applied Mathematics 352 (2019) 146–164

163

Let we call P0 = P, and Pi = Ci Ci−1 · · · C1 (P) for i = 1, . . . , m. As Ck1 and Ck1 −1 are of different type, by Lemma 9 2(1+µ) asp(Pk1 ) = asp(Ck1 Ck1 −1 (Pk1 −2 )) ≤ (1−µ)2 . Let k2 the index of the next type change after k1 , that is, the operators of Ck2 −1 Ck2 −2 · · · Ck1 +2 Ck1 +1 are of the same type, different of that of Ck2 . For example two chains of the types HVVVHVVV or VHHHHVVV have both k1 = 4 and k2 = 8. Note that in this way the subchain Ck2 −1 Ck2 −1 · · · Ck1 +2 Ck1 +1 has as least two operators. Applying the reciprocal of Lemma 7 to the chain Ck2 −1 · · · Ck1 +1 (Pk1 ), as the operators are of the same type, we have that the number of operators that it contains, that is k2 − k1 − 1, is lesser or equal than 1 + 2(1+µ) log2 (1−µ)2 . 1−log2 (1+µ)

then k2 − k1 − 1 ≤ 1 +

=

1

1−log2 (1+µ)

. As asp(Pk1 ) ≤

2(1+µ) , (1−µ)2

That implies

2(1+µ) (1−µ)2

1 + log2 (1 + µ) − 2 log2 (1 − µ) =2+ 1 − log2 (1 + µ) 1 − log2 (1 + µ) 2 − 2 log2 (1 + µ) + 1 + log2 (1 + µ) − 2 log2 (1 − µ)

k2 − k1 ≤ 2 +

=

log2

log2 (asp(Pk ))

1 − log2 (1 + µ) 3 − log2 (1 + µ) − 2 log2 (1 − µ) 1 − log2 (1 + µ)

=1+2

1 − log2 (1 − µ) 1 − log2 (1 + µ)

= L1 .

Hence the length of the second subchain verifies k2 − k1 ≤ L1 . The next type change index k3 after k2 also verifies k3 − k2 ≤ L1 , by a similar reasoning (the reciprocal of Lemma 7 applied to the chain Ck3 −1 · · · Ck2 +1 (Pk2 )). Hence the third subchain also has length lesser or equal than L1 . We repeat the reasoning for the following type change indexes, concluding that, if there are t type changes, of indexes k1 , k2 , . . . kt , we have k1 ≤ L0 and ki+1 − ki ≤ L1 for i = 1, . . . , t − 1. We also have that the subchain after the last type change, Cm Cm−1 · · · Ckt +1 , is composed of operators of the same type, hence by the reciprocal of log2 (asp(Pk ))

2(1+µ)

t Lemma 7 the number of operators that it contains is m − kt ≤ 1 + 1−log (1+µ , that operating with asp(Pkt ) ≤ (1−µ)2 implies ) 2 m − kt ≤ L1 − 1 and indeed m − kt ≤ L1 . As the t subchains that are not the final one are each composed of two or more operators, all of the same type except the last one, then there are at least t operators of each type, one in each such subchain. If tH is the number of operators of horizontal type in Cm Cm−1 · · · C1 , and tV of vertical type, we have the upper bound t ≤ min(tH , tV ). We will bound t from below accounting the number of operators (the length) of the first subchain, above upper-bounded by L0 , and of the other subchains, upper-bounded by L1 . The length of the whole chain, m, is the sum of the lengths of the subchains, then we have

m−L

m−L

m ≤ L0 + tL1 , that is, L 0 ≤ t. Hence L 0 ≤ min(tH , tV ), and, as 1 1 in the inequality of Lemma 6 we have that

( dmR (Cm · · · C1 (P)) ≤

1+µ 2

)min(tH ,tV )

( dmR (P) ≤

1+µ 2

1+µ 2

≤ 1, then

( 1+µ )min(tH ,tV ) 2



( 1+µ ) m−L L0 2

1

. By substitution

) m−L L0 1

dmR (P)

as we wanted to show. □ References [1] Anthony Ralston, Philip Rabinowitz, A first course in numerical analysis, Second, McGraw-Hill Book Co., New York, 1978, p. xix+556, International Series in Pure and Applied Mathematics. [2] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, D. Sorensen, LAPACK Users’ Guide, third ed., Society for Industrial and Applied Mathematics, Philadelphia, PA, 1999. [3] B.T. Smith, J.M. Boyle, J.J. Dongarra, B.S. Garbow, Y. Ikebe, V.C. Klema, C.B. Moler, Matrix eigensystem routines—EISPACK guide, Second, Springer-Verlag, Berlin, 1976, p. ix+551, Lecture Notes in Computer Science, Vol. 6. [4] Steven Fortune, An iterated eigenvalue algorithm for approximating roots of univariate polynomials, J. Symbolic Comput. 33 (5) (2002) 627–646, Computer algebra (London, ON, 2001). [5] J.H. Wilkinson, The Algebraic Eigenvalue Problem, Clarendon Press, Oxford, 1965, p. xviii+662. [6] Victor Y. Pan, Solving a polynomial equation: some history and recent progress, SIAM Rev. 39 (2) (1997) 187–220. [7] Dario A Bini, Leonardo Robol, Solving secular and polynomial equations: A multiprecision algorithm, J. Comput. Appl. Math. 272 (2014) 276–292. [8] Victor Y Pan, Elias Tsigaridas, Accelerated approximation of the complex roots and factors of a univariate polynomial, Theoret. Comput. Sci. 681 (2017) 138–145. [9] Chee K Yap, Michael Sagraloff, A simple but exact and efficient algorithm for complex root isolation, in: Proceedings of the 36th International Symposium on Symbolic and Algebraic Computation, ACM, 2011, pp. 353–360. [10] J.L. García Zapata, J.C. Díaz Martín, Finding the number of roots of a polynomial in a plane region using the winding number, Comput. Math. Appl. 67 (3) (2014) 555–568. [11] J.L. García Zapata, J.C. Díaz Martín, A geometric algorithm for winding number computation with complexity analysis, J. Complexity 28 (2012) 320–345. [12] Peter Henrici, Applied and Computational Complex Analysis. vol. 1, in: Wiley Classics Library, John Wiley & Sons Inc., New York, 1988, p. xviii+682, Power series—integration—conformal mapping—location of zeros, Reprint of the 1974 original, A Wiley-Interscience Publication. [13] Peter Henrici, Methods of search for solving polynomial equations, J. Assoc. Comput. Mach. 17 (1970) 273–283. [14] Aaron Herman, Hoon Hong, Elias Tsigaridas, Improving root separation bounds, J. Symbolic Comput. 84 (2018) 25–56. [15] A.N. Kolmogorov, S.V. Fom¯ın, Introductory real analysis, Dover Publications Inc., New York, 1975, p. xii+403, Translated from the second Russian edition and edited by Richard A. Silverman, Corrected reprinting. [16] R.M. Burstall, Proving properties of programs by structural induction, Comput. J. 12 (1) (1969) 41–48. [17] John E. Hopcroft, Rajeev Motwani, Jeffrey D. Ullman, Introduction to Automata Theory, Languages, and Computation, second ed., Addison Wesley, 2000.

164

J.L. García-Zapata, J.C. Díaz Martín and Á. Cortés Fácila / Journal of Computational and Applied Mathematics 352 (2019) 146–164

[18] Xingren Ying, I. Norman Katz, A reliable argument principle algorithm to find the number of zeros of an analytic function in a bounded domain, Numer. Math. 53 (1–2) (1988) 143–163. [19] James Renegar, On the worst-case arithmetic complexity of approximating zeros of polynomials, J. Complexity 3 (2) (1987) 90–113. [20] Arnold Schönhage, The fundamental theorem of algebra in terms of computational complexity, Technical report, Mathematisches Institut Universität Tübingen, 1982. [21] C. Andrew Neff, John H. Reif, An efficient algorithm for the complex roots problem, J. Complexity 12 (2) (1996) 81–115. [22] Victor Y. Pan, Univariate polynomials: Nearly optimal algorithms for numerical factorization and rootfinding, J. Symbolic Comput. 33 (2001) 2002. [23] Donald Ervin Knuth, Axioms and Hulls, springer-Verlag Berlin, 1992.