Counting Solutions of a Polynomial System Locally and Exactly
CCounting Solutions of a Polynomial SystemLocally and Exactly
Ruben Becker Michael Sagraloff
Abstract
We propose a symbolic-numeric algorithm to count the number of solutions of a poly-nomial system within a local region. More specifically, given a zero-dimensional system f = · · · = f n = 0 , with f i ∈ C [ x , . . . , x n ] , and a polydisc ∆ ⊂ C n , our method aimsto certify the existence of k solutions (counted with multiplicity) within the polydisc. Incase of success, it yields the correct result under guarantee. Otherwise, no informationis given. However, we show that our algorithm always succeeds if ∆ is sufficiently smalland well-isolating for a k -fold solution z of the system.Our analysis of the algorithm further yields a bound on the size of the polydisc forwhich our algorithm succeeds under guarantee. This bound depends on local parameterssuch as the size and multiplicity of z as well as the distances between z and all othersolutions. Efficiency of our method stems from the fact that we reduce the problem ofcounting the roots in ∆ of the original system to the problem of solving a truncatedsystem of degree k . In particular, if the multiplicity k of z is small compared to the totaldegrees of the polynomials f i , our method considerably improves upon known completeand certified methods.For the special case of a bivariate system, we report on an implementation of ouralgorithm, and show experimentally that our algorithm leads to a significant improvement,when integrated as inclusion predicate into an elimination method. In this paper, we propose a randomized but certified (i.e. Las-Vegas type) algorithm, denoted
PolySol , to count the number of solutions of a zero-dimensional polynomial system F within a given polydisc ∆ ⊂ C n . Let F : f ( x ) = . . . = f n ( x ) = 0 , with f i ∈ C [ x ] = C [ x , . . . , x n ] for all i = 1 , . . . , n, (1)be a zero-dimensional polynomial system. We further assume that each of the coefficients c i,α of the polynomials f i = (cid:88) α =( α ,...,α n ) c i,α · x α = (cid:88) α =( α ,...,α n ) c i,α · x α · · · x α n n can be approximated to any desired precision. That is, for any given non-negative integer(precision) ρ , we can ask for a dyadic approximation ˜ c i,α ∈ − ρ · ( Z + i · Z ) of c i,α with | ˜ c i,α − c i,α | < − ρ for the cost of reading the approximations. There are only finitely many solution in complex projective n -space. a r X i v : . [ c s . S C ] D ec iven a polydisc ∆ = ∆ r ( m ) = { z ∈ C n : (cid:107) z − m (cid:107) ∞ < r } of radius r centered at m ,we aim to compute the number of solutions of F = 0 in ∆ . Here, solutions are counted withmultiplicity. As input, our algorithm PolySol receives (arbitrary good approximations of)the coefficients of F , the polydisc ∆ , and an integer K ∈ { , , . . . , d F } , where d F := max i d i is defined as the maximum of the degrees d i of the polynomials f i . As output, it returnsan integer k ∈ N ∪ {− } . If k = − , nothing can be said, that is, the algorithm fails toprovide an answer to our request. Otherwise, k equals the number of solutions of F = 0 in ∆ . In this case, we say that the method succeeds. We further show that our method alwayssucceeds if (1) r is small enough, (2) K ≥ k , and (3) the smaller polydisc ∆ (cid:48) := ∆ r (cid:48) ( m ) , with r (cid:48) := r n ( K +1) n , contains a k -fold solution of F . We also derive a bound on the size of r thatguarantees success of our method if the other two requirements are fulfilled. The given boundis adaptive in the sense that it does not only depend on global parameters such as the degreeand the size of the coefficients of the polynomials f i , but also on solution-specific parameters,that is, the multiplicity and the size of z as well as the distances between z and the othersolutions of F . Here, we state our main result for the special case, where F is defined overthe integers. For a more general statement, see Theorem 88. Theorem 1.
Suppose that z is a k -fold solution of a polynomial system F as in (11) withpolynomials f i ∈ Z [ x ] of total degree d i and with integer coefficients c i,α of bit-size less than τ F . Then, for any K ≥ k , there exists an L ∗ ∈ N with L ∗ = ˜ O (cid:18) D F · max i =1 ,...,n d F + τ F d i + D F · log( z ) + log( ∂ ( z , F ) − ) + ( K + 1) n · log( σ ( z , F ) − ) (cid:19) such that, with probability at least / , the algorithm PolySol ( F , ∆ , K ) returns k forany disc ∆ = ∆ r ( m ) with r ≤ − L ∗ and (cid:107) m − z (cid:107) ∞ < r . Here, we use the definitions log( x ) := max(1 , log max(1 , (cid:107) x (cid:107) ∞ )) , and d F := max i d i , D F := (cid:89) ni =1 d i ,σ ( z i , F ) := min j (cid:54) = i (cid:107) z i − z j (cid:107) ∂ ( z i , F ) := (cid:89) j (cid:54) = i (cid:107) z i − z j (cid:107) µ ( z j , F ) , where z , . . . , z N denote the distinct solutions of F and µ ( z i , F ) the multiplicity of z i . Notice that our method never yields the exact multiplicity of a solution, even in the casewhere there is a well separated k -fold solution z in ∆ . Instead, we only obtain the sum ofthe multiplicities of all solutions contained in ∆ . However, in the considered computationalmodel, where only approximations of the coefficients of the input polynomials are known, itis simply not possible to achieve a stronger result. This is due to the fact that arbitrary smallperturbations of the input already destroy the multiplicity structure of non-simple roots.We see a series of applications of our method. For instance, our method can be used toverify correctness of the result provided by a numerical (non-certified) method such as homo-topy (e.g. [Ver99Ver99; BHS+13BHS+13]) or subdivision methods (e.g. [MP09MP09; BCG+08BCG+08]). Correspondingimplementations of such methods (e.g. Bertini, PHCpack, axel) are available and have provento be efficient and reliable in practice. Suppose that such a method returns an approxima-tion ζ of a k -fold solutions z such that (cid:107) ζ − z (cid:107) ∞ < − L , however, without any guarantee on he correctness of the result . Now, in order to show correctness, we may run the algorithm PolySol with input F , K = k , and ∆ = ∆ n ( k +1) n · − L ( ζ ) . According to the abovetheorem, the method returns k if the claimed result is actually correct and L is large enough.Hence, we eventually succeed if the numerical solver provides a sufficiently good approxima-tion of z together with the correct multiplicity. Again, we remark that the method does notprovide a proof that there is exactly one root of multiplicity k , but only a proof that there k roots counted with multiplicity in ∆ .For polynomial systems that are defined over the integers, there exist complete and certifiedmethods (e.g. [Rou99Rou99; Laz09Laz09; BS16BS16]) to compute isolating regions for all solutions togetherwith the corresponding multiplicities, however, their possible application is limited in practice.In particular, if the polynomials f i are of large degree, the running time for the necessarysymbolic computations (e.g. that of a Gröbner Basis or resultants) becomes prohibitive.Combining our method with a numerical solver may instead yield a certified result on theexistence of solutions in a certain region.In Section 55, we report on preliminary implementation of our method for the specialcase of a bivariate system. That is, we integrated an implementation of our method in Bi-solve [BEK+13BEK+13; KS15KS15], a highly efficient algorithm for isolating the solutions of a bivariatepolynomial systems with integer coefficients. There, it serves as an inclusion predicate toverify the existence of a k -fold solution of the system. Compared to the original approach in Bisolve , we observe a considerable improvement with respect to running time and precisiondemand.
Overview of the Algorithm.
There exists a simple method, also known as Pellet’s The-orem, to count the number of roots of a univariate polynomial f ∈ C [ x ] in a disc D r ( m ) = { x ∈ C : | x − m | ≤ r } of radius r centered at a point m ∈ C . The method works as follows:We first compute the Taylor-expansion f [ m ]( x ) := f ( m + x ) = (cid:88) i ≥ c i · x i = (cid:88) i ≥ f ( i ) ( m ) i ! · x i at m and then check whether | c k | · r k > (cid:80) i (cid:54) = k | c i | · r i for some k . Notice that the latterinequality implies that the part c k · x k of f [ m ] of degree k dominates the remaining parts onthe boundary of the disc D r (0) . If this is the case, then D r ( m ) contains exactly k roots of f ,which follows directly from Rouché’s Theorem applied to f [ m ] and its degree k -part c k · x k .In [BSS+15BSS+15], we give sufficient conditions on r and the locations of the roots with respect to m such that the above inequality is fulfilled. In particular, for m being a k -fold root of f ,we give a bound r in terms of the degree of f and the separation of m such that Pellet’sTheorem applies for any r < r ; see Lemma 99 for details.Our algorithm PolySol can be considered as an extension of Pellet’s Theorem topolynomial systems. Similar as in the one-dimensional case, we make crucial use of the factthat, for a sufficiently small neighborhood ∆ of a k -fold solution z of F , the system F [ z ] : f ( x + z ) = · · · = f n ( x + z ) = 0 obtained by shifting each of the polynomials f i by z is dominated by terms of degree k orless. Hence, in order to study the local behavior of F at z , it should suffice to consider thetruncation F [ z ] ≤ k of F [ z ] , where we only consider the part f i [ z ] ≤ k = (cid:80) α : | α |≤ k c (cid:48) i,α · x α of each f i [ z ] = f ( x + z ) = (cid:80) α c (cid:48) i,α · x α that is of degree k or less. In fact, in Corollary 33, we prove3hat, for any K ≥ k , the system F [ z ] ≤ K has a k -fold solution at the origin, and we give abound on its separation in terms of the separation of z as a solution of the original system F . In Theorem 77, we even show that if K ≥ k , and if (cid:107) m − z (cid:107) ∞ < − L for a sufficientlylarge L , then we can work with F [ m ] ≤ K instead of F [ z ] . Namely, in this case, F [ m ] ≤ K has k solutions of norm less than · − L , whereas all remaining solutions have considerably largernorm, that is, larger than some value that does not depend on L .We now provide an overview of our approach. For the sake of simplicity, we omit technicaldetails and only give the main ideas. Also, we do not treat any special cases, which consider-ably simplifies the approach when compared to the actual algorithm as given in Section 33. Wefirst define L := (cid:100) log r n ( K +1) n (cid:101) such that r n ( K +1) n ≤ − L ≤ r n ( K +1) n = r (cid:48) . Obviously, wecannot check in advance whether the above requirements on m and L are fulfilled, however,we can check whether F [ m ] ≤ K has a cluster of solutions near the origin. For this, we usea complete and certified algorithm to compute isolating regions of all solutions of F [ m ] ≤ K that are contained in the polydisc ∆ = ∆ r (0) . Notice that if K is small compared to thedegrees of the polynomials f i , then the cost for computing the solutions of F [ m ] ≤ K is muchlower than solving the original system directly. In particular, for K = 1 , the truncated system F [ m ] ≤ K becomes a linear system in n variables. Now, suppose that ∆ contains k (cid:48) solutionsof F [ m ] ≤ K ( k (cid:48) does not have to be equal to k ) that are well separated from the remainingsolutions, then we are left to show that F [ m ] contains the same number of solutions in ∆ .For this, we use a generalization of Rouché’s Theorem that applies to analytic functions in n -dimensional complex space; see Theorem 66. This approach requires to compute a lowerbound LB for (cid:107)F [ m ] ≤ K ( x ) (cid:107) ∞ := max i | f i ( x ) | on the boundary of ∆ as well as a correspond-ing upper bound UB on the error (cid:107)F [ m ]( x ) ≤ K − F [ m ]( x ) (cid:107) ∞ that occurs when passing from F [ m ] to the truncated system F [ m ] ≤ K . While the computation of UB is straightforward (see(1111) in Section 33), the computation of LB is more involved. Namely, we first compute thehidden-variable resultant R (cid:96) := Res( F [ m ] ≤ K , x (cid:96) ) ∈ Q [ x (cid:96) ] with respect to each of the variables x (cid:96) ; see Section 22 for details on the hidden variable approach. The roots of R (cid:96) are the projec-tions of the solutions of F [ m ] ≤ K on the x (cid:96) -axis, and R (cid:96) is contained in the ideal given by thepolynomials f i [ m ] ≤ K , that is, there exist g (cid:96), , . . . , g (cid:96),n ∈ Q [ x ] with R (cid:96) = g (cid:96), · f i [ m ] ≤ K + · · · + g (cid:96), · f i [ m ] ≤ K . (2)Using a recent result [DKS13DKS13] on the arithmetic Nullstellensatz, we derive upper bounds on theabsolute value of the coefficients of the polynomials g (cid:96), ; see Corollary 22 and (1313) in Section 33.In addition, we use our results on Pellet’s Theorem from [BSS+15BSS+15] to derive a lower boundfor | R (cid:96) | on the boundary of the disc D r (0) ⊂ C , which is the projection of the polydisc ∆ into one-dimensional space; see Lemma 99. Combining the latter two bounds then yields LB .Finally, we check whether LB > UB , in which case we conclude from Rouché’s Theorem that F [ m ] has the same number of solution in ∆ as the truncated system F [ m ] ≤ K . If UB < LB ,we return − .In the analysis of our algorithm, we show that if (cid:107) m − z (cid:107) ∞ < r n ( K +1) n for a sufficientlysmall r , then LB approximately scales like c · r k for some constant C , whereas UB scales like C (cid:48) · r − ( K +1) L for some constant C (cid:48) . Thus, in this case, our algorithm eventually succeeds if K ≥ k . As already mentioned, we omitted many details in the above description. In particular,for completeness, we needed to address certain special cases. In particular, this comprises thecase where F [ m ] ≤ k has distinct solutions whose projections on one of the coordinate axis are(almost) equal or solutions at infinity that yield roots of the hidden variable resultant. We4how how to handle such situations by means of a random rotation of the coordinate systemwithout harming the claimed complexity bounds. Implementation for the Bivariate Case.
For the special case of a polynomial system F : f ( x , x ) = f ( x , x ) = 0 in two variables, with f , f ∈ Z [ x , x ] , we implementedour algorithm in Sage . As an oracle for computing an arbitrary good approximation of asolution z of F , we used a subroutine of the so-called Bisolve algorithm from [BEK+13BEK+13;KS15KS15], which currently constitutes one of the fastest exact and complete algorithm for solvingbivariate systems.
Bisolve is a classical elimination approach that projects the solutions ofthe system on each of the two coordinate axis in a first step by means of resultant computationand root isolation. This yields a set of points on a two-dimensional grid that are all possiblecandidates for the solutions of the system. Also, the candidates can be approximated toan arbitrary precision using root refinement for univariate polynomials. Then, in a secondstep, in order to check whether a certain candidate is a solution or not,
Bisolve combinesinterval arithmetic and an inclusion test based on bounds on the cofactors g and g in therepresentation R = g · f + g · g of the resultant polynomials R as an element in the ideal (cid:104) f , f (cid:105) . This inclusion test is similar to our approach proposed in this paper, however, notruncation of the original system is considered. Also, it is tailored to the bivariate case anddoes not yield the multiplicity of a solution. In our experiments, we replaced the originalinclusion test in the Bisolve algorithm by
PolySol and compared the precision demandand the running time to that of the original variant. We observed that, for a multiplicity k of z that is small in comparison to the degrees of the input polynomials, our novel approachoutperforms the original variant. At least, for the considered instances, we observed a sub-linear dependency of the needed precision on the degrees of the input polynomials. Noticethat this is not in line with the derived bounds on the precision demand, which suggest atleast a quadratic dependency. However, we remark that the given bounds are just worst-case bounds. In addition, our experiments can only be considered as preliminary at thecurrent time, nevertheless we are confident that future work on this topic will support ourfirst impressions. Related Work.
The literature on solving zero-dimensional polynomial systems is vast andwe can only give an incomplete overview. A historical summary and an overview of knowntechniques can be found in [Laz09Laz09] and [DE06DE06], respectively.There are roughly two different classes of methods – numeric and symbolic methods.To the best of our knowledge, all existing complete and certified algorithms are based onelimination techniques. Using Gröbner bases [Buc06Buc06; Fau02Fau02] or resultants, they reduce theproblem of solving a multivariate system to the problem of computing the roots of a univariatepolynomial. Such methods further allow us to compute the coordinates of all solutions in termsof rational functions in the roots of a univariate polynomial (also called Rational UnivariateRepresentation). A corresponding implementation [Rou99Rou99] has proven to be quite efficient forsystems of moderate size. Also, these methods are well understood in theory and correspondingcomplexity bounds are available [BS16BS16]. The major drawback of these methods is that the costfor the considered symbolic operations becomes prohibitive for larger systems. In contrast,numerical methods, e.g. based on subdivision techniques or homotopy continuation, oftenallow us to compute good approximations of the solutions. Unfortunately, they typically failto give guarantees on the correctness of the computed results.5ne classical numeric approach is Newton’s method, see [Rum10Rum10, Section 13] for a gen-eral description and an approach that uses Newton’s iteration with interval arithmetic. Shuband Smale introduced α -theory [Blu98Blu98], where they provide conditions on a simple solutionsuch that Newton iteration is guaranteed to yield quadratic convergence. Recent work [HL17HL17]uses Newton iteration and α -theory to verify the existence of simple solutions of systems ofpolynomial-exponential equations, however, the approach does not extend to multiple solu-tions. In [Zhi17Zhi17], an extension of α -theory is introduced that allows us to also certify multiplesolutions of a polynomial system in a “numerical fashion” as studied in this paper.Another very popular numeric approach are homotopy continuation methods. There hasbeen also quite some implementation effort, see PHCpack [Ver99Ver99] and Bertini [BHS+13BHS+13]. Inparticular, we want to mention the work by Verschelde and Haegemans [VH94VH94]. From a high-level point of view, their approach is similar to ours as it is also based on Rouché’s theorem.Their method relies on finding a sparse part of the polynomial system that dominates the restof the system on the border of a considered region and can be used as a better starting systemfor homotopy based techniques. The main differences to our approach are the following. First,we use our technique to directly certify the existence of a zero, not only in order to constructa starting system for a numerical method. Moreover, the system that we use in order toapproximate the input system is of lower degree, more precisely our “dominating part” isalways of degree k if k is the multiplicity of the zero in the given region. In contrast to theirresult, we also show that the precision that is needed in order to do so directly depends onthe arrangements of the zeros of the system. Van der Hoeven [Hoe11Hoe11] describes methods fortracking homotopy paths in a certified manner. Using an analytic variant of the geometricresolution method [GHM+95GHM+95].Subdivision methods [MP09MP09; BCG+08BCG+08] are usually incomplete in the sense that they onlyprovide exclusion predicates and lack inclusion predicates. Thus they can be used in order tocompute regions that are guaranteed to be free of solutions to the system but cannot ultimatelyguarantee that a region contains a zero. We want to stress that our work now provides aninclusion predicate that could be included in these approaches in order to turn these methodsinto complete methods.
We start by introducing frequently used notation and important definitions.1. For a point x = ( x , . . . , x n ) ∈ C , we define the norm (cid:107) x (cid:107) of x to be the ∞ -norm bydefault, that is, (cid:107) x (cid:107) := (cid:107) x (cid:107) ∞ = max i ∈ [ n ] | x i | . In addition, we define M ( x ) := max(1 , (cid:107) x (cid:107) ) and log( x ) = M (log( M ( x ))) .2. For a polynomial f = (cid:80) α c α x α ∈ C [ x ] , we define deg( f ) := max α =( α ,...,α n ): c α (cid:54) =0 α + . . . + α n Note that it is not strictly necessary to know this parameter k , since a binary search for k can find a goodenough approximation.
6o be the (total) degree of f . The norm of f is defined as (cid:107) f (cid:107) := (cid:107) f (cid:107) ∞ = max α | c α | . We further define τ f := (cid:100) log (cid:107) f (cid:107)(cid:101)
3. For a polynomial system F = ( f , . . . , f n ) , with f i ∈ C [ x ] of total degree d i , we define d F := max i ∈ [ n ] d i , D F := (cid:89) i ∈ [ n ] d i , and τ F := max i ∈ [ n ] τ f i .D F is also called the Bézout bound in the literature. It constitutes an upper bound onthe total number of solutions (counted with multiplicities) of a zero-dimensional system F . For a system F with generic coefficients, it actually equals the number of solutions.4. We further say that a polynomial f = (cid:80) α c α x α ∈ Z [ x ] with integer coefficients has magnitude ( d, τ ) if d f ≤ d and τ f ≤ τ . A system F = ( f , . . . , f n ) with f i ∈ Z [ x ] hasmagnitude ( d, τ ) if each polynomial has magnitude ( d, τ ) .5. For a polynomial f = (cid:80) α c α x α ∈ C [ x ] and a positive integer κ , we say that φ = (cid:80) α ˜ c α x α is an (absolute) κ -bit approximation of f if each ˜ c α is a dyadic number of the form ( m + m (cid:48) · i ) · − ( κ +1) ∈ Q + i · Q , with m, m (cid:48) ∈ Z , and (cid:107) f − φ (cid:107) ≤ − κ . In other words,each ˜ c α approximates c α to κ bits after the binary point.6. For z ∈ C n and a polynomial f ∈ C [ x ] , we define f [ z ]( x ) := f ( x + z ) = (cid:88) α c α ( x + z ) α to be the shift of f to z . For a system F = ( f , . . . , f n ) , we define the shift of F to z as F [ z ] = ( f [ z ] , . . . , f n [ z ]) .7. For k ∈ [ d ] , we denote with f ≤ k := (cid:80) α : | α |≤ k c α x α the truncation of f of degree k .For a system F = ( f , . . . , f n ) , we define the truncation of F of degree k as F ≤ k =( f ≤ k , . . . , f n ≤ k ) . We first collect some bounds on the size of | f ( z ) | and (cid:107) f [ z ] (cid:107) depending on the modulus ofsome point z ∈ C n and the norm (cid:107) f (cid:107) of some polynomial f ∈ C [ x ] . We also give bounds onthe error that occurs when computing f ( z ) or f [ z ] not exactly at z but at a nearby point ζ . Lemma 1.
Let f ( x ) = (cid:80) α c α x α ∈ C [ x ] = C [ x , . . . , x n ] of total degree d and with (cid:107) f (cid:107) ≤ τ .Moreover, let k ∈ [ d ] = { , . . . , d } , z ∈ C n , and ζ be an approximation of z with (cid:107) ζ − z (cid:107) < − L ,then it holds:(a) | f ≤ k ( z ) | ≤ (cid:0) n + kk (cid:1) · τ · M ( z ) k , and in particular | f ( z ) | < (cid:0) n + dd (cid:1) · τ · M ( z ) d .(b) If (cid:107) z (cid:107) ≤ , then | f ( z ) − f ≤ k ( z ) | ≤ (cid:107) z (cid:107) k +1 · (cid:0) n + dd (cid:1) · τ .(c) | f ( ζ ) − f ( z ) | ≤ τ − L · (cid:0) n + dd (cid:1) · M ( z ) d . d) (cid:107) f [ z ] (cid:107) < d n · (cid:0) n + dd (cid:1) · τ · M ( z ) d and (cid:107) f [ ζ ] − f [ z ] (cid:107) < τ − L · d n · (cid:0) n + dd (cid:1) · M ( z ) d .Proof. Part (aa) and (bb) follow immediately from the fact that f ≤ k has at most (cid:0) n + kk (cid:1) coefficientsand each occurring term c α · z α has absolute value bounded by τ · (cid:107) z (cid:107) | α | . Part (cc) is a directconsequence of [MOS11MOS11, Theorem 12], which provides general bounds on the error whenevaluating a multivariate polynomial using floating point computation. For the last claim,notice that f [ z ] = (cid:88) α ∈ Z n ≥ ∂ α f ( z ) α ! x α and f [ ζ ] = (cid:88) α ∈ Z n ≥ ∂ α f ( ζ ) α ! x α , where α = ( α , . . . , α n ) , α ! := α ! · · · α n ! , and ∂ α f := ∂ | α | f∂x α ··· ∂x αnn . The polynomials ∂ α fα ! havetotal degree bounded by d and their norm is upper bounded by τ · d n = 2 τ + n log d . Hence,Part (aa) implies the first part of (dd). The second part follows from Part (cc) because, for any α , it holds that | ∂ α f ( z ) α ! − ∂ α f ( ζ ) α ! | ≤ τ − L · d n · (cid:18) n + dd (cid:19) · M ( z ) d . We further provide the following lemma that investigates the influence of considering onlyan approximation of a polynomial f when looking at shift and truncation. Lemma 2.
Let f ( x ) ∈ C [ x ] be a polynomial of total degree d with norm (cid:107) f (cid:107) ≤ τ , and let z ∈ C n and ζ such that (cid:107) ζ − z (cid:107) < − L . Furthermore, let φ be an approximation of f [ ζ ] ≤ k oftotal degree at most k , with k ∈ [ d ] , such that (cid:107) φ − f i [ ζ ] ≤ k (cid:107) ≤ − ( k +1) L . Then, for any x with (cid:107) x (cid:107) ∈ [2 − L , , it holds | φ ( x ) − f [ ζ ]( x ) | ≤ (cid:107) x (cid:107) k +1 · d n τ +1 [ M ( ζ ) · ( n + d ) ] d . Proof.
We first observe that using the triangle inequality, simple bounds on the number ofmonomials of lower ( ≤ k ) and higher ( ≥ k + 1 ) degree, and the fact that (cid:107) x (cid:107) ≤ yields | φ ( x ) − f [ ζ ]( x ) | ≤ | f [ ζ ]( x ) − f [ ζ ] ≤ k ( x ) | + | φ ( x ) − f [ ζ ] ≤ k ( x ) |≤ (cid:107) x (cid:107) k +1 · (cid:107) f [ ζ ] − f [ ζ ] ≤ k (cid:107) · (cid:18) n + dd (cid:19) + (cid:107) φ − f [ ζ ] ≤ k (cid:107) · (cid:18) n + kk (cid:19) . Then, applying Lemma 11 part (dd) to the left summand and the condition on the approximation (cid:107) φ − f [ ζ ] ≤ k (cid:107) ≤ − ( k +1) L , we conclude that | φ ( x ) − f [ ζ ]( x ) | ≤ (cid:107) x (cid:107) k +1 · d n · (cid:18) n + dd (cid:19) · τ · M ( ζ ) d + 2 − ( k +1) L · (cid:18) n + kk (cid:19) ≤ (cid:107) x (cid:107) k +1 · [ d n · τ · ( n + d ) d · M ( ζ ) d + ( n + d ) d ] ≤ (cid:107) x (cid:107) k +1 · d n τ +1 [ M ( ζ ) · ( n + d ) ] d , where the second to last inequality follows from (cid:107) x (cid:107) ≥ − L .In our algorithm, we will consider a transformation of the coordinate system induced by arotation x (cid:55)→ S · x , where S ∈ SO( n ) is a rotation matrix with rational entries. The followinglemma quantifies the impact of such a rotation on the bit-size of the coefficients of a givenpolynomial f . 8 emma 3. Let f = (cid:80) α c α · x α ∈ C [ x ] be a polynomial of total degree d and S ∈ SO( n ) be arotation matrix. Then, f ∗ := f ◦ S , it holds that (cid:107) f ∗ (cid:107) ≤ τ F · (cid:0) n + dd (cid:1) .Proof. Notice that each of the entries a r,s of the rotation matrix S = ( a r,s ) r,s has absolutevalue at most . Thus, f ∗ ( x ) = f ◦ S ( x ) = (cid:80) α : c α c α · [( a x + · · · + a ,n x n ) α · · · ( a n x + · · · + a n,n x n ) α n ] has coefficients of absolute value bounded by τ F · (cid:0) n + dd (cid:1) as, when expandingthe product ( a x + · · · + a ,n x n ) α · · · ( a n x + · · · + a n,n x n ) α n for a fixed α , there can be atmost (cid:0) n + dd (cid:1) terms contributing to a specific monomial x α (cid:48) . Let us assume that an arbitrary zero-dimensional system F = ( f , . . . , f n ) as in (11) is given.That is, f i has total degree d i , (cid:107) f i (cid:107) < τ i for all i , and it is assumed that the total numberof solutions of F = 0 , also at “infinity” (see the considerations below for an explanation), isfinite. We now briefly describe the so-called hidden-variable approach that allows us to projectthe zeros of the system on an arbitrary coordinate axis. For more details, we recommend theexcellent textbook [CLO05CLO05] by Cox, Little, and O’Shea.In a first step, we consider a homogenization of the system, that is, we introduce anadditional (homogenizing) variable x n +1 and multiply each occurring term in each f i witha suitable power of x n +1 such that the so obtained polynomials f hi ∈ C [ x , . . . , x n +1 ] arehomogenous and of total degree d i , respectively; see also the example below. Notice that eachsolution ( x , . . . , x n ) ∈ C of F = 0 yields a solution ( x , . . . , x n , of the homogenized system F h : f h ( x , . . . , x n +1 ) = . . . = f hn ( x , . . . , x n +1 ) = 0 (3)In addition, if ( x , . . . , x n +1 ) ∈ C n +1 is a solution of F h = 0 , then ( t · x , . . . , t · x n +1 ) is a zeroof F h for all t ∈ C . In particular, if x n +1 (cid:54) = 0 , we can set t = 1 /x n +1 , which yields the solution ( x /x n +1 , . . . , x n /x n +1 ) of F = 0 . It is thus preferable to consider the set S of solutions of theabove homogenized system as a set of points in the n -dimensional projective space P n . Theset S then decomposes into the set S < ∞ = { ( x : . . . : x n +1 ) ∈ S : x n +1 = 1 } of so-called affinesolutions , for which x n +1 = 1 , and the set S ∞ = { ( x : . . . : x n +1 ) : x n +1 = 0 } of solutionsat infinity , for which x n +1 = 0 . Notice that there is a one-to-one correspondence between theaffine solutions of the homogenized system and the solutions of the original system (11).As mentioned above, we aim to compute the projections of the solutions of F = 0 on oneof the coordinate axis, say w.l.o.g., x = x . For this, suppose that we fix some value ξ for x .Plugging x = ξ into the initial system then yields the specialized system F [ ξ ] : f [ ξ ]1 ( x , . . . , x n ) = . . . = f [ ξ ] n ( x , . . . , x n ) = 0 , with f [ ξ ] i ( x , . . . , x n ) := f i ( ξ, x , . . . , x n ) and the corresponding homogenized system ( F [ ξ ] ) h : ( f [ ξ ]1 ) h ( x , . . . , x n +1 ) = . . . = ( f [ ξ ] n ) h ( x , . . . , x n +1 ) = 0 , (4)where ( f [ ξ ] i ) h denotes the homogenization of f [ ξ ] i . Notice that, in general, ( f [ ξ ] i ) h does notequal ( f hi ) [ ξ ] = f hi ( ξ, x , . . . , x n +1 ) , that is, we cannot deduce the system in (44) from plugging ξ into the homogenized system in (33). The reason is that the total degree of f i may becomesmaller for certain values for ξ , and thus homogenization does not commute with specialization.9 xample. For f := x x − x + x x + x and ξ = 2 , we have f h = x x − x x + x x x + x x , f [ ξ ] = f (2 , x , x ) = 2 x + x , and ( f [ ξ ] ) h = 2 x x + x , which does not equal ( f h ) [ ξ ] = f h (2 , x , x , x ) = 2 x − x x + 2 x x + x x .You may notice that (44) is a polynomial system consisting of n homogenous polynomialsin n variables. If the initial homogenized system had a solution with x = ξ , then this wouldyield a solution of (44) and vice versa. In other words, ξ would be the projection of a solution ofthe initial system. The following important result now gives a necessary and sufficient criteriato check whether this is actually the case. Theorem 2 ([CLO05CLO05], Chapter 3, Theorems 2.3 and 3.1) . Let G be a system of n homogeneouspolynomials in n variables of total degrees d , . . . , d n . Then, there is a unique polynomial Res( G ) = Res d ,...,d n ∈ Z [ u ] in the coefficients u of G if and only if G ( x ) = 0 has a non-trivialsolution x ∈ P n − . Res( G ) is homogeneous in the variables of f i of degree d · · · d i − · d i +1 · · · d n and its total degree equals (cid:80) ni =1 d · · · d i − · d i +1 · · · d n .Example. The homogeneous system G : ax + bx x + cx = dx + ex = 0 (with generalcoefficients a, b, c, d , and e ) has a solution in P if and only if the involved coefficients fulfillthe equality Res( G ) = Res , = ae − be + cd = 0 .For an arbitrary polynomial system G consisting of n + 1 (not necessarily homogenous)polynomials in C [ x , . . . , x n ] , we simply define Res( G ) = Res( G h ) . Since G has the samecoefficients as G h , it still holds that Res( G ) is a polynomial in the coefficients of G . In addition,since there is a one-to-one correspondence between the solutions of G and the affine solutionsof G h , it follows that Res( G ) = 0 if and only if G h = 0 has a solution in P n .Now, in order to compute all values ξ such that there exists a solution ( x , . . . , x n ) of ourinitial system F = 0 with x = ξ , we aim to apply the above theorem to the system as definedin (44), however we now consider ξ as an indeterminate (so called hidden variable ) rather thana fixed value. There are some subtleties with this approach. In particular, the degrees of thepolynomials f [ ξ ] i may be different for certain values for ξ , which is crucial as the definitionof the resultant polynomial Res strongly depends on the degrees of the given polynomials.However, we can avoid such critical situations if we assume that the given polynomials f i fulfill some mild prerequisites. Lemma 4.
Suppose that each polynomial f i contains a term of total degree d i that does notdepend on x and write f i = (cid:88) α =( α ,...,α n ) c i,α ( x ) · x α · · · x α n n ∈ C [ x ][ x , . . . , x n ] as a polynomial in x , . . . , x n with coefficients c i,α ∈ C [ x ] . Furthermore, let F i := (cid:88) α =( α ,...,α n ) c i,α ( x ) · x α · · · x α n n · x d i − α − ... − α n n +1 be its corresponding homogenization (with respect to the variables x , . . . , x n ), then it holds:(a) For all ξ ∈ C , we have ( f [ ξ ] i ) h = F [ ξ ] i and ( f [ ξ ] i ) h has total degree d i . We remark that
Res only depends on the actual degrees of the polynomials. b) Each root x = ξ ∈ C of R ( x ) := Res( F , . . . , F n ) yields a solution ( x , . . . , x n ) ∈ C n of F = 0 with x = ξ and vice versa.Proof. Part (a) follows directly from the fact that the total degree of f [ ξ ] i is equal to d i forall ξ as there exists a term of degree d i that does not depend on ξ . For (b), we first remarkthat the resultant of the polynomials F i is a polynomial in the coefficients of the F i , and thusa polynomial in x . Since the degree of each f i does not depend on the choice of x = ξ , wealso have R ( ξ ) = Res( F | x = ξ , . . . , F n | x = ξ ) . Now, let x = ξ be a complex root of R , thenaccording to Theorem 22, there must exist a solution ( ξ : . . . : ξ n +1 ) ∈ P n − of the system ( F i | x = ξ ) i =1 ,...,n . In order to prove that this solution is an affine solution (i.e. a solution of F ),we assume for contradiction that ξ n +1 = 0 . Plugging x n +1 = 0 into the polynomials F i yields F i | x n +1 =0 = (cid:88) α =( α ,...,α n ): α + ... + α n = d i c i,α ( x ) · x α · · · x α n n . Hence, each of the terms c i,α ( x ) occurring in the above sum is a constant that does notdepend on x . Since ( ξ : ξ : . . . : ξ n ) is a solution of the system F i | x n +1 =0 , we concludethat ( x : ξ : . . . : ξ n ) is a solution of F i | x n +1 =0 for any x . This contradicts our assumptionthat F has only finitely many solutions. It follows that ( ξ /ξ n +1 , . . . , ξ n /ξ n +1 , is a solutionof ( F i | x = ξ ) i =1 ,...,n , and thus ( ξ, ξ /ξ n +1 , . . . , ξ n /ξ n +1 ) is a solution of F = 0 . For the otherdirection, let ( ξ , . . . , ξ n ) be a solution of F = 0 , then ( ξ : . . . : ξ n : 1) is an affine solutionof the corresponding homogenized system, and thus ( ξ : . . . : ξ n : 1) a solution of the system ( F i | x = ξ ) i =1 ,...,n = 0 . This implies that R ( ξ ) = 0 .Obviously, the above considerations apply for any coordinate (hidden-variable) x k ontowhich we aim to project the solutions. The corresponding resultant polynomial Res( F , x k ) ∈ C [ x k ] is called the hidden-variable resultant with respect to x k . The following theorem [BS16BS16]bounds the cost for computing the hidden-variable resultant in the special case where thepolynomials f i have integer coefficients. The technique is based on a method due to Emirisand Pan [EP05EP05] and an asymptotically fast algorithm for determinant computation due toStorjohann [Sto05Sto05]. Theorem 3 ([BS16BS16, Prop. 1]) . Let F = ( f i ) i =1 ,...,n be a polynomial system with integerpolynomials f i ∈ Z [ x , . . . , x n ] of magnitude ( d, τ ) . There is a Las-Vegas algorithm to compute Res( F , x k ) in an expected number of bit operations bounded by ˜ O ( n ( n − ω +1) ( d + τ ) d ( ω +2) n − ω − ) . We further remark that a root ξ of Res( F , x k ) might origin from several solutions z =( z , . . . , z n ) of F = 0 sharing the same x k -coordinate x k = ξ . Under the requirements fromLemma 44, it holds that the multiplicity of ξ as a root of Res( F , x k ) equals the sum of themultiplicities of all these solutions z . Also, the roots of Res( F , x k ) are exactly the projectionsof the finite solutions onto the x k -coordinate, and vice versa. Furthermore, if there no solutionat infinity, then Res( F , x k ) has degree D F as the system has exactly D F solutions (countedwith multiplicity), which are all finite, and the roots of Res( F , x k ) are exactly the projectionsof these solutions onto the x k -coordinate. ω denotes the exponent in the complexity of matrix multiplication. The current record bound for ω is ω ≤ . according to [Gal14Gal14] emma 5. Let F = ( f i ) i =1 ,...,n , with f i = (cid:80) α : | α |≤ d i c i,α · x α ∈ C [ x ] of total degree d i , be apolynomial system in n variables x = ( x , . . . , x n ) with general coefficients c i,α . Consider thedecomposition f i ( x ) = (cid:88) α : | α | = d i c i,α · x α (cid:124) (cid:123)(cid:122) (cid:125) := f i,di ( x ) + (cid:88) α : | α | 0) = f i,d i , and a generic system of n homogenouspolynomials in n variables has no solution. Thus, there exists no solution at infinity, whichalso rules out the possibility of the system being non zero-dimensional. Now, suppose that thecoefficients are generically chosen such that all solutions are finite. Then, the total number ofsolutions equals the Bézout number D F and the degree of Res( F , x k ) equals D F . Accordingto Theorem 22, LC(Res( F , x k )) ∈ Z [ c i,α ] is a polynomial in the coefficients c i,α . Now, if LC(Res( F , x k )) would depend on some coefficient c i,α with | α | < d i , then, for generic choiceof all other coefficients, we could choose such a c i,α in a way such that the leading coefficientbecomes zero, and thus deg Res( F , x k ) < D F , a contradiction. This shows that, for genericchoice of the coefficients c i,α , the leading coefficient LC(Res( F , x k )) does not depend on thecoefficients of the polynomials f i, Let F = ( f i ) i =1 ,...,n be an arbitrary polynomial system as in Lemma 55 with d i = d for all i , and let ¯ F : ¯ f i := n (cid:88) α : | α | = d +1 c i,α · x α + f i , with c i,α ∈ Z for all α with | α | = d + 1 , be the system obtained by adding polynomials of the form (cid:80) nα : | α | = d +1 c i,α · x α to each f i . If ¯ F does not have any solution at infinity (which is the case for generic choice of the coefficients c i,α ), then it holds that LC(Res( ¯ F , x k )) ∈ Z (cid:54) =0 . Proof. If ¯ F has no solution at infinity, then ¯ F is zero-dimensional and, in addition, Res( ¯ F , x k ) has degree D ¯ F = ( d + 1) n . From Lemma 55, we further conclude that LC(Res( ¯ F , x k )) onlydepends on the coefficients c i,α of the degree ( d + 1) -parts ¯ f i,d +1 of the polynomials ¯ f i . Hence,we have LC(Res( ¯ F , x k )) ∈ Z (cid:54) =0 . Example: Let ¯ f i = (cid:80) nj =1 a ij · x d +1 j + f i with f i ∈ C [ x ] ≤ d polynomials of total degree at most d . Then, it holds that Res( ¯ F , x k ) = ± det( a ij ) · x ( d +1) n k + · · · det( a i,j ) (cid:54) = 0 , then ¯ F has no solution at infinity as each such solution would yielda non-trivial solution of the linear system (cid:80) nj =1 a i,j · X j = 0 . Thus, ¯ F is zero-dimensionalin this case and Res( ¯ F , x k ) has degree D ¯ F = ( d + 1) n . From Lemma 55, we further concludethat LC(Res( ¯ F , x k )) only depends on the coefficients a i,j of the degree ( d + 1) -parts ¯ f i,d +1 ofthe polynomials ¯ f i . Hence, we have LC(Res( ¯ F , x k )) = LC(Res( ¯ f , = d +1 , . . . , ¯ f n, = d +1 , x k )) , andusing Theorem 2.3 and Theorem 3.5 in [CLO05CLO05] further shows that Res( ¯ f , = d +1 , . . . , ¯ f n, = d +1 , x k ) = det( a i,j ) ( d +1) n · Res(( x d +11 , . . . , x d +1 n ) , x k )= ± det( a i,j ) ( d +1) n · x ( d +1) n k . It is also well known (e.g. this follows from Theorem 44 below) that Res( F , x k ) is containedin the ideal I := (cid:104) f , . . . , f n (cid:105) defined by the polynomials f , . . . , f n . In particular, for polyno-mials f i ∈ Z [ x , . . . , x n ] with integer coefficients, this guarantees the existence of an integer λ ,with λ (cid:54) = 0 , and polynomials g i ∈ Z [ x , . . . , x n ] with λ · Res( F , x k ) = g · f + · · · g n · f n . (5)Recent work [DKS13DKS13] allows us to bound the magnitude of the polynomials g i as well as thesize of λ . For this, we first write f i = (cid:80) α c i,α ( x k ) x α (cid:54) = k as a polynomial in x (cid:54) = k with coefficients c i,α ∈ Z [ x k ] , where x (cid:54) = k denotes all but the k ’th variable. We further introduce a variable u i,α for every coefficient polynomial c i,α . Let u i = ( u i,α ) α be the variables corresponding tothe polynomial f i , and let u = ( u , . . . , u n ) denote the variables for all polynomials. Then, F can be considered as a system consisting of n polynomials in n − variables x (cid:54) = k withcoefficients u . Thus, its resultant Res( F ) is a polynomial in Q [ u ] , which is further containedin the ideal (cid:104) f , . . . , f n (cid:105) ⊂ Q [ u , x (cid:54) = k ] . The following theorem, which is a consequence ofTheorem 4.28 in [DKS13DKS13] (see also [DKS13DKS13, pp. 6]), gives bounds on the degree and height ofthe polynomials in the cofactor-representation of Res( F ) in this ideal. Theorem 4 ([DKS13DKS13] Consequence of Theorem 4.28) . Given a polynomial system Φ =( ϕ , . . . , ϕ n +1 ) with ϕ i = (cid:80) α u i,α x α ∈ Z [ u , x ] of total degree d i in x = ( x , . . . , x n ) . Then,for any k ∈ [ n ] , there exists a λ ∈ Z (cid:54) =0 and polynomials γ i ∈ Z [ u , x ] such that λ · Res(Φ) = (cid:88) i ∈ [ n +1] γ i · ϕ i , deg u j ( γ i ϕ i ) ≤ (cid:89) (cid:96) (cid:54) = j d (cid:96) and τ λ Res(Φ ,x k ) , τ γ i ≤ (6 n + 10) log( n + 3) D Φ for j ∈ [ n ] and i ∈ [ n + 1] , where τ p denotes the bit-size of a polynomial p ∈ Z [ u , x ] . We can now derive bounds on the degree and the bit-sizes of the polynomials g i as well ason the bit-size of λ in (55) from the above theorem: Corollary 2. Given a zero-dimensional polynomial system F = ( f , . . . , f n ) with polynomials f i ∈ C [ x ] , we can explicitly compute (see (66) and (77)) positive integers A F and B F , with F = ˜ O (cid:0) nD F (cid:1) and B F = ˜ O ( n · D F + τ F · max i D F d i ) , such that there exists an integer λ ∈ Z (cid:54) =0 and polynomials g i ∈ C [ x ] with | λ | ≤ A F , deg g i ≤ D F , τ g i ≤ B F , and λ · Res( F , x k ) = n (cid:88) i =1 g i · f i . If all polynomials f i have only integer coefficients, then we may further assume that the poly-nomials g i have only integers coefficients as well.Proof. For each i ∈ [ n ] , write f i ( x ) = (cid:80) α u i,α x α (cid:54) = k as a polynomial in the variables x (cid:54) = k andwith coefficients u i,α ∈ C [ x k ] . Theorem 44 now guarantees the existence of a positive λ ∈ Z (cid:54) =0 and polynomials g i ∈ Z [ u , x (cid:54) = k ] with λ · Res( F , x k ) = (cid:80) i ∈ [ n ] g i · f i . Notice that since u only depends on x k , we may consider each g i as an element in C [ x ] . In addition, we have deg u j ( g i f i ) ≤ (cid:81) (cid:96) (cid:54) = j d (cid:96) , and thus deg x ( g i ) ≤ deg x ( g i f i ) ≤ D F = d · · · d n as each u j,α hasdegree bounded by d j . We can now write each polynomial g i as g i = (cid:80) α : | α |≤ D F P i,α ( u ) · x α (cid:54) = k with polynomials P i,α = (cid:80) β =( β ,...,β N ): | β |≤ D F /d i c i,α,β · u β . From Theorem 44, we conclude that c i,α,β are integers of absolute value | c i,α,β | < A F , where A F := (cid:100) (6 n + 4) log( n + 2) D F (cid:101) = O (cid:0) nD F log( n ) (cid:1) . (6)In addition, N ≤ (cid:80) ni =1 (cid:0) d i + nd i (cid:1) ≤ n · (cid:0) d F + nd F (cid:1) ≤ n ( d F + n ) n denotes the number of distinctcoefficients u i,α . Further notice that, for each β , u β is a product of at most D F univariatepolynomials in C [ x k ] , each of degree at most d F and of norm bounded by τ F . Hence, it canbe written as a sum of at most ( d F + 1) D F terms, each of absolute value at most τ F · D F /d i .We conclude that the norm of g i is bounded by (cid:18) D F + ND F (cid:19) · ( d F + 1) D F · A F · τ F · D F /d i ≤ [( n ( d F + n ) n + D F ) · ( d F + 1)] D F · A F · τ F · D F /d i ≤ [( n + 1)( d F + n ) n +1 ] D F · A F · τ F · D F /d i ≤ [( d F + n ) n +4 · A F · τ F · D F /d i ≤ B F , where we define B F := 2 · (cid:100) D F · (6 n + 4) log( d F + n ) (cid:101) + τ F · max i D F d i = ˜ O ( n · D F + τ F · max i D F d i ) . (7)The final claim follows from the fact that f i ∈ Z [ x ] for all i implies that u i,α ∈ Z [ x k ] for all i, α , and thus g j ∈ Z [ x ] for all j . In the previous subsection, we have outlined how to project the solutions of a polynomial ontoone of the coordinate axis. One subtlety of the approach was that certain mild conditions onthe input polynomials need to be fulfilled in order to guarantee that the roots of the hiddenvariable resultant are exactly the projections of the (finite) solutions of the initial system; seeLemma 44. Another drawback of the approach is that distinct solutions might be projected14nto the same point or onto two very nearby points on the coordinate axis, that is, the actualdistance between distinct solutions is no longer preserved after the projection. We will showhow to address these issues by using a random rotation of the coordinate system. We firststart with the special case of dimension . Lemma 6. Let p (cid:96) = ( x (cid:96) , y (cid:96) ) ∈ C be N points such that (cid:107) p (cid:96) (cid:107) (cid:54) = 0 for all (cid:96) = 1 , . . . , N . Let k be chosen uniformly at random from [2 L ] . Then, with probability at least − N L , for eachpoint p (cid:48) (cid:96) = (cid:18) x (cid:48) (cid:96) y (cid:48) (cid:96) (cid:19) := S k ( L ) · (cid:18) x (cid:96) y (cid:96) (cid:19) , with S k ( L ) := − ( k · − L ) k · − L ) − · ( k · − L )1+( k · − L ) · ( k · − L )1+( k · − L ) − ( k · − L ) k · − L ) ∈ SO(2) , it holds that min( | x (cid:48) (cid:96) | , | y (cid:48) (cid:96) | ) > − ( L +2) · (cid:107) p (cid:96) (cid:107) for all (cid:96) .Proof. Notice that each matrix S k ( L ) is a rotation matrix with respect to the angle φ k ∈ [0 , π/ with cos φ k = − ( k · − L ) k · − L ) and sin φ k = · ( k · − L )1+( k · − L ) . We further note that the function h ( t ) = ( − t t , t t ) describes the trace of a point on the quarter-circle. Moreover, we have ˙ h ( t ) = ( − t (1+ t ) , − t +2(1+ t ) ) , and since | ˙ h ( t ) | = t is a decreasing function in t , it follows that thedifference between two consecutive angles φ k +1 and φ k is decreasing in k . We thus concludethat all differences are lower bounded by φ L − φ L − = (cid:82) − − L t dt ≥ (cid:82) − − L dt = 2 − L .Now, let L k ⊂ R be the line passing through the origin and the point (cos φ k , sin φ k ) , and let L ⊥ k ⊂ R be line that passes through the origin and is orthogonal to L k . In addition, for eachpoint p (cid:96) = ( (cid:60) ( x (cid:96) ) + i · (cid:61) ( x (cid:96) ) , (cid:60) ( y (cid:96) ) + i · (cid:61) ( y (cid:96) )) , we define ¯ p (cid:96) = (cid:40) ( (cid:60) ( x (cid:96) ) , (cid:60) ( y (cid:96) )) if (cid:60) ( x (cid:96) ) + (cid:60) ( y (cid:96) ) ≥ (cid:61) ( x (cid:96) ) + (cid:61) ( y (cid:96) ) ( (cid:61) ( x (cid:96) ) , (cid:61) ( y (cid:96) )) otherwise . . Then, ¯ p (cid:96) is a point in R with (cid:107) ¯ p (cid:96) (cid:107) ≥ (cid:107) p (cid:96) (cid:107) / √ . Let ∆ (cid:96) ⊂ R be the disc centered at ¯ p (cid:96) ofradius r (cid:96) = 2 − L − · (cid:107) p (cid:96) (cid:107) . Let q, r ∈ ∆ (cid:96) be any two points in ∆ (cid:96) and α be the angle at theorigin of the triangle given by the origin and the points q and r . Then, it holds that α ≤ · arctan (cid:18) − L − · (cid:107) p (cid:96) (cid:107)(cid:107) p (cid:96) (cid:107) / √ (cid:19) < − L − ) < · − L − ≤ − L . Since the angle between any two distinct lines L k and L k (cid:48) is lower bounded by − L , it thusfollows that there can be at most one k such that L k or L ⊥ k intersects ∆ (cid:96) . Hence, if we picka k ∈ { , . . . , L } uniformly at random and choose L k and L ⊥ k as the axis of the coordinatesystem obtained by rotating the initial system by φ k , then, with probability at least − N L ,the new coordinates (¯ x (cid:48) (cid:96) , ¯ y (cid:48) (cid:96) ) of each point ¯ p (cid:96) will meet the condition that min( | ¯ x (cid:48) (cid:96) | , | ¯ y (cid:48) (cid:96) | ) > − L − · (cid:107) p (cid:96) (cid:107) . Hence, the same holds true for the points S k ( L ) · p (cid:96) .We now turn to the general n -dimensional case. For integers k and L and distinct indices15 , j ∈ { , . . . , n } , we define S [ ij ] k ( L ) := · · · · · · · · · ... ... ... ... ... ... ... · · · − ( k · − L ) k · − L ) · · · − · ( k · − L )1+( k · − L ) · · · ... ... ... ... ... ... ... · · · · ( k · − L )1+( k · − L ) · · · − ( k · − L ) k · − L ) · · · ... ... ... ... ... ... ... · · · · · · · · · ∈ SO( n ) , (8)to be a rotation matrix that operates on the i -th and j -th coordinate only. We further definethe set of rotation matrices S N := (cid:89) i,j ∈ [ n ] : i Let N be a positive integer and p (cid:96) ∈ C n be N (cid:48) , with N (cid:48) ≤ N , points such that (cid:107) p (cid:96) (cid:107) (cid:54) = 0 for all (cid:96) = 1 , . . . , N (cid:48) . S N and L are defined as in (99). Then, it holds(a) Choosing integers k ij ∈ [2 L ] for every pair i, j uniformly at random yields, with proba-bility at least / , a rotation matrix S ∈ S N such that, for each point p (cid:48) (cid:96) := S ( L ) · p (cid:96) , itholds that min i | p (cid:48) (cid:96),i | ≥ (2 n N ) − n · (cid:107) p (cid:96) (cid:107) .(b) There is an integer λ of bit-size ˜ O ( n log N ) such that the entries of λS and λS − areinteger numbers of bit-size ˜ O ( n log N ) as well.Proof. The proof follows almost immediately from Lemma 66. Namely, with probability at least − N/ L , both entries p (cid:48) (cid:96),i and p (cid:48) (cid:96),j of each point p (cid:48) (cid:96) := S [ ij ] k ij ( L ) · p (cid:96) will have absolute value atleast − ( L +2) · max( | p (cid:96),i | , | p (cid:96),j | ) . Since at least one of the coordinates of p (cid:96) has absolute value (cid:107) p (cid:96) (cid:107) , we conclude that, with probability (1 − N/ L )( n ) > (1 − N/ L ) n > − n / L /N > / ,each coordinate of each point p (cid:48) (cid:96) = S ( L ) · p (cid:96) has absolute value at least − n · ( L +2) · (cid:107) p (cid:96) (cid:107) ≥ − n · (4(log(2 n N )+1)+2) · (cid:107) p (cid:96) (cid:107) ≥ (2 16 log(2 n N ) ) − n · (cid:107) p (cid:96) (cid:107) = (2 n N ) − n · (cid:107) p (cid:96) (cid:107) . It remains to show the existence of an integer λ of bit-size ˜ O ( n log N ) such that the entries of λS and λS − are of that bit-size as well. Each entry of a matrix S [ ij ] k ij ( L ) is rational number withdenominator L + k ij of bit-size O ( L ) . The matrix S is a product of O ( n ) many such matrices,thus for λ = (cid:81) i,j ∈ [ n ] : i Let F = ( f , . . . , f n ) be a polynomial system as in (11) , S ∈ S D F be a randomlychosen matrix, and let F (cid:48) := F ◦ S − be the corresponding rotated system. Then, with proba-bility larger than / , it holds:(a) For each k ∈ [ n ] , each of the polynomials f i ◦ S − ∈ F (cid:48) contains a monomial of degree d i that does not depend on x k .(b) For each solution z ∈ C n \ of F = 0 , it holds that min i | z (cid:48) i | ≥ (2 n D F ) − n · (cid:107) z (cid:107) , where z (cid:48) = ( z (cid:48) , . . . , z (cid:48) n ) := S · z is the corresponding (rotated) solution of F (cid:48) = 0 .Proof. Since D F constitutes an upper bound on the number of solutions of F = 0 , it followsfrom Lemma 77 (with N = D F ) that, with probability at least / , the inequality in (b) isfulfilled. It thus suffices to prove that, with probability larger than / , the condition in (a)is fulfilled for each coordinate x k . For this, let f ( x ) = (cid:88) α c α x α ∈ C [ x , . . . , x n ] be a polynomial of total degree d , and let S ( k ) − = ( a rs ( k )) rs be the matrix depending on thevalues k := ( k ij ) i,j . Notice that each entry a rs ( k ) is a rational function in k with numeratorsand denominators of total degree (in k ) at most n . Further notice that S ( k ) − maps thepoint (1 , , . . . , to the first column of S ( k ) − and that a full-dimensional subset T of thestrictly positive part { x ∈ R n : (cid:107) x (cid:107) = 1 and x > } of the ( n − -dimensional sphere S n − ⊂ R n is reached via a suitable choice of k ∈ R n . Composing f and S ( k ) − now yields F ( x , k ) = (cid:88) α =( α ,...,α n ) c α · ( a ( k ) · x + · · · + a n ( k ) · x n ) α · · · ( a n ( k ) · x + · · · + a nn ( k ) · x n ) α n , and the coefficient C ( k ) of the monomial x d is thus given by C ( k ) = (cid:88) α : | α | = d c α · a ( k ) α · · · a n ( k ) α n . We first argue that C ( k ) does not vanish identically. Let ˆ f := (cid:80) α : | α | = d c α · x α · · · x α n n be thecorresponding homogenous polynomial of degree d such that ˆ f ( a ( k ) , . . . , a n ( k )) = C ( k ) .Assume that C ( k ) = 0 for all k , then this implies that ˆ f vanishes on each point in T . Sincethe vanishing set of any non-zero homogenous polynomial in n variables has dimension atmost n − , we conclude that ˆ f is the zero-polynomial, and thus c α = 0 for all coefficients of ˆ f . This contradicts our assumption on f . 17ence, it follows that C ( k ) is a non-zero rational function in k . In addition, each term c α · a ( k ) α · · · a n ( k ) α n has a numerator of total degree at most n d in k and a denominatorof the form (cid:81) i,j (2 L + k i,j ) e i,j , with e i,j ∈ N , of degree at most n d in k . This showsthat C ( k ) can be written as a rational function in k of total degree n d + n d ≤ n d as (cid:81) ni,j =1: i For a given polynomial system F we say that a rotationmatrix S ∈ S D F is admissible with respect to F if the statements (aa) and (bb) from Lemma 88hold. We further denote by S F := { S ∈ S D F : S is admissible with respect to F } ⊂ S D F the set of admissible matrices with respect to F . Notice that, even though it is difficult (probably as difficult as computing all solutions of F ) to determine whether a certain matrix in S D F is admissible with respect to F , the previouslemma shows that at least half of the matrices in S D F are admissible. We first sketch our algorithm PolySol and then prove its correctness. We refer the readerto the pseudo-code in Algorithm 11 for details regarding PolySol . The algorithm can beroughly split into 3 main steps: Step 1: Shifting and Truncation. Given a polynomial system F := ( f , . . . , f n ) , a polydisc ∆ = ∆ r ( m ) , and an integer K ∈ { , . . . , d F } , we define a “precision” L := (cid:100) log r n ( K + 1) n (cid:101) . lgorithm 1: PolySol ( F , ∆ , K ) Input : Zero-dimensional system F = ( f , . . . , f n ) , polydisc ∆ = ∆ r ( m ) , and aninteger K ∈ { , . . . , d F } . Output: An integer k ∈ N ∪ {− } . If k ≥ , the polydisc ∆ contains exactly k solutions of F (counted with multiplicity). If k = − , nothing can be said. // * Shift and Truncation * // L := (cid:100) log r n ( K +1) n (cid:101) Compute a ( K + 1) · L -bit approximation Φ (cid:48) = ( φ (cid:48) , . . . , φ (cid:48) n ) of F [ m ] ≤ K = ( f [ m ] ≤ K , . . . , f n [ m ] ≤ K ) . // * Adding a degree ( K + 1) -perturbation ; as mentioned, this step seems to be onlynecessary in theory. In practice, we recommend to directly proceed with Φ := Φ (cid:48) . * // Φ( x ) := ( φ , . . . , φ n ) , with φ i := x K +1 i + φ (cid:48) i // * Solving the truncated system * // Compute a list ( ∆ , k ) , . . . , ( ∆ (cid:96) , k (cid:96) ) of disjoint polydiscs ∆ i = ∆ r i ( m i ) of radius atmost − L and corresponding multiplicities k i such that each ∆ i contains exactly k i solutions of Φ , and each solution of Φ is contained within one ∆ i . k := (cid:80) i : (cid:107) m i (cid:107) < r n k i k + := (cid:80) i : (cid:107) m i (cid:107) < nr k i if k = k + then // * Projection step * // Pick S ∈ S D Φ uniformly at random and compute the rotated system Φ ∗ := Φ ◦ S − = ( φ i ◦ S − ) i . for (cid:96) = 1 , . . . , n do ( b − (cid:96) , k − (cid:96) , LB − (cid:96) ) := T ∗ (∆ r √ n (0) , Res(Φ ∗ , x (cid:96) ))( b + (cid:96) , k + (cid:96) , LB + (cid:96) ) := T ∗ (∆ √ nr (0) , Res(Φ ∗ , x (cid:96) )) . if (cid:86) (cid:96) ∈ [ n ] b − (cid:96) ∧ (cid:86) (cid:96) ∈ [ n ] b + (cid:96) then // * Bound Computation and Comparison * // UB( m , r ) := r K +1 · d n F τ F +2 [ M ( m ) · ( n + d F ) ] d F B Φ ∗ := 2 · (cid:100) D Φ ∗ · (6 n + 4) · log( d Φ ∗ + n ) (cid:101) + τ Φ ∗ · D Φ ∗ k +1 LB( m , r ) := min (cid:96) min(LB − (cid:96) , LB + (cid:96) ) · (cid:16) n · (cid:0) D Φ ∗ + nD Φ ∗ (cid:1) · B Φ ∗ (cid:17) − if UB( m , r ) ≤ LB( m , r ) then return k return − Then, in a first step, we compute a ( K + 1) · L -bit approximation Φ (cid:48) ( x ) := ( φ (cid:48) , . . . , φ (cid:48) n ) , with φ (cid:48) i ∈ Q [ x ] ≤ K , of F [ m ] ≤ K ( x ) , i.e., we compute φ (cid:48) i such that (cid:107) φ (cid:48) i − f i [ m ] ≤ K (cid:107) < − ( K +1) L and ( K +1) L · φ (cid:48) i ∈ Z [ x ] i . Recall that the centered polynomial system F [ m ]( x ) was defined as F [ m ]( x ) := ( f [ m ]( x ) , . . . , f n [ m ]( x )) = ( f ( m + x ) , . . . , f n ( m + x )) , for m ∈ C n and that the truncation F [ m ] ≤ K ( x ) := ( f [ m ] ≤ K , . . . , f n [ m ] ≤ K ) , of the centered system F [ m ] , is defined by simply omitting all terms of f i [ m ]( x ) of total degreemore than K , as defined in Section 2.12.1. We further define Φ( x ) := ( φ , . . . , φ n ) , with φ i := x K +1 i + φ (cid:48) i , the system obtained by adding the term x K +1 i of degree K + 1 to the polynomial φ (cid:48) i . This stepseems to be odd at first sight, however, it ensures certain properties of Φ . In particular, Φ isguaranteed to have no zeros at infinity (and thus being zero-dimensional as well) according toCorollary 11 and our considerations in the corresponding example. This further implies that Φ = 0 has exactly D Φ = ( K + 1) n finite solutions counted with multiplicity. Also, our choiceof Φ allows us to bound the leading coefficient of Res(Φ , x (cid:96) ) for all (cid:96) = 1 , . . . , n , which turnsout to be useful in the analysis of our approach. Remark. In practice, the latter step does not seem to be necessary in most cases, and thuswe recommend to simply proceed with Φ := Φ (cid:48) and to check Φ (cid:48) for being zero-dimensional.Also, when implementing our algorithms, we observed that proceeding with Φ (cid:48) instead of Φ only improves the overall performance. Step 2: Solving Φ . We will later prove that, under the assumption that L is sufficientlylarge (or equivalently ∆ is sufficiently small), and z is a k -fold solution of the initial systemwith (cid:107) m − z (cid:107) < − L , the system Φ (as well as Φ (cid:48) for generic choice of its coefficients) yieldsa cluster of k (not necessarily distinct) solutions with norm less than · − L , whereas allother solutions have norm larger than δ (cid:29) − L . Here, δ is a constant that depends on thepolynomial system but not on L ; see Theorem 77 for the exact definition of δ and furtherdetails. We first check whether there exists a cluster of solutions of Φ near the origin that iswell separated from all other solutions of Φ . For this, we use a certified method (e.g. [BS16BS16])to compute all solutions of Φ . Here, by computing all solutions, it is meant to compute a setof disjoint discs, each of size less than − L , together with the number of solutions contained ineach disc such that the union of all discs contains all complex solutions. For the more involvedproblem of computing isolating regions of comparable size, the following theorem applies. Theorem 5. [BS16BS16, Thm. 9, 10] There is a Las Vegas algorithm to compute isolating regionsof size less than − ρ for all complex solutions of a zero-dimensional polynomial system F =( f , . . . , f n ) , with integer polynomials f i ∈ Z [ x ] , using ˜ O ( n ( n − ω +1)+1 ( nd F + τ F ) d F ( ω +2) n − ω − + n · d F n · ρ ) bit operations in expectation. Since ( K +1) L · φ i is a polynomial of degree K + 1 with integer coefficients of magnitude τ = O ( KL + n + τ F + d F log( m )) , we conclude from the above theorem that the cost forsolving the system Φ = 0 is bounded by ˜ O ( n ( n − ω +1)+1 ( nK + KL + τ + d F log( m )) · ( K + 1) ( ω +2) n − ω − ) (10)20it operations in expectation. Finally, we check whether the polydisc ∆ − := ∆ r/ (2 n ) (0) contains the same number k (cid:48) of solutions of Φ as the enlarged polydisc ∆ + := ∆ nr (0) .Notice that, from the above remark, this holds true if m is an L -bit approximation of a k -foldzero of F for large enough L as then · − L < r/ (2 n ) and δ > nr . If the zeros of Φ do notfulfill the latter condition, we return − . Otherwise, we proceed. Remark. We remark that computing the solutions of Φ is typically much more affordable thancomputing the solutions of the initial system F directly, in particular, in the case where n issmall and K (cid:28) d . Notice that, for n of constant size, the cost for solving the initial systemdirectly scales like d ( ω +2) n − ω − τ F , whereas the cost for solving the truncated system scaleslike ( KL + τ F + d log( m )) · ( K + 1) ( ω +2) n − ω − . Hence, for L and m of moderate size, therunning times might differ by factor of size ≈ ( d/ ( K + 1)) ( ω +2) n − ω − . Step 3: Passing from Φ to F . In the final step, we aim to certify that F [ m ] has the samenumber of zeros (i.e. k counted with multiplicity) in ∆ := ∆ r (0) as Φ . In order to do so, weaim to apply the following generalization of Rouché’s Theorem to Φ and F [ m ] , see [VH94VH94,Thm. 2.1] or [Llo75Llo75, Thm. 1] for a proof. Theorem 6 (Multidimensional Rouché) . Let F = ( f , . . . , f n ) and G = ( g , . . . , g n ) , with f i , g i ∈ C [ x ] for all i , define polynomial mappings from C n to C n . If, for a given boundeddomain D ⊂ C n , we have (cid:107)F ( x ) − G ( x ) (cid:107) < (cid:107)F ( x ) (cid:107) for all x ∈ ∂ D , where ∂ D is the boundary of D , then F and G have finitely many zeros in D and the numberof zeros (counted with multiplicities) of F and G in D is the same. In order to apply the above theorem to F := Φ and G := F [ m ] , we derive an upper bound UB( m , r ) on the absolute error sup x ∈ C n : (cid:107) x (cid:107) = r (cid:107)F [ m ]( x ) − Φ( x ) (cid:107) = sup x ∈ C n : (cid:107) x (cid:107) = r max i ∈ [ n ] | f i [ m ]( x ) − φ i ( x ) |≤ sup x ∈ C n : (cid:107) x (cid:107) = r max i ∈ [ n ] ( | f i [ m ]( x ) − φ (cid:48) i ( x ) | + | x i | K +1 ) ≤ r K +1 + sup x ∈ C n : (cid:107) x (cid:107) = r max i ∈ [ n ] | f i [ m ]( x ) − φ (cid:48) i ( x ) | when passing from Φ to F [ m ] as well as a lower bound LB( m , r ) on the norm of Φ( x ) on theboundary of the polydisc ∆ . The construction of UB( m , r ) is rather straightforward usingLemma 22. That is, we may choose UB( m , r ) := r K +1 · d F n · τ F +2 · [ M ( m ) · ( n + d F ) ] d F (11)In contrast, the construction of LB( m , r ) is more involved: We already mentioned that if L is large enough, then there are k zeros z , . . . , z k of Φ that have norm less than · − L ,whereas all other zeros have norm δ (cid:29) − L . Hence, under this assumption, picking arandom rotation matrix S from S D Φ = S ( K +1) n and considering a corresponding rotation of In practice, we recommend to consider S = id n and thus Φ ∗ = Φ as the initial choice as this turns out tobe sufficient in most cases. / , that theprojection of any zero of the “rotated system” Φ ∗ = ( φ ∗ , . . . , φ ∗ n ) := Φ ◦ S − on any coordinate axis, except for the k solutions S · z i , yields a value that is large comparedto r . Hence, in this case, the hidden-variable resultant R ∗ (cid:96) := Res(Φ ∗ , x (cid:96) ) of Φ ∗ has k rootsof absolute value less than r , whereas all other roots of R ∗ (cid:96) have absolute value (cid:29) r . Noticethat each φ ∗ j ∈ Q [ x ] is a polynomial of degree K + 1 with rational coefficients, and accordingto Lemma 11 and Lemma 33, we have (cid:107) φ ∗ j (cid:107) = 2 ˜ O ( n + τ F + d F log( m )) . Lemma 77 further yields the existence of an integer λ of absolute value ˜ O ( n log K ) with λ · S − ∈ Z n × n . Hence, we conclude that each term of degree K + 1 of λ K +1 · φ ∗ j has integer coefficients,and [CLO05CLO05, Theorem 3.1] further yields that Res( λ K +1 · Φ ∗ , x (cid:96) ) = λ n ( K +1) n · Res(Φ ∗ , x (cid:96) ) . It thus follows that | LC(Res(Φ ∗ , x (cid:96) )) | = λ − n ( K +1) n · | LC(Res( λ K +1 Φ ∗ , x (cid:96) )) | ≥ λ − n ( K +1) n = 2 − ˜ O (( K +1) n ) , where we use Corollary 11 to show that | LC(Res( λ K +1 Φ ∗ , x (cid:96) )) | is a positive integer, hence largerthan or equal to . Since Res(Φ ∗ , x (cid:96) ) is contained in the ideal generated by the polynomials φ ∗ j , we may write Res(Φ ∗ , x (cid:96) ) = g (cid:96), · φ ∗ + · · · + g (cid:96),n · φ ∗ n (12)with polynomials g (cid:96),j ∈ Q [ x ] of total degree bounded by D Φ ∗ = ( K + 1) n . Corollary 22 furtheryields the following upper bound on the size of the coefficients of the g (cid:96),j ’s: log (cid:107) g (cid:96),j (cid:107) ≤ B Φ ∗ = ˜ O ( D Φ ∗ · n + τ Φ ∗ · D Φ ∗ K + 1 )) = ˜ O (( K + 1) n + ( K + 1) n − · ( τ F + d log( m ))) . Using Lemma 11, part aa, this further yields a corresponding upper bound γ := (cid:18) n + D Φ ∗ D Φ ∗ (cid:19) · B Φ ∗ = 2 ˜ O (( K +1) n +( K +1) n − · ( τ F + d log( m ))) (13)such that max (cid:96),j sup x : (cid:107) x (cid:107)≤ | g (cid:96),j ( x ) | ≤ γ . Remark. The reader might wonder why we do not compute the above cofactor representation(1212) directly and then derive bounds on the size of max i,j sup x : (cid:107) x (cid:107) =1 | g (cid:96),j ( x ) | using intervalarithmetic, but instead use Corollary 22? The simple reason is that, at least in practice,computing the polynomials g i,j turns out to be considerably more costly than computing theresultant polynomials R ∗ (cid:96) ( x ) = Res(Φ ∗ , x (cid:96) ) only. In contrast, our approach of computingthe bound γ does not require to compute the polynomials g i,j , and thus comes at almost noadditional cost. We further remark at this point that we will use the bounds from Corollary 22in our complexity analysis of the algorithm.In the next step, we compute lower bounds LB − (cid:96) and LB + (cid:96) for | R (cid:96) ( x ) ∗ | on the boundaryof the two discs D − := ∆ r/ √ n (0) ⊂ C and D + := ∆ r ·√ n (0) ⊂ C , respectively. For this, weuse the so-called T k -test, an approach that has recently been proposed in an algorithm forcomplex root isolation [BSS+15BSS+15]. 22 emma 9 ([BSS+15BSS+15]) . Let f ∈ C [ x ] be a uni-variate polynomial of degree d and let ∆ :=∆ r (0) ⊂ C be the disc with radius r centered at . The so-called T k -test returns a pair T k (∆ , f ) = ( b, LB) = (cid:16) | f ( k ) (0) | r k k ! − · (cid:88) i (cid:54) = k | f ( i ) (0) | r i i ! > , · | f ( k ) (0) | r k k ! (cid:17) . (14) If b = True , we say that T k (∆ , f ) succeeds. If T k (∆ , f ) succeeds, ∆ contains exactly k rootscounted with multiplicity and · LB > | f ( x ) | > LB for all x ∈ ∂ ∆ r (0) . In addition, if ∆ r/ (16 d ) (0) as well as ∆ d r (0) contain exactly k roots, then T k (∆ , f ) succeeds.We further define T (cid:63) (∆ , f ) = (cid:40)(cid:16) True , k, · | f ( k ) (0) | r k k ! (cid:17) if T k (∆ , f ) succeeds for some k, ( False , − otherwise. (15)Now, suppose that T ∗ ( D − , R ∗ (cid:96) ) = ( b − (cid:96) , k − (cid:96) , LB − (cid:96) ) as well as T ∗ ( D + , R ∗ (cid:96) ) = ( b + (cid:96) , k + (cid:96) , LB + (cid:96) ) succeed for all (cid:96) , then min(LB − (cid:96) , LB + (cid:96) ) constitutes a lower bound for | R ∗ (cid:96) | on the boundary of D − as well as D + . From (1212), (1313), and the definition of LB( m , r ) , we now conclude that (cid:107) Φ ∗ ( x ) (cid:107) > LB( m , r ) := min (cid:96) min(LB − (cid:96) , LB + (cid:96) ) nγ , for all x with (cid:107) x (cid:107) = r √ n or (cid:107) x (cid:107) = √ n · r. (16)Since the maximum and minimum of a holomorphic function (in several variables) on abounded domain is taken at its boundary, we further conclude that the above inequality holdsfor any x with r/ √ n ≤ (cid:107) x (cid:107) ≤ √ n · r . Notice that the rotation of the system by means of therotation matrix maintains the -norm (cid:107) . (cid:107) of any point. Thus, the the norm of any point x differs from the norm of the rotated point S · x by a factor that is lower and upper boundedby r/ √ n and √ n · r , respectively. Hence, from the above bound on (cid:107) Φ ∗ (cid:107) , we conclude that (cid:107) Φ( x ) (cid:107) > LB( m , r ) for any x with (cid:107) x (cid:107) = r. (17)Now in order to apply Rouché’s Theorem to Φ and F [ m ] , it suffices to check whether LB( m , r ) > UB( m , r ) , in which case we have shown that Φ and F [ m ] have the same numberof roots in ∆ r (0) . Hence, we return True in this case. Otherwise, the algorithm returns False .In the next section, we will show that, if m is a sufficiently good approximation (i.e. for largeenough L ) of a k -fold solution of F , our algorithm succeeds. Here, we only give an informalargument: Notice that, for large L , the bound UB( m , r ) scales like C · r K +1 for some constant C . The bound γ does not depend on L , hence LB( m , r ) scales like min (cid:96) min(LB − (cid:96) , LB + (cid:96) ) for large enough L . However, in this situation, each R ∗ (cid:96) has a cluster of k roots near theorigin that is well separated from all of its remaining roots, and thus min(LB − (cid:96) , LB + (cid:96) ) scaleslike [ | R ∗ (cid:96) ( k ) (0) | / ( √ n k k !)] · r k . Hence, we conclude that LB( m , r ) scales like C (cid:48) · r k for someconstant C (cid:48) , which implies that LB( m , r ) must be smaller than UB( m , r ) for large enough L .We remark that the precise argument is slightly more involved as many subtleties need to beaddressed. In particular, we need to show that | R ∗ (cid:96) ( k ) (0) | / ( √ n k k !) does not depend on r if r is small enough, even though the definition of R ∗ strongly depends on the choice of m , r , andthe rotation matrix S . We will give details in the next section.23 Analysis We start by introducing some further notation. For a zero-dimensional polynomial system F = ( f , . . . , f n ) in n variables, let z , . . . , z N denote its zeros. We define σ ( z i , F ) := min j (cid:54) = i (cid:107) z i − z j (cid:107) and ∂ ( z i , F ) := (cid:89) j (cid:54) = i (cid:107) z i − z j (cid:107) µ ( z j , F ) to be the separation of z i with respect to F and the geometric derivative of F at z i , respectively.We remark that these terms are derived from the interpretation of these quantities in theunivariate case, where the separation of a root z of a polynomial f ∈ C [ x ] is defined inexactly the same way, and the first non-vanishing derivative (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∂ µ ( z ,f ) f∂z µ ( z ,f ) ( z ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = LC( f ) · (cid:89) z (cid:54) = z : f ( z )=0 | z − z | µ ( z,f ) of f at z can be expressed as a product involving the leading coefficient of f and the distancesbetween z and the other roots. We first provide some bounds on (cid:107) z i (cid:107) , σ ( z i , F ) , and ∂ ( z i , F ) for the special case where each f i has only integer coefficients. For similar bounds that arealso adaptive with respect to the sparseness of the given system, we refer to [EMT10EMT10]. Lemma 10. Let F = ( f i ) i =1 ,...,n be a zero-dimensional system with integer polynomials f i ,and let z , . . . , z N denote the zeros of F . Then it holds: N (cid:88) i =1 µ ( z i , F ) · log( z i ) = ˜ O ( n · B F ) = ˜ O ( n · [ n · D F + τ F · max i D F d i ]) , with B F as in (77) | log σ ( z i , F ) | = ˜ O ( D F · B F ) , log( ∂ ( z i , F ) − ) = ˜ O ( n · D F · [ B F + log( z i )]) = ˜ O ( n · D F · B F ) , and log( σ ( z i , F ) − ) = ˜ O ( D F · log( z i ) + n · B F + log( ∂ ( z i , F ) − )) = ˜ O ( n · D F · B F ) . Proof. From Corollary 22, we conclude that Res( F , x (cid:96) ) is an integer polynomial of magnitude ( D F , B F ) for all (cid:96) = 1 , . . . , n . Since the (cid:96) -th coordinate z i,(cid:96) of each solution of F = 0 is a rootof multiplicity at least µ ( z i , F ) of Res( F , x (cid:96) ) and since the Mahler measure Mea(Res( F , x (cid:96) )) = LC(Res( F , x (cid:96) )) · (cid:89) z ∈ C :Res( F ,x (cid:96) )( z )=0 M ( z ) µ ( z, Res( F ,x (cid:96) )) of Res( F , x (cid:96) ) is upper bounded by its -norm (cid:107) Res( F , x (cid:96) )) (cid:107) ≤ √ D F · B F (e.g. see [Yap00Yap00]),it follows that N (cid:88) i =1 µ ( z i , F ) · log( z i ) ≤ D F + n (cid:88) (cid:96) =1 log(Mea(Res( F , x (cid:96) ))) = ˜ O ( nB F ) . For the second claim, notice that σ ( z i , F ) ≥ σ ( z i,(cid:96) , Res( F , x (cid:96) )) for at least one (cid:96) (as twodistinct solutions must differ in at least one coordinate), and that the separation of an integerpolynomial of magnitude ( D F , B F ) is lower bounded by − ˜ O ( D F · B F ) ; e.g. see [MSW15MSW15] for aproof. For the bound on ∂ ( z i , F ) , notice that ∂ ( z i , F ) ≥ (cid:81) n(cid:96) =1 ∂ ( z i,(cid:96) , Res( F , x (cid:96) )) (cid:81) n(cid:96) =1 (cid:81) z (cid:54) = z i,(cid:96) :Res( F ,x (cid:96) )( z )=0 M ( z − z i,(cid:96) ) µ ( z, Res( F ,x (cid:96) )) . ∂ ( z , f ) = 2 − ˜ O ( dL ) for any root z ofa polynomial f ∈ Z [ x ] of magnitude ( d, L ) . This shows that ∂ ( z i,(cid:96) , Res( F , x (cid:96) )) = 2 − ˜ O ( D F B F ) for all (cid:96) . It remains to derive an upper bound on the denominator in the above fraction. Forthis, we define R (cid:96) := Res( F , x (cid:96) )[ z i,(cid:96) ] . Then, it holds that (cid:89) n(cid:96) =1 (cid:89) z (cid:54) = z i,(cid:96) :Res( F ,x (cid:96) )( z )=0 max(1 , | z − z i,(cid:96) | ) µ ( z, Res( F ,x (cid:96) )) = n (cid:89) (cid:96) =1 Mea( R (cid:96) )LC( R (cid:96) ) . According to Lemma 11, R (cid:96) is a polynomial of magnitude ( D F , ˜ O ( B F + D F · log( z i,(cid:96) ))) , and, inaddition, it has the same leading coefficient as Res( F , x (cid:96) ) . In particular, its leading coefficientis a non-zero integer, and thus of absolute value larger than or equal to . Thus, we have Mea( R (cid:96) )LC( R (cid:96) ) ≤ (cid:107) R (cid:96) (cid:107) = 2 ˜ O ( B F + D F · log( z i,(cid:96) )) = 2 ˜ O ( B F + D F · log( z i )) , which shows that log( ∂ ( z i , F ) − ) = ˜ O ( nD F ( B F + log( z i ))) .For the last claim, notice that σ ( z i , F ) appears as one of the factors in the definition of ∂ ( z i , F ) . Since the product of all remaining factors is upper bounded by (cid:89) z j (cid:54) = z i M ( z j − z i ) µ ( z j , F ) ≤ (cid:89) z j (cid:54) = z i [2 · M ( z j ) · M ( z i )] µ ( z j , F ) ≤ D F · M ( z i ) D F · N (cid:89) j =1 M ( z j ) µ ( z j , F ) , the claim follows directly from the bound on (cid:80) Ni =1 µ ( z i , F ) · log( z i ) and on log( ∂ ( z i , F ) − ) .We are now ready to derive one of our main results in this paper. More specifically,the following theorem shows that, in a sufficiently small neighborhood (which we will alsoquantify) of a k -fold solution z = z i of F = 0 , (cid:107)F ( x ) (cid:107) scales like c · (cid:107) x (cid:107) k with c a constant.We further argue that this implies that a sufficiently good approximation Φ of the shifted andtruncated system F [ z ] ≤ K , with arbitrary K ≥ k , has a cluster of k solutions near the origin,whereas all remaining solutions are well separated from this cluster. We also give bounds onthe approximation error that involve the quantities σ ( z , F ) and ∂ ( z , F ) that are intrinsic tothe hardness of the given polynomial system. Theorem 7. Let F be a zero-dimensional system, z a zero of F of multiplicity k , and m be anapproximation of z with (cid:107) m − z (cid:107) < − L . Let K ≥ k , and Φ (cid:48) = ( φ (cid:48) i ) i =1 ,...,n be a ( K + 1) · L -bitapproximation of F [ m ] ≤ K with polynomials φ (cid:48) i of degree at most K , and let a , . . . , a n ∈ C bearbitrary complex values of magnitude ≤ | a i | ≤ for all i . Then, the polynomial system Φ := ( φ (cid:48) i + a i · x K +1 i ) i =1 ,...,n is zero-dimensional, and there exists an L ∈ N such that, for any L ≥ L , Φ has exactly k zeros (counted with multiplicity) of norm smaller than · − L , whereas all other zeros havenorm larger than δ := σ ( z , F )(2 n D F ) n . In the special case, where each polynomial in F has onlyinteger coefficients, it holds that L = ˜ O ( n · D F · [ n + max i τ F + d F d i + log( z )] + log( ∂ ( z , F ) − ))) . roof. We denote z , . . . , z N , with z = z i , the zeros of F . Let S ∈ S D F be an admissiblerotation matrix with respect to F as well as with respect to the shifted system F [ z ] . Noticethat such a matrix exists as more than half of the matrices in S D F are admissible withrespect to F and more than half of the matrices are admissible with respect to F [ z ] . Let F ∗ := F ◦ S − be the corresponding “rotation” of F and z ∗ , . . . , z ∗ N be the zeros of F ∗ suchthat z ∗ j = ( z ∗ j, , . . . , z ∗ j,n ) = S · z j . Since S is admissible with respect to F [ z ] , Lemma 88 yieldsthat | z ∗ j,(cid:96) − z ∗ i,(cid:96) | ≥ (2 n D F ) − n · (cid:107) z j − z (cid:107) ≥ σ ( z , F )(2 n D F ) n for all (cid:96) and j (cid:54) = i . In addition, since S is also admissible with respect to F , Lemma 44 andLemma 88 guarantees that each root of the resultant polynomial Res( F ∗ , x (cid:96) ) is the projectionof a finite zero of F ∗ on the x (cid:96) -coordinate. Thus, Res( F ∗ , x (cid:96) ) has a k -fold root at z ∗ i,(cid:96) , whereasall other roots z ∗ j,(cid:96) of Res( F ∗ , x (cid:96) ) have distance at least σ ( z , F )(2 n D F ) n to z ∗ i,(cid:96) . Now, applyingLemma 99 to a disc with center z ∗ i,(cid:96) and arbitrary radius smaller than r ∗ := σ ( z , F )16 D F · (2 n D F ) n , yields that | Res( F ∗ , x (cid:96) )( x ) | > k ! · (cid:12)(cid:12)(cid:12)(cid:12) ∂ k Res( F ∗ , x (cid:96) ) ∂x k ( z ∗ i,(cid:96) ) (cid:12)(cid:12)(cid:12)(cid:12) · | x − z ∗ i,(cid:96) | k for all x ∈ C with | x − z ∗ i,(cid:96) | < r ∗ . Denoting LC (cid:96) := LC(Res( F ∗ , x (cid:96) )) , this further yields (cid:12)(cid:12)(cid:12)(cid:12) ∂ k Res( F ∗ , x (cid:96) ) ∂x k ( z ∗ i,(cid:96) ) (cid:12)(cid:12)(cid:12)(cid:12) = | LC (cid:96) | · (cid:89) j (cid:54) = i | z ∗ i,(cid:96) − z ∗ j,(cid:96) | µ ( z ∗ j,(cid:96) , Res( F ∗ ,x (cid:96) )) ≥ | LC (cid:96) | · (cid:89) j (cid:54) = i (cid:18) (cid:107) z − z j (cid:107) (2 n D F ) n (cid:19) µ ( z j , F ) ≥ | LC (cid:96) | (2 n D F ) nD F · ∂ ( z , F ) , and thus it follows that | Res( F ∗ , x (cid:96) )( x ) | > | LC (cid:96) | k ! · (2 n D F ) nD F · ∂ ( z , F ) · | x − z ∗ i,(cid:96) | k > | LC (cid:96) | D F ! · (2 n D F ) nD F · ∂ ( z , F ) · | x − z ∗ i,(cid:96) | k > | LC (cid:96) | (4 n D F ) nD F · ∂ ( z , F ) · | x − z ∗ i,(cid:96) | k if | x − z ∗ i,(cid:96) | < r ∗ . (18)Furthermore, Res( F ∗ , x (cid:96) ) is contained in the ideal spanned by the polynomials F ∗ =( f ∗ , . . . , f ∗ n ) , that is, there exist polynomials g (cid:96),j ∈ C [ x ] with Res( F ∗ , x (cid:96) ) = (cid:80) nj =1 g (cid:96),j f ∗ j .According to Corollary 22, we may assume that log (cid:107) g (cid:96),j (cid:107) ≤ B F ∗ = ˜ O ( n · D F ∗ + τ F ∗ · max i D F d i ) = ˜ O ( n · D F + ( τ F + d F ) · max i D F d i ) (cid:96), j , where the last inequality follows from Lemma 33. Using Lemma 11 then implies that γ F ∗ := (cid:18) n + D F D F (cid:19) · B F∗ · ( M ( z ) + 1) D F ≥ sup x : (cid:107) x − z ∗ i (cid:107)≤ | g i,j ( x ) | for all (cid:96), j. (19)Now, combining (1818) and (1919) yields (cid:107)F ∗ ( x ) (cid:107) ≥ (cid:104) min (cid:96) | LC (cid:96) | n · γ F ∗ · (4 n D F ) nD F · ∂ ( z , F ) (cid:105) · (cid:107) x − z ∗ i (cid:107) k for all x ∈ C n with (cid:107) x − z ∗ i (cid:107) < r ∗ .So what can we conclude about our initial (non-rotated) system? Since a rotation main-tains the Euclidean distance and since the max-norm differs from the Euclidean norm by afactor of at most √ n , it follows that a point x of max-norm (cid:107) x (cid:107) is rotated via S (or S − )onto a point x (cid:48) of max-norm (cid:107) x (cid:48) (cid:107) ≤ (cid:107) x (cid:48) (cid:107) = (cid:107) x (cid:107) ≤ √ n · (cid:107) x (cid:107) . Hence, it holds that everypoint x ∗ = S x with (cid:107) x − z (cid:107) < r := r ∗ / √ n satisfies (cid:107) x ∗ − z ∗ i (cid:107) < r ∗ . Thus, for all x with (cid:107) x − z (cid:107) < r , it holds that (cid:107)F ( x ) (cid:107) ≥ (cid:104) min (cid:96) | LC (cid:96) | n · √ n k · γ F ∗ · (4 n D F ) nD F · ∂ ( z , F ) (cid:105)(cid:124) (cid:123)(cid:122) (cid:125) =: c ·(cid:107) x − z (cid:107) k (20)Now, suppose that L ≥ log(8 /r ) and thus − L < r . Since m is an approximation of z with (cid:107) z − m (cid:107) < − L , F [ m ] has exactly one solution (namely, ˆz := z − m ) of multiplicity k in ∆ r / (0) and (cid:107)F [ m ]( x ) (cid:107) = (cid:107)F ( x + m ) (cid:107) > c k · (cid:107) x (cid:107) k for all x ∈ C n with − L +1 < (cid:107) x (cid:107) < r , as, for such x , it holds that (cid:107) x (cid:107) / < (cid:107) x + m − z (cid:107) < r . Applying Lemma 22 to each φ i andusing the fact that M ( m ) ≤ M ( z ) then shows that, for all x with − L +1 < (cid:107) x (cid:107) < r / , itholds that (cid:107) Φ( x ) − F [ m ]( x ) (cid:107) ≤ − ( K +1) L + (cid:107) x (cid:107) K +1 · d F n · τ F + d F · [ M ( z ) · ( n + d F ) ] d F ≤ (cid:107) x (cid:107) K +1 · d n F · τ F + d F +1 · [ M ( z ) · ( n + d F ) ] d F . Notice that, due to the construction of Φ and Corollary 11, Φ is zero-dimensional. Hence,Rouché’s Theorem applied to F [ m ] and Φ shows that the polydisc ∆ ρ (0) contains the samenumber of solutions of Φ and F [ m ] if − L +1 < ρ < r / and if, in addition, ρ fulfills thefollowing inequality ρ K +1 · d n F · τ F + d F +1 · [ M ( z ) · ( n + d F ) ] d F < c k · ρ k . Equivalently, we must have − L +1 < ρ < r / and ρ K +1 − k < c k · d n F · τ F + d F +1 [ M ( z ) · ( n + d F ) ] d F . L > L := max (cid:20) log 8 r , log 2 k · d F n τ F + d F +1 [ M ( z ) · ( n + d F ) ] d F c (cid:21) = ˜ O ( D F · ( n + max i τ F + d F d i + log( z )) + log( ∂ ( z , F ) − ) + log( σ ( z , F ) − ) − | log min (cid:96) | LC (cid:96) || ) . (21)each polydisc ∆ ρ (cid:48) ( z ) , with arbitrary radius ρ (cid:48) ∈ (4 · − L , r / , contains exactly k zeros of Φ .Since δ < r / , this proves the first part of the theorem.It remains to prove the claim bound on L for the special case, where F is a polynomialsystem defined over the integers. For this, we need to estimate the size of the leading coefficientof Res( F ∗ , x (cid:96) ) . Notice that there exists an integer λ of size ˜ O ( n log d F ) with λ · S − ∈ Z n × n ,and thus F (cid:48) : ( f (cid:48) , . . . , f (cid:48) n ) = ( λ deg f ∗ · f ∗ , . . . , λ deg f ∗ n · f ∗ n ) is a polynomial system with integer coefficients, which shows that | LC(Res( F (cid:48) , x (cid:96) )) | ≥ .Using [CLO05CLO05, Thm. 2.3 and 3.5] then shows that | LC (cid:96) | = λ − nD F · | LC(Res( F (cid:48) , x (cid:96) )) | ≥ λ − nD F = 2 − O ( n D F ) . Hence, the bound follows from (2121) and the bound for log σ ( z , F ) − from Lemma 1010.From the previous Theorem, we now immediately obtain the following result by setting m := z and Φ := F [ z ] ≤ K for an arbitrary K ≥ k . Corollary 3. Let z be a k -fold zero of a zero dimensional system F and K ≥ k . Then, F [ z ] ≤ K has a k -fold zero at the origin, and all other zeros have norm larger than δ := σ ( z ,F )(2 n D F ) n . We can now show that Algorithm 11 terminates and yields a correct result assuming that L is large enough and the oracle, which provides an approximation m of the solution z , returnsa correct answer. Theorem 8. If PolySol ( F , ∆ , K ) returns an integer k ≥ , then the polydisc ∆ = ∆ r ( m ) contains exactly k solutions of F counted with multiplicity. Vice versa, suppose that z is a solution of F of multiplicity k and K ≥ k , then there exists a positive integer L ∗ of size L ∗ = ˜ O ( L + ( K + 1) n · log 1 δ + d F · log( z )) . with L and δ as in Theorem 77, such that PolySol ( F , ∆ , K ) returns k with probabilityat least / if r ≤ − L ∗ and (cid:107) z − m (cid:107) < r n ( K +1) n . If F has only integer coefficients, it holds: L ∗ = ˜ O ( D F · max i d F + τ F d i + D F · log( z ) + log( ∂ ( z , F ) − ) + ( K + 1) n · log( σ ( z , F ) − ))= ˜ O (( K + 1) n · D F · [ D F + τ F · max (cid:96) D F d (cid:96) ]) . We remark that F [ z ] ≤ K does not necessarily have to be zero-dimensional. Rouché’s Theorem only guar-antees that the polydisc ∆ δ (0) contains exactly k zeros of F [ z ] ≤ K . roof. For the first part, we proceed similarly as in the proof of Theorem 77, however, wework with the system Φ instead of the initial system F . From Line 77 in the algorithm, wealready know that Φ has exactly k solutions with ( (cid:107) . (cid:107) -) norm less than rn , whereas all othersolutions have norm at least nr . Now, when considering a random rotation matrix S ∈ S D Φ ,the corresponding rotated system Φ ∗ = Φ ◦ S − has exactly k solutions with norm less than r √ n , whereas all other solutions have norm at least √ n · r . We may now further write Res(Φ ∗ , x (cid:96) ) = γ (cid:96), · φ ∗ + · · · + γ (cid:96),n · φ ∗ n (22)with polynomials γ (cid:96),j ∈ C [ x ] . From Part dd of Lemma 11 and Corollary 22 we conclude that (cid:18) D Φ ∗ + nD Φ ∗ (cid:19) · B Φ ∗ ≥ sup x : (cid:107) x (cid:107)≤ | γ (cid:96),j ( x ) | for all (cid:96), j. In addition, since T ∗ (∆ r √ n (0) , Res(Φ ∗ , x (cid:96) )) = ( True , k − (cid:96) , LB − (cid:96) ) and T ∗ (∆ √ nr (0) , Res(Φ ∗ , x (cid:96) )) =( True , k + (cid:96) , LB + (cid:96) ) for all (cid:96) (Line 1010), we conclude from Lemma 99 that | Res(Φ ∗ , x (cid:96) )( x ) | > min(LB − (cid:96) , LB + (cid:96) ) for all i and all x with | x (cid:96) | = r √ n or | x (cid:96) | = √ nr. Hence, using (2222), this shows that (cid:107) Φ ∗ ( x ) (cid:107) ≥ min (cid:96) =1 ,...,n min(LB − (cid:96) , LB + (cid:96) ) n · (cid:0) D Φ ∗ + nD Φ ∗ (cid:1) · B Φ ∗ for all x with (cid:107) x (cid:107) = r √ n or (cid:107) x (cid:107) = √ nr. Since a holomorphic mapping cannot take its minimum or maximum in the interior of somedomain, it thus follows that the above inequality even holds for any x with r √ n ≤ (cid:107) x (cid:107) ≤ √ nr .It thus follows that (cid:107) Φ( x ) (cid:107) ≥ min (cid:96) =1 ,...,n min(LB − (cid:96) , LB + (cid:96) ) n · (cid:0) D Φ ∗ + nD Φ ∗ (cid:1) · B Φ ∗ for all x with (cid:107) x (cid:107) = r. According to (1111), UB( m , r ) constitutes an upper bound on the error (cid:107)F [ m ]( x ) − Φ( x ) (cid:107) forany x with (cid:107) x (cid:107) ≤ . Hence, in particular, we also have (cid:107)F [ m ]( x ) − Φ( x ) (cid:107) ≤ UB( m , r ) for all x with x with (cid:107) x (cid:107) = r. Hence, using Rouché’s Theorem, we conclude that F [ m ] and Φ have the same number ofsolutions in the polydisc ∆ r (0) .It remains to prove the second claim. For this, suppose that F has a k -fold solution at z with (cid:107) m − z (cid:107) < r n ( K +1) n and that L := (cid:100) log 32 n ( K + 1) n r (cid:101) > L := max( L , log 1 δ + log[2048 · n · (2 n D F ( K + 1) n ) n ])= ˜ O ( L + log 1 δ + n log D F ) , (23)with δ = σ ( z , F )(2 n D F ) n as defined in Theorem 77. Let Φ be an approximation of F [ m ] ≤ K asdefined in Theorem 77. Then, Φ has k solutions z , . . . , z k of norm (cid:107) z i (cid:107) < · − L < r/ (2 n ) ,29hereas all remaining solutions, denoted by ¯ z , . . . , ¯ z m , have norm (cid:107) ¯ z j (cid:107) > δ > nr . We thusconclude that the if-condition is satisfied in Line 77. Now, when choosing a random rotationmatrix S ∈ S D Φ , the solutions z i near the origin are mapped to solutions z ∗ i := S ◦ z i of Φ ∗ = Φ ◦ S − with norm (cid:107) z ∗ i (cid:107) < √ n · − L ≤ r √ n ( K +1) n , whereas the remaining solutions aremapped to solutions ¯ z ∗ j of Φ ∗ with norm (cid:107) ¯ z ∗ j (cid:107) > δ / √ n > √ nr . In addition, with probabilitymore than / , we have | z ∗ j,(cid:96) | ≥ (cid:107) ¯ z j (cid:107) · (2 n D Φ ∗ ) − n > δ · (2 n ( K + 1) n ) − n > · √ n · ( K + 1) n r for all j ∈ [ m ] and all (cid:96) ∈ [ n ] . This implies that each of the resultant polynomials Res(Φ ∗ , x (cid:96) ) has k roots of absolute value less than r K +1) n √ n , whereas all remaining roots are of absolutevalue larger than · √ n · ( K + 1) n r . Notice that each polynomial Res(Φ ∗ , x (cid:96) ) has degree ( K + 1) n , and thus Lemma 99 guarantees success in Line 1010 of the algorithm. Now, recallthe lower bounds LB − (cid:96) and LB + (cid:96) for | Res(Φ ∗ , x (cid:96) ) | on the boundary of ∆ r/ √ n (0) and ∆ √ nr (0) ,respectively, as computed in Line 1111. For arbitrary x ∈ C with | x | = r/ √ n , we have LB − (cid:96) ≥ | Res(Φ ∗ , x (cid:96) )( x ) | > | LC(Res(Φ ∗ , x (cid:96) )) | · (cid:18) r √ n (cid:19) k · (cid:18) δ · (2 n ( K + 1) n ) − n (cid:19) ( K +1) n − k > | LC(Res(Φ ∗ , x (cid:96) )) | · (2 n ( K + 1) n ) − n ( K +1) n · r k · δ ( K +1) n . Using the fact LB + (cid:96) ≥ | Res(Φ ∗ ,x (cid:96) )( x ) | for all x with | x | = 2 √ nr , an analogous computationshows that LB + (cid:96) fulfills the same bound, that is, LB + (cid:96) ≥ | LC(Res(Φ ∗ , x (cid:96) )) | · (2 n ( K + 1) n ) − n ( K +1) n · r k · δ ( K +1) n . From Lemma 11 and our construction of Φ ∗ , the leading coefficient of each polynomial Res(Φ ∗ , x (cid:96) ) is a non-zero integer, hence we obtain that LB( m , r ) = min (cid:96) =1 ,...,n min(LB − (cid:96) , LB + (cid:96) ) n · (cid:0) D Φ ∗ + nD Φ ∗ (cid:1) · B Φ ∗ ≥ δ ( K +1) n n (2 n ( K + 1) n ) n ( K +1) n · (cid:0) D Φ ∗ + nD Φ ∗ (cid:1) · B Φ ∗ (cid:124) (cid:123)(cid:122) (cid:125) =: C · r k . Notice that, for small r , LB( m , r ) scales like C · r k , with a constant C that does not dependon r . The upper bound UB( m , r ) = r K +1 · d F n τ F +2 [ M ( m ) · ( n + d F ) ] d F (cid:124) (cid:123)(cid:122) (cid:125) =: C (cid:48) scales like C (cid:48) · r K +1 , and thus our algorithm succeeds if r fulfills the condition in (2323) (i.e. (cid:100) log n ( K +1) n r (cid:101) ≥ L ) and r K − k +1 < CC (cid:48) . Both condition are fulfilled if log(1 /r ) ≥ L ∗ := max( L , log C (cid:48) C ) = ˜ O ( L + ( K + 1) n · log 1 δ + d F · log( z )) . The claimed bound on L ∗ for the special case where F is defined over the integers followsdirectly from the corresponding bound on L from Theorem 77 and our bounds on log( z ) , log( σ ( z , F ) − ) , and log( ∂ ( z , F ) − ) from Lemma 1010.30 Application: Computing the Zeros of a Bivariate System In this section, we report on an application of our technique in the context of eliminationmethods for the bivariate case. More precisely, we incorporate the algorithm PolySol as an inclusion predicate in the Bisolve algorithm [BEK+13BEK+13; KS15KS15]. Comparing Sage implementations of the original Bisolve algorithm and its modified variant, we empiricallyshow that the idea of truncating the original system with respect to the multiplicity of thesolution yields a considerable performance improvement. Bisolve is a classical eliminationmethod for computing the real [BEK+13BEK+13] or complex [KS15KS15] zeros within a given polydisc ∆ = ∆ × ∆ ⊂ C of a bivariate system F : f ( x , x ) = f ( x , x ) = 0 , with polynomials f , f ∈ Z [ x , x ] . It achieves the best known complexity bound (i.e. ˜ O ( d F + d F · τ F ) bit operations for computingall complex solution) that is currently known for this problem, and its implementation showssuperior performance when compared to other complete and certified methods. As we aim tomodify the Bisolve algorithm at some crucial steps, we start with a brief description of theoriginal version. Bisolve in a Nutshell. In an initial projection phase , Bisolve computes a set C of candidate regions using resultant computation and univariate root finding. More specifically,we first compute the hidden-variable resultants R (cid:96) ( x ) := Res( F , x (cid:96) ) for (cid:96) = 1 , . Then, for eachroot z (cid:96),i in ∆ (cid:96) , we compute an isolating disc ∆ (cid:96),i such that T (cid:63) (∆ (cid:96),i , R (cid:96) ) = ( True , k (cid:96),i , LB (cid:96),i ) .That is, the T (cid:63) -test succeeds and yields the multiplicity of z (cid:96),i as a root of R (cid:96) as well as a lowerbound for | R (cid:96) | on the boundary of ∆ (cid:96),i . By taking the pairwise product of any two discs ∆ ,i and ∆ ,j , we obtain a set C of polydiscs ∆ i,j := ∆ ,i × ∆ ,j in C . Notice that each solution z in ∆ of F must be one of the candidate solutions z i,j := ( z ,i , z ,j ) as each coordinate of z is a root of the corresponding polynomial R (cid:96) . Hence, each solutions must be contained in oneof the candidate regions, even though most candidate regions do not contain any solution. Inaddition, each candidate region ∆ i,j contains at most one solution, which must be z i,j In the validation phase , the algorithm checks for every candidate region ∆ i,j whether itcontains a solution or not. In other words, we check whether the corresponding candidatesolution z = z i,j is actually a solution or not. The approach used in Bisolve shares manysimilarities to the algorithm PolySol as proposed in this paper. That is, we write R (cid:96) = g (cid:96), · f + g (cid:96), · f with g (cid:96), , g ,(cid:96) ∈ Z [ x, y ] and (cid:96) = 1 , and compute an upper bound UB for | g (cid:96),i ( x ) | for (cid:96) = 1 , , i = 1 , , and arbitrary x ∈ ∆ i,j .Similar as in PolySol , this is achieved without actually computing the polynomials g (cid:96), and g (cid:96), , but by exploiting the fact that these polynomials can be written as determinantsof “Sylvester-like” matrices ; see [KS15KS15] for details. Together with the lower bounds LB ,i and LB ,j as computed above this yields a lower bound LB ∗ = min(LB ,i , LB ,j )2 UB for (cid:107)F (cid:107) ∞ =max( | f | , | f | ) on the boundary of ∆ i,j . Notice that this is one crucial point, where our novel approach differs from Bisolve . Namely, for PolySol , we use the results from [DKS13DKS13] on the arithmetic Nullstellensatz to derive corresponding boundson the cofactors g (cid:96),j . This was necessary for generalizing the method to arbitrary dimension. Another crucialdifference is that no truncation of the system is considered in Bisolve . z as a solution, Bisolve proceed in rounds, wherea m -bit approximation ζ of z is computed at the beginning of the m -th round. As an exclusion predicate , interval arithmetic is used in order to compute a superset (cid:3) f (cid:96) ( ∆ − m ( ζ )) of f (cid:96) ( ∆ − m ( ζ )) for (cid:96) = 1 , . If we can show that either f or f does not vanish, the candidateis discarded. As an inclusion predicate the above lower bound LB ∗ on the boundary of ∆ i,j is compared to the values that f and f take at the approximation ζ of the candidate z .More specifically, if max( | f ( ζ ) | , | g ( ζ ) | ) < LB ∗ , then ∆ i,j contains a solution; see Theorem 4in [BEK+13BEK+13]. If neither the exclusion nor the inclusion predicate applies, we proceed with thenext round. The BisolvePlus routine. Notice that, even though Bisolve computes the set Z := { z i,j ∈ ∆ : f ( z i,j ) = f ( z i,j ) = 0 } of all solutions of F within ∆ , it does not reveal the multiplicity k of a specific solution z = z i,j = ( z ,i , z ,j ) ∈ Z . However, due to the properties of the resultant polynomials, itholds that k = µ ( z ,i , R ) if the following two conditions are both fulfilled: deg x f = deg f and deg x f = deg f (24) ∀ x ∈ C \ { z ,j } : f ( z ,i , x ) (cid:54) = 0 or f ( z ,i , x ) (cid:54) = 0 (25)The first condition guarantees that there is no solution of F at infinity above any z ∈ C ,whereas the second condition guarantees that there is no other finite (complex) solution of F that shares the first coordinate with z . We remark that it is easy to check the first condition,however, checking the second condition is more difficult. This is due to the fact that z mightbe the only solution in ∆ of F with x = z ,i , but there is a another solution of F with x = z i, that is not contained within ∆ . We aim to address this problem by the followingapproach (see also Algorithm 22):Let L ∈ N be fixed non-negative integer. In a first step, we check whether (2424) is fulfilled.If this is not the case, we return False , otherwise, we proceed. Now, for each solution z i,j ∈ Z and each (cid:96) ∈ { , } , we use a complex root finder to compute a set of pairwise disjoint discs D (cid:96),j of radius less than − ρ such that each disc contains at least one root and the union ofall discs D (cid:96),j contains all complex roots of f (cid:96) ( z ,i , x ) ∈ C [ x ] . Then, we determine all discs D ,j , . . . , D ,j s that have a non-empty intersection with one of the discs D ,j (cid:48) . It followsthat each common root of f ( z ,i , x ) and f ( z ,i , x ) must be contained in one of the discs D ,j s (cid:48) . Hence, if each of these discs is contained in ∆ ,i , then x = z ,j is the unique solutionof f ( z ,i , x ) = f ( z ,i , x ) = 0 , and thus (2525) is fulfilled. In this case, we may concludethat µ ( z ,i , R ) equals the multiplicity of z . If we succeed in computing the multiplicities forall solutions in Z , we return the solutions together with their corresponding multiplicities.Otherwise, we return False .Obviously, the above approach cannot succeed if one of the above conditions is not fulfilled.However, even if both conditions are fulfilled, it may still fail due to the fact that ρ has notbeen chosen large enough. Lemma 11. Suppose that both conditions (2424) and (2525) are fulfilled. Then, there exists a L ∈ N such that Algorithm 22 succeeds for all L > L . Each of the methods in [BSS+16BSS+16; BSS+15BSS+15; MSW15MSW15] applies to polynomials with arbitrary complex coef-ficients. Also, for computing only approximations of the roots, there are no restrictions on the multiplicitiesof the roots. In our implementation, we use a strongly simplified variant of the algorithm from [BSS+16BSS+16]. lgorithm 2: BisolvePlus Input : Zero-dimensional bivariate system F : f ( x , x ) = f ( x , x ) = 0 with f , f ∈ Q [ x ] , polydisc ∆ r ( m ) Output: Z + such that Z + = { ( z , k ) : f ( z ) = f ( z ) = 0 , k = µ ( z , F ) } . for ρ := 1 , , , . . . . do Choose a matrix S ∈ S D F and compute F ∗ = ( f ∗ , f ∗ ) = F ◦ S − if deg x f ∗ = deg f ∗ and deg x f ∗ = deg f ∗ then Call Bisolve with input F ∗ and ∆ r ( m ) to compute (for (cid:96) = 1 , ): • Discs ∆ (cid:96),i , i = 1 , . . . , i (cid:96) , that isolate the roots z (cid:96),i of R (cid:96) := Res( F ∗ , x (cid:96) ) . • The multiplicity k (cid:96),i = µ ( z (cid:96),i ) of z (cid:96),i as a root of R (cid:96) . • The set Z := { z i,j = ( z ,i , z ,j ) : f ∗ ( z i,j ) = f ∗ ( z i,j ) = 0 } of all solutions of F ∗ . for each solution z = z i,j ∈ Z dofor (cid:96) = 1 , do Compute disjoint discs D (cid:96), , . . . , D (cid:96),s (cid:96) of radius less than − ρ such that • Each disc D (cid:96),s contains at least one root of f ∗ (cid:96) ( z i , x ) ∈ C [ x ] . • (cid:83) s (cid:96) s =1 D (cid:96),s contains all roots of f ∗ (cid:96) ( z i , x ) ∈ C [ x ] .Determine the set D ∗ := { D ,s : ∃ D ,s (cid:48) with D ,s ∩ D ,s (cid:48) (cid:54) = ∅} of all discs D ,s that have non-empty intersection with one of the discs D ,s (cid:48) . if for all D ,s ∈ D ∗ it holds that D ,s ⊂ ∆ i then Set b i,j = True if (cid:86) i,j b i,j thenreturn { ( S − ◦ z i,j , k i ) : z i,j ∈ Z and z i,j ∈ ∆ r ( m ) } Proof. Let (cid:15) be a lower bound on the distance between any distinct roots of f ( z ,i , x ) and f ( z ,i , x ) . Now, if − L < (cid:15)/ , then two discs D ,j (cid:48) and D ,j (cid:48)(cid:48) can only intersect if theycontain a common root of f ( z i, , x ) and f ( z i, , x ) . Since z ,j is the only common root, wethus conclude that each of the discs D ,j s (cid:48) must contain z ,j . Hence, if L is large enough, then ∆ contains D ,j s (cid:48) .The problem with this approach is that we do neither know in advance whether the con-dition (2525) is fulfilled nor do we know whether ρ has been chosen sufficiently large. In orderto overcome this issue, we consider a rotation of the system by means of a rotation matrix S ∈ S D F . Then, with probability at least / , both conditions (2424) and (2525) are fulfilled forthe rotated system F ∗ := F ◦ S − . We now proceed in rounds (numbered by m ), where,in each round, we choose a matrix S ∈ S D F at random and run Algorithm 22 with input F ∗ = F ◦ S − and ρ := 2 m . Since there are only finitely many different choices for S and since33 lgorithm 3: Validate Input : Zero-dimensional system F = ( f , f ) with polynomials f , f ∈ Z [ x ] ofdegrees d and d , respectively, polydisc ∆ = ∆ × ∆ ⊂ C , candidate z = ( z , z ) and multiplicities k (cid:96) s.t. the multiplicity of z (cid:96) as a root of R (cid:96) ( x ) := Res( F , x (cid:96) ) is k (cid:96) for (cid:96) ∈ { , } . Output: k ∈ N . If k ≥ , then there is a unique solution z of F = 0 of multiplicity k within the polydisc ∆ . Otherwise, ∆ contains no solution. for L = 1 , , , . . . do Compute L -bit approximation ζ of z if ∆ − L +2 (cid:100) log(min( k ,k (cid:101) +6 ( ζ ) ⊂ ∆ thenif / ∈ (cid:3) f ( ∆ − L ( ζ )) or / ∈ (cid:3) f ( ∆ − L ( ζ )) thenreturn for K = 1 , , , . . . , (cid:100) log(min( k ,k )) (cid:101) doif PolySol ( F , ∆ − L ( ζ ) , K ) = k ≥ thenreturn k the conditions (2424) and (2525) are fulfilled for at least the half of the systems F ∗ , Lemma 1111guarantees that, for sufficiently large ρ , Algorithm 22 returns the solutions of F ∗ in ∆ togetherwith the corresponding multiplicities with probability at least / . New Validation Phase. We are now ready to modify Bisolve by using an inclusionpredicate based on the algorithm PolySol . More specifically, let C := { z i,j } i,j be the setof candidate solutions and let ∆ i,j = ∆ ,i × ∆ ,j be the corresponding candidate regions ascomputed in the projection phase of Bisolve . The validation routine, see also Algorithm 33,that is called for each candidate solution z = z i,j = ( z ,i , z ,j ) again works in rounds, wherein round m , we compute a L = 2 m -bit approximation ζ = ( ζ , ζ ) of z such that (cid:107) ζ − z (cid:107) < − L and ∆ 64 max( d ,d ) · − L ( ζ ) ⊂ ∆ ,j . The exclusion predicate is identical to the original Bisolve routine, i.e., we check whether we can guarantee that f or f does not vanish on ∆ − L ( ζ ) by evaluating interval extensions (cid:3) f (cid:96) ( ∆ − L ( ζ )) for (cid:96) = 1 , using interval arithmetic.The inclusion predicate now works as follows. There is still one tiny detail that preventsus from directly plugging in PolySol as an inclusion predicate. Namely, even if thecandidate solution would actually turn out to be a solution of the system, we do not knowthe multiplicity k of this solution. Thus, we have to search for the multiplicity k . As wehave seen in the previous section that PolySol ( F , ∆ − L ( ζ ) , K ) actually succeeds for any K ≥ k , we can use exponential search for k by calling PolySol ( F , ∆ − L ( ζ ) , K ) for K =1 , , , . . . , (cid:100) log(min( k ,k )) (cid:101) , where k (cid:96) for (cid:96) ∈ { , } is the multiplicity of z (cid:96) as a root of R (cid:96) ( x ) :=Res( F , x (cid:96) ) . In the calls to PolySol , we use the above described BisolvePlus -routine inorder to implement the computation of the solutions of the truncated system in Line 44 ofAlgorithm 11. We remark that our actual implementation of the new inclusion predicate differsslightly from the description of PolySol in one more detail. For efficiency reasons, weconsider a partial change of order of the three considered steps Solving the truncated system , Projection step , and Bound Computation and Comparison .34 .1 Setting We performed experiments on a compute server with 48 Intel (R) Xeon (R) CPU E5-2680 v3@ 2.50GHz cores and a total of 256 GB RAM running Debian GNU/Linux 8. All code wasimplemented in SageMath version 7.6, release date 2017-03-25. The instances on which we compared the implementations are generated as follows. Givena trivariate polynomial P ∈ Z [ x, y, z ] . There are several different ways of obtaining twobivariate polynomials f, g from P that have solutions of higher multiplicity. The differentways are encoded by the strings , , , x0y , y0x in the file names. The followingtable summarizes the meaning of these abbreviations. We denote p v = ∂ v p for any polynomial p ∈ R [ v ] for some ring R . f = Res( P, P z , z ) g = f x · f x f = Res( P, P z , z ) g = f x · f y f = Res( P, P z , z ) g = f y · f y x0y f = Res( P, P z , z ) · f x g = f y y0x f = Res( P, P z , z ) · f y g = f x From the resulting system f, g , we construct the sheared system f, g ← f ( ax + by, cx + dy ) , g ( ax + by, cx + dy ) with integers a, b, c, d drawn uniformly at random from [ − , . Thisis done in order to make degenerate situations where multiple solutions share the same x or y -value less likely. We create an even larger set of instances by renaming the variables of P from x, y, z to x, z, y or y, z, x (or equivalently considering P x and P y instead of P z ). Weabbreviate this choice with xyz , xzy , and yzx . Now, let z be a solution of such a system f, g of multiplicity k . We pick random polynomials p, q of increasing degrees and considerthe systems f · p, g · q . This results in systems f d , g d of increasing degrees d that have thesame solution z of multiplicity k . For each degree d , we create three such system f d , g d bymultiplying f, g with different random polynomials.There are two different classes of instances that we consider depending on how the initialtrivariate polynomial P is chosen. In the first class, called herwig_hauser , we pick thepolynomial P from the set of polynomials given as three dimensional surfaces in the HerwigHauser Classics gallery [HauHau]. In the second class, called random , we pick P randomly. In thefirst class called herwig_hauser we let d = 10 , , . . . , , whereas in the second class random ,we let d = 16 , , . . . , . We note that in the latter case we pick the random polynomialwith which we multiply f, g in order to get f d , g d as sparse polynomials as otherwise evaluating f, g already becomes non-trivial.The generated instances can be found on the project page. A folder corresponding toa candidate contains one file called orig.cnd , which refers to the polynomials f, g . Theremaining files correspond to the polynomials f d , g d as described above. Every file containsfour lines, the first two contain the system, while the third and fourth contain the boundaries x − r, x + r and y − r, y + r such that the solution is contained within this range. http://resources.mpi-inf.mpg.de/systemspellet/http://resources.mpi-inf.mpg.de/systemspellet/ .11.010.0 10 20 30 40 degree t i m e Method, k=1 standardtruncate 100010000 10 20 30 40 degree p r e c i s i on Method, k=1 standardtruncate0.11.010.0 10 20 30 40 degree t i m e Method, k=2 standardtruncate 100010000 10 20 30 40 degree p r e c i s i on Method, k=2 standardtruncate0.11.010.0 10 20 30 40 degree t i m e Method, k=4 standardtruncate 100010000 10 20 30 40 degree p r e c i s i on Method, k=4 standardtruncate0.11.010.0 10 20 30 40 degree t i m e Method, k=8 standardtruncate 100010000 10 20 30 40 degree p r e c i s i on Method, k=8 standardtruncate Figure 1: Evaluation for validation of k -fold roots, for k = 1 , , , . In all plots the degree is on thehorizontal axis. On the left the validation time is on the vertical axis and on the right the precisiondemand is on the vertical axis. The red dots correspond to the validation method of the original Bisolve routine, called standard . The blue dots correspond to the validation method that uses thenew inclusion predicate, called truncate . 12 10 20 30 40 degree p r e c i s i on multiplicity k 12 48 100010000 100 1000 degree p r e c i s i on multiplicity k 12 4 Figure 2: Comparison of the dependence of the precision demand on the degree for different values of k . For herwig_hauser -instances on the left, and for random -instances on the right. In the first experiment, we compare the running time as well as the precision demand of thetwo respective validation methods called standard for the method included in the original Bisolve routine and truncate for the method using the new inclusion predicate on theinstance class herwig_hauser . In Figure 11, we can see the evaluation for validating k -foldroots for k = 1 , , , . The measurements are repeated three times, for each method andsystem. This results in 9 measurements (3 different random polynomials, 3 different runs)per degree per method. On the left, the running times are on the vertical logarithmic axis,whereas the degree of the systems is on the linear horizontal axis. On the right, the precisiondemand is on the vertical logarithmic axis, whereas the degree of the systems is on the linearhorizontal axis. The error bars indicate 95%-confidence intervals.We can see a clear advantage for our new method truncate . On average over all instancesof degree 40, we obtain an improvement of a factor of . , . , . , . for k = 1 , , , in the precision demand. In Figure 22 on the left, we can see the precision demand for the herwig_hauser instances for different k = 1 , , , for the truncate method. We can see thatthe precision demand increases with k in a comparable amount as the theoretical worst-casebounds predict, namely, we can roughly see a quadratic dependence between the precisiondemand and the multiplicity k in Figure 22 on the left.In Figure 22 on the right, we can see results for the same experiment for the random instances. In this experiment, we only include the truncate method as the original methoddoes not scale well enough for solving instances of that degree. Here both axis are logarithmicand the degree goes up to 4096. Fitting a linear model to the data points leads an estimatefor the exponent of . ± . 05 0 . ± . , and . ± . for k = 1 , , . The coefficientsof determination lie above . in all three cases that is roughly 94% of the variance of thedata can be explained by the fitted power model. Thus, we may conjecture that the precisiondemand depends at most a linearly on d . We remark that the plot suggests that the impactof the degree d dominates over the impact of k for very large d as we cannot see a differencebetween the curves for different values of k for large d . We remark that the impact of k forsmall d explains the smaller exponent in the fitted linear model for k = 4 compared to k = 1 , . The source code, the statistical data underlying the plots, the instances, and the script used for enchmarking are available for download on the project page. References [AZG16] Sergei A. Abramov, Eugene V. Zima, and Xiao-Shan Gao, eds. Proceedings ofthe ACM on International Symposium on Symbolic and Algebraic Computation,ISSAC 2016, Waterloo, ON, Canada, July 19-22, 2016 . ACM, 2016. isbn : 978-1-4503-4380-0.[BCG+08] Michael A. Burr, Sung Woo Choi, Benjamin Galehouse, and Chee-Keng Yap.“Complete subdivision algorithms, II: isotopic meshing of singular algebraic curves”.In: Symbolic and Algebraic Computation, International Symposium, ISSAC 2008,Linz/Hagenberg, Austria, July 20-23, 2008, Proceedings . 2008, pp. 87–94 (cit. onpp. 22, 66).[BEK+13] Eric Berberich, Pavel Emeliyanenko, Alexander Kobel, and Michael Sagraloff.“Exact symbolic-numeric computation of planar algebraic curves”. In: Theor.Comput. Sci. 491 (2013), pp. 1–32 (cit. on pp. 33, 55, 3131, 3232).[BHS+13] D.J. Bates, J.D. Hauenstein, A.J. Sommese, and C.W. Wampler. NumericallySolving Polynomial Systems with Bertini: Software, Environments, and Tools.Society for Industrial and Applied Mathematics, 2013. isbn : 9781611972696 (cit.on pp. 22, 66).[BS16] Cornelius Brand and Michael Sagraloff. “On the Complexity of Solving Zero-Dimensional Polynomial Systems via Projection”. In: Proceedings of the ACM onInternational Symposium on Symbolic and Algebraic Computation, ISSAC 2016,Waterloo, ON, Canada, July 19-22, 2016 . Ed. by Sergei A. Abramov, Eugene V.Zima, and Xiao-Shan Gao. ACM, 2016, pp. 151–158. isbn : 978-1-4503-4380-0(cit. on pp. 33, 55, 1111, 2020).[BSS+15] Ruben Becker, Michael Sagraloff, Vikram Sharma, and Chee-Keng Yap. “A Sim-ple Near-Optimal Subdivision Algorithm for Complex Root Isolation based onthe Pellet Test and Newton Iteration”. In: CoRR abs/1509.06231 (2015) (cit. onpp. 33, 44, 2222, 2323, 3232).[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 2016, Waterloo, ON, Canada, July 19-22, 2016 . Ed. bySergei A. Abramov, Eugene V. Zima, and Xiao-Shan Gao. ACM, 2016, pp. 71–78. isbn : 978-1-4503-4380-0 (cit. on p. 3232).[Blu98] L. Blum. Complexity and Real Computation . Springer New York, 1998. isbn :9780387982816 (cit. on p. 66).[Buc06] Bruno Buchberger. “Bruno BuchbergerâĂŹs PhD thesis 1965: An algorithm forfinding the basis elements of the residue class ring of a zero dimensional poly-nomial ideal”. In: Journal of Symbolic Computation issn : 0747-7171 (cit. on p. 55). http://resources.mpi-inf.mpg.de/systemspellet/http://resources.mpi-inf.mpg.de/systemspellet/ Using Algebraic Geometry . Graduate Textsin Mathematics. Springer New York, 2005. isbn : 9780387207063 (cit. on pp. 99,1010, 1313, 2222, 2828).[DE06] A. Dickenstein and I.Z. Emiris. Solving Polynomial Equations: Foundations,Algorithms, and Applications . Algorithms and Computation in Mathematics.Springer Berlin Heidelberg, 2006. isbn : 9783540273578 (cit. on p. 55).[DKS13] Carlos D’Andrea, Teresa Krick, and Martín Sombra. “Heights of varieties in mul-tiprojective spaces and arithmetic Nullstellensätze”. eng. In: Annales scientifiquesde l’École Normale Supérieure Symbolic and Algebraic Compu-tation, International Symposium, ISSAC 2010, Munich, Germany, July 25-28,2010, Proceedings . Ed. by Wolfram Koepf. ACM, 2010, pp. 243–250. isbn : 978-1-4503-0150-3 (cit. on p. 2424).[EP05] Ioannis Z. Emiris and Victor Y. Pan. “Improved algorithms for computing de-terminants and resultants”. In: J. Complexity Proceedings of the 2002 International Sym-posium on Symbolic and Algebraic Computation . ISSAC ’02. Lille, France: ACM,2002, pp. 75–83. isbn : 1-58113-484-3 (cit. on p. 55).[GHM+95] M. Giusti, J. Heintz, J. E. Morais, and L. M. Pardo. “When polynomial equa-tion systems can be “solved” fast?” In: Applied Algebra, Algebraic Algorithmsand Error-Correcting Codes: 11th International Symposium, AAECC-11 Paris,France, July 17–22, 1995 Proceedings . Ed. by Gérard Cohen, Marc Giusti, andTeo Mora. Berlin, Heidelberg: Springer Berlin Heidelberg, 1995, pp. 205–231. isbn : 978-3-540-49440-9 (cit. on p. 66).[Gal14] François Le Gall. “Powers of tensors and fast matrix multiplication”. In: Inter-national Symposium on Symbolic and Algebraic Computation, ISSAC ’14, Kobe,Japan, July 23-25, 2014 . Ed. by Katsusuke Nabeshima, Kosaku Nagasaka, FranzWinkler, and Ágnes Szántó. ACM, 2014, pp. 296–303. isbn : 978-1-4503-2501-1(cit. on p. 1111).[HL17] Jonathan D. Hauenstein and Viktor Levandovskyy. “Certifying solutions to squaresystems of polynomial-exponential equations”. In: Journal of Symbolic Compu-tation issn :0747-7171 (cit. on p. 66).[Hau] Herwig Hauser. https://imaginary.org/users/herwig-hauser (cit. on p. 3535).[Hoe11] J. van der Hoeven. Reliable homotopy continuation . Tech. rep. http://hal.archives-ouvertes.fr/hal-00589948/fr/ . HAL, 2011 (cit. on p. 66).[KS15] Alexander Kobel and Michael Sagraloff. “On the complexity of computing withplanar algebraic curves”. In: J. Complexity Symbolic and Algebraic Computation, International Sym-posium, ISSAC 2010, Munich, Germany, July 25-28, 2010, Proceedings . ACM,2010. isbn : 978-1-4503-0150-3.[Laz09] Daniel Lazard. “Thirty years of Polynomial System Solving, and now?” In: Jour-nal of Symbolic Computation issn : 0747-7171 (cit. on pp. 33, 55).[Llo75] N. G. Lloyd. “On Analytic Differential Equations”. In: Proceedings of the LondonMathematical Society s3-30.4 (1975), pp. 430–444. issn : 1460-244X (cit. on p. 2121).[MOS11] Kurt Mehlhorn, Ralf Osbild, and Michael Sagraloff. “A general approach to theanalysis of controlled perturbation algorithms”. In: Comput. Geom. J. Symb. Comput. J. Symb. Comput. 66 (2015), pp. 34–69 (cit. on pp. 2424, 2525, 3232).[Rou99] Fabrice Rouillier. “Solving Zero-Dimensional Systems Through the Rational Uni-variate Representation”. In: Applicable Algebra in Engineering, Communicationand Computing issn : 1432-0622 (cit. on pp. 33, 55).[Rum10] Siegfried M. Rump. “Verification methods: rigorous results using floating-pointarithmetic”. In: Symbolic and Algebraic Computation, International Symposium,ISSAC 2010, Munich, Germany, July 25-28, 2010, Proceedings . Ed. by WolframKoepf. ACM, 2010, pp. 3–4. isbn : 978-1-4503-0150-3 (cit. on p. 66).[Sto05] Arne Storjohann. “The shifted number system for fast linear algebra on integermatrices”. In: J. Complexity Theor. Comput. Sci. ACM Trans. Math. Softw. Fundamental Problems of Algorithmic Algebra . Oxford UniversityPress, 2000. isbn : 9780195125160 (cit. on p. 2424).[Zhi17] Lihong Zhi. “Computing Multiple Zeros of Polynomial Systems: Case of BreadthOne (Invited Talk)”. In: Computer Algebra in Scientific Computing: 19th Interna-tional Workshop, CASC 2017, Beijing, China, September 18-22, 2017, Proceed-ings . Ed. by Vladimir P. Gerdt, Wolfram Koepf, Werner M. Seiler, and Evgenii V.Vorozhtsov. Cham: Springer International Publishing, 2017, pp. 392–405. isbnisbn