Certified evaluations of Hölder continuous functions at roots of polynomials
Parker B. Edwards, Jonathan D. Hauenstein, Clifford D. Smyth
CCertified evaluations of H¨older continuousfunctions at roots of polynomials
Parker B. Edwards − − − ,Jonathan D. Hauenstein − − − , andClifford D. Smyth − − − Department of Applied and Computational Mathematics and Statistics, Universityof Notre Dame, Notre Dame, IN 46556 { parker.edwards,hauenstein } @nd.edu ∼ jhauenst Department of Mathematics and Statistics, University of North Carolina atGreensboro, Greensboro, NC 27402 [email protected] ∼ cdsmyth/ Abstract.
Various methods can obtain certified estimates for roots ofpolynomials. Many applications in science and engineering additionallyutilize the value of functions evaluated at roots. For example, critical val-ues are obtained by evaluating an objective function at critical points. Foranalytic evaluation functions, Newton’s method naturally applies to yieldcertified estimates. These estimates no longer apply, however, for H¨oldercontinuous functions, which are a generalization of Lipschitz continuousfunctions where continuous derivatives need not exist. This work developsand analyzes an alternative approach for certified estimates of evaluatinglocally H¨older continuous functions at roots of polynomials. An imple-mentation of the method in
Maple demonstrates efficacy and efficiency.
Keywords:
Roots of polynomials · H¨older continuous functions · Cer-tified evaluations.
For a univariate polynomial p ( x ), the Abel-Ruffini theorem posits that the rootscannot be expressed in terms of radicals for general polynomials of degree atleast 5. A simple illustration of this is that the solutions of the quintic equation p ( x ) = x − x − α -theory [9]. Kantorovich’s approach is based on bounds for a twice-differentiable function in an open set while Smale’s approach only uses localestimates at one point coupled with the analyticity of the function. Certified ap-proximations of roots of polynomials can also be obtained using interval methodssuch as [2,6,8] along with the Krawczyk operator [4,5]. a r X i v : . [ c s . S C ] J a n P.B. Edwards et al.
Although computing certified estimates for roots of polynomials is impor-tant, many applications in science and engineering utilize the roots in furthercomputations. As an illustrative example, consider the optimization problemmin { x − x − x + 3 : x ∈ R } . (2)The global minimum is the minimum of the critical values which are obtainedby evaluating the objective function g ( x ) = 21 x − x − x + 3 at its criticalpoints, i.e. at the real roots of g (cid:48) ( x ) = 168 x ( x − x − g (cid:48) ( x ), only approximations of the roots of g (cid:48) can be computed. Onemust translate these approximate roots to certified evaluations of the objectivefunction g ( x ) evaluated at the roots of the polynomial g (cid:48) ( x ) to obtain certifiedbounds on the global minimum of (2).One approach for computing a certified evaluation of f ( x ) at roots of a poly-nomial p ( x ) is via certified estimates of solutions to the multivariate system (cid:20) p ( x ) y − f ( x ) (cid:21) = 0 . (3)For sufficiently smooth f , approaches based on Newton’s method generate cer-tified estimates. The downside is that this requires f to be differentiable. Al-ternatively, one can follow a two-stage procedure: develop certified bounds of aroot of p ( x ) and then use interval evaluation methods, e.g., see [6, Chap. 5], todevelop certified bounds on f ( x ) evaluated at the root. This two-stage approachdoes not allow direct control on the precision of the certified evaluation.The approach in this paper considers certified evaluations of locally H¨oldercontinuous functions at roots of polynomials and links the desired output of thecertified evaluation with the error in the approximation of the root. H¨older con-tinuous functions are a generalization of Lipschitz functions which are indeedcontinuous, but they need not be differentiable anywhere. E.g. Weierstrass func-tion is locally H¨older. Satisfying the local H¨older continuity condition does notguarantee that a function can be evaluated exactly for, say, rational input. Ourapproach also incorporates numerical evaluation error into the certified boundsto address this issue.The rest of the paper is organized as follows. Section 2 describes the neces-sary analysis of locally H¨older continuous functions, with a particular focus onpolynomials and rational functions. Section 3 summarizes the approach used fordeveloping certified bounds on roots of polynomials. Section 4 combines the cer-tification of roots and evaluation bounds on H¨older continuous functions yieldingour approach for computing certified evaluations. Section 5 presents informationregarding the implementation in Maple along with several examples demonstrat-ing its efficacy and efficiency. Section 6 applies the techniques developed forcertified evaluations to prove non-negativity of coefficients arising in a seriesexpansion of a rational function. The paper concludes in Section 7.
The following describes the collection of functions under consideration. ertified evaluations 3
Definition 1.
A function f : C → C is locally H¨older continuous at a point x ∗ ∈ C if there exist positive real constants (cid:15), C, α such that | f ( x ∗ ) − f ( y ) | ≤ C · | x ∗ − y | α ≤ C · (cid:15) α (4) for all y ∈ B ( x ∗ , (cid:15) ) where B ( x ∗ , (cid:15) ) = { z ∈ C : | z − x ∗ | ≤ (cid:15) } . In this case, f ( x ) issaid to have H¨older constant C and H¨older exponent α at x ∗ . Moreover, if α = 1 ,then f ( x ) is said to be Lipschitz continuous at x ∗ with Lipschitz constant C . Functions which are locally H¨older continuous at a point are clearly con-tinuous at that point and the error bound provided in (4) will be exploited inSection 4 to provide certified evaluations. Every function f ( x ) which is contin-uously differentiable in a neighborhood of x ∗ is locally H¨older continuous with α = 1, i.e., locally Lipschitz continuous. For n ≥ f ( x ) = n (cid:112) | x | is continuousbut not differentiable at x ∗ = 0. It is locally H¨older continuous at x ∗ = 0 withH¨older constant C = 1 and H¨older exponent α = 1 /n .A computational challenge is to determine a H¨older constant C and H¨olderexponent α for f ( x ) on B ( x ∗ , (cid:15) ) given f ( x ), x ∗ , and (cid:15) >
0. Sections 2.1 and 2.2describe a strategy for polynomials and rational functions, respectively.
Since every polynomial f ( x ) is continuously differentiable, we can take the H¨olderexponent to be α = 1 at any point x ∗ . However, the H¨older constant C dependsupon x ∗ and (cid:15) . The Fundamental Theorem of Calculus shows that one just needs C ≥ max y ∈ B ( x ∗ ,(cid:15) ) | f (cid:48) ( y ) | . (5)Although one may attempt to compute this maximum, Taylor series expansionof f (cid:48) ( x ) at x ∗ provides an easy to compute upper bound. If d = deg f , f (cid:48) ( x ) = d (cid:88) i =1 f ( i ) ( x ∗ )( i − x − x ∗ ) i − , so that the triangle inequality yields C := d (cid:88) i =1 | f ( i ) ( x ∗ ) | ( i − (cid:15) i − ≥ max y ∈ B ( x ∗ ,(cid:15) ) | f (cid:48) ( y ) | . (6) The added challenge with a rational function f ( x ) is to ensure that it is definedon B ( x ∗ , (cid:15) ). One may attempt to compute the poles of f ( x ) and ensure thatnone are in B ( x ∗ , (cid:15) ), the implementation in Section 5 is based on the followinglocal approach that also enables computing local upper bounds on | f (cid:48) ( x ) | . For P.B. Edwards et al. f ( x ) = a ( x ) /b ( x ), one can prove b ( y ) (cid:54) = 0 for all y ∈ B ( x ∗ , (cid:15) ) by showing that | b ( x ∗ ) | > | b ( y ) − b ( x ∗ ) | for all y ∈ B ( x ∗ , (cid:15) ). If d b = deg b , then | b ( y ) − b ( x ∗ ) | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) d b (cid:88) i =1 b ( i ) ( x ∗ ) i ! ( y − x ∗ ) i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ d b (cid:88) i =1 | b ( i ) ( x ∗ ) | i ! (cid:15) i . Therefore, a certificate that f ( x ) is continuously differentiable on B ( x ∗ , (cid:15) ) is | b ( x ∗ ) | > d b (cid:88) i =1 | b ( i ) ( x ∗ ) | i ! (cid:15) i in which case min y ∈ B ( x ∗ ,(cid:15) ) | b ( y ) | ≥ | b ( x ∗ ) | − d b (cid:88) i =1 | b ( i ) ( x ∗ ) | i ! (cid:15) i > . (7)When b ( x ∗ ) (cid:54) = 0, it is clear that one can always take (cid:15) small enough to satisfy (7).When f ( x ) is continuously differentiable on B ( x ∗ , (cid:15) ), then one can take theH¨older exponent α = 1 and the H¨older constant C as in (5). Hence,max y ∈ B ( x ∗ ,(cid:15) ) | f (cid:48) ( x ) | ≤ max y ∈ B ( x ∗ ,(cid:15) ) | a (cid:48) ( y ) | min y ∈ B ( x ∗ ,(cid:15) ) | b ( y ) | + max y ∈ B ( x ∗ ,(cid:15) ) | a ( y ) | · max y ∈ B ( x ∗ ,(cid:15) ) | b (cid:48) ( y ) | min y ∈ B ( x ∗ ,(cid:15) ) | b ( y ) | where the maxima can be upper bounded similar to (6) and the minimum canbe lower bounded using (7). The initial task of determining certified evaluation bounds at roots of a givenpolynomial is to compute certified bounds of the roots. From a theoretical per-spective, we assume that we know the polynomial p ( x ) exactly. From a computa-tional perspective, we assume that p ( x ) has rational coefficients, i.e., p ( x ) ∈ Q [ x ].The certification of roots of p ( x ) can thus be performed using RealRootIsolate based on [1,7,10,11,12] in
Maple as follows.Since p ( x ) is known exactly, we first reduce down to the irreducible case withmultiplicity 1 roots by computing an irreducible factorization of p ( x ), say p ( x ) = p ( x ) r · · · p s ( x ) r s where p , . . . , p s are irreducible with corresponding multiplicities r , . . . , r s ∈ N .For p ( x ) ∈ Q [ x ], factor in Maple computes the irreducible factors in Q [ x ], i.e.,each p i ( x ) ∈ Q [ x ]. If z ∈ C is a root of p j ( x ), then z has multiplicity 1 withrespect to p j ( x ), i.e., z is a simple root of p j ( x ) with p j ( z ) = 0 and p (cid:48) j ( z ) (cid:54) = 0.In contrast, z has multiplicity r j with respect to p ( x ). ertified evaluations 5 For each irreducible factor q := p j , we transform the domain C into R via q ( x + iy ) = q r ( x, y ) + i · q i ( x, y ) where i = √− . (8)Therefore, solving q = 0 on C corresponds with solving q r = q i = 0 on R . Ap-plying RealRootIsolate with an optional absolute error bound abserr that willbe utilized later guarantees as output isolating boxes for every real solution to q r = q i = 0 on R . Therefore, looping over the irreducible factors of p , one ob-tains certified bounds for every root z of p ( x ) in C of the form a ≤ real( z ) ≤ a and b ≤ imag( z ) ≤ b where a , a , b , b ∈ Q .One can enhance this certification to provide additional information aboutthe roots. For example, suppose one aims to classify the roots of p ( x ) in C based on their modulus being less than 1, equal to 1, and greater than 1.The RealRootIsolate command can be applied to q r ( x, y ) = q i ( x, y ) = 0 for( x, y ) ∈ R together with x + y − < x + y − x + y − > Combining information on H¨older continuous functions from Section 2 and cer-tification of roots of polynomials from Section 3 yields the following approach todevelop certified evaluations. With input a polynomial p ( x ), a H¨older continuousfunction f ( x ) which is defined at each root of p ( x ), and an error bound (cid:15) > f ( z ) within (cid:15) for each root z of p ( x ). Since one may not be able to evaluate f ( x ) exactly, we incorporate anevaluation error of δ ∈ (0 , (cid:15) ). Typically, δ can be decreased by utilizing higherprecision computations. For rational input, rounding to produce finite decimalrepresentations constitutes the only source of representation error in Maple . Ourimplementation utilizes enough digits to have δ = (cid:15) .The first step is to utilize Section 3 to determine initial certified bounds foreach root of p ( x ). For an initial error bound on the roots, we start with γ = (cid:15)/ z is approximated by x ∗ with | z − x ∗ | < γ . Certified evaluationsare then obtained root by root since the H¨older constants are dependent uponlocal information near each root. In particular, the next step is to compute aH¨older exponent α and H¨older constant C that is valid on the ball B ( x ∗ , γ ). Ifthis is not possible, e.g., if f ( x ) cannot be certified to be defined on B ( x ∗ , γ ),one can simply reduce γ , e.g., by replacing γ by γ/
2, and repeat the processusing a newly computed certified approximation of z . Since f ( x ) is defined ateach root of p ( x ), this loop must terminate.The final step is to utilize local information to compute a new approximationof root z that will produce a certified evaluation within (cid:15) . Consider µ such that0 < µ ≤ min (cid:40) γ, α (cid:114) (cid:15) − δC (cid:41) and z ∗ an approximation of z such that z ∈ B ( z ∗ , µ ). Since | x ∗ − z ∗ | ≤ γ ,we have B ( z ∗ , µ ) ⊂ B ( x ∗ , γ ) so that all of the H¨older constants are valid on P.B. Edwards et al. B ( z ∗ , µ ). Hence, all that remains is to compute a certified approximation of f ( z ∗ ), say f ∗ , within the evaluation error of δ since | f ∗ − f ( z ) | ≤ | f ∗ − f ( z ∗ ) | + | f ( z ∗ ) − f ( z ) | ≤ δ + C · | z ∗ − z | α ≤ δ + C · µ α ≤ (cid:15). The certified evaluation procedure has been implemented as a
Maple pack-age entitled
EvalCertification available at https://github.com/P-Edwards/EvalCertification along with
Maple notebooks for the examples. The main ex-port is the procedure
EstimateRootsAndCertifyEvaluations with the follow-ing high level signature.
Input:–
Univariate polynomial p ∈ Q [ x ]. – List of locally H¨older continuous functions f , . . . , f m with which to certifi-ably estimate evaluations at the roots of p ( x ). – List of procedures specifying how to compute local H¨older constants andexponents for f , . . . , f m . (See Section 5.2 for example of the syntax). – Desired accuracy (cid:15) ∈ Q > . Main output:–
Complex rational root approximations z ∗ , . . . , z ∗ s , one for each of the distinctroots z , . . . , z s of p ( x ), such that | z j − z ∗ j | ≤ (cid:15) . – For each f i and x j , a complex decimal number f ∗ ij with | f i ( x j ) − f ∗ ij | ≤ (cid:15) .The EvalCertification package is formatted in a .mpl file which can be readinto a notebook with: read("EvalCertification.mpl")with(EvalCertification)
This lists the package’s following four exports: the main function and three builtin procedures for determining local H¨older constants and exponents for commonclasses of H¨older functions.
EstimateRootsAndCertifyEvaluations, HolderInformationForExponential,HolderInformationForPolynomial, HolderInformationForRationalPolynomial
The following examples highlight the specific
Maple types of these inputs andoutputs as well as other interface details.
As a first example, consider (2) by certifiably evaluating f ( x ) = 21 x − x − x + 3at the roots of p ( x ) = f (cid:48) ( x ) = 168 x ( x − x − ertified evaluations 7 f_polynomial := 21x^8 - 42x^4 - 56x^3 + 3;f_derivative := diff(f_polynomial, x);EstimationPrecision := 1/10^14 The main call to
EstimateRootsAndCertifyEvaluations is subsequently: solutions_information :=EstimateRootsAndCertifyEvaluations(f_derivative,[f_polynomial, f_derivative],HolderInformationForPolynomial,EstimationPrecision);
The first argument provides the polynomial to solve and the second a list ofpolynomials to evaluate. For illustration, we include evaluating the polynomialto solve in the evaluation list. The third argument is a procedure for computingH¨older constants which, in this case, uses the procedure that implements theestimates in Section 2.1 for polynomials. Notice that we need only provide theprocedure once since all functions for evaluation fall into the same class of H¨olderfunctions, namely polynomials. The last argument is the final error bound.The output solutions information is formatted as a Record. Certifiablyestimated roots are stored in a list as illustrated. solutions_information:-root_values =[0,2691619717901426047/2305843009213693952,26745188167908553113/147573952589676412928 -19995423894655642147*I/18446744073709551616, ...]
Evaluations are also stored in lists, one list for each function to evaluate withone entry for each root of p . Estimates are ordered so that the estimate at index i in its list corresponds to the root at index i in the roots list. solutions_information:-evaluations_functions_1 =[3., -91.6600084778015707, ...];solutions_information:-evaluations_functions_2 =[0, -6.692143197043304*10^(-16),...]; Therefore, the solution to (2) is − . − . Polynomial and rational functions can utilize the built-in procedures for comput-ing local H¨older constants. One more feature of
EvalCertification is the abilityto extend the certification procedures to new classes of functions by specifyinghow to compute local H¨older constants. As an example, consider f ( x ) = | x | r where r >
1. As noted in Section 2, f is globally H¨older continuous with H¨olderconstant C = 1 and H¨older exponent α = r . The following shows the necessaryformat for defining a new procedure to compute the local H¨older information. P.B. Edwards et al. alpha := 1/r;ExponentialInfo := proc(InputFunction,point,radius)return(Record("exponent"=alpha,"constant"=1,"domain_estimate"=false));end proc;
All custom H¨older information procedures must follow the same signature as thisexample. The exponent α and constant C in the output Record should satisfy theH¨older conditions for InputFunction on the ball B ( point , radius ) ⊂ C . Thedomain estimate lists points within radius of points missing from the inputfunction’s domain, or false if defined everywhere.To illustrate, consider f ( x ) = | x | / and f ( x ) = − x − x + 18 x + 88 x − x + 30 x + 16 x − x − x . The following commands produce certified evaluations of f ( x ) and f ( x ) at theroots of p ( x ) = 21 x − x − x +3 from Section 5.1 using the above procedure ExponentialInfo for r = 3. f_polynomial := 21x^8 - 42x^4 - 56x^3 + 3;f1 := abs(x)^(1/3);f2:=(-6x^5-8x^4+18x^3+88x^2-54)/(6x^5+30x^4+16x^3-18x^2-44x);EstimationPrecision := 1/10^14;solutions_information :=EstimateRootsAndCertifyEvaluations(f_polynomial,[f1,f2],[ExponentialInfo,HolderInformationForRationalPolynomial],EstimationPrecision); The third argument is a list instead of single procedure. This is necessary whenevaluating functions that require different approaches for computing H¨older in-formation. The built-in procedure
HolderInformationForExponential takesas input a single number α and outputs the procedure ExponentialInfo forthat α . We could equivalently replace ExponentialInfo with the procedure
HolderInformationForExponential(1/3) in this example.
The dominant computational cost is in estimating roots and computing localH¨older constants. Suppose that R ( p, (cid:15) ) is the complexity of approximating rootsof p within (cid:15) , H ( f , . . . , f n , p, (cid:15) ) is the minimum complexity of computing H¨olderconstants at one root, and A ( p, f , . . . , f n , (cid:15) ) is the number of repetitions requiredto find an accuracy γ ≤ (cid:15) where local H¨older constants can be calculated. Then A ( p, f , . . . , f n , (cid:15) )( R ( p, (cid:15) ) + n deg( p ) H ( f , . . . , f n , p, (cid:15) )) + R ( p, (cid:15) ) ertified evaluations 9 (a) (b) (c) (d) Fig. 1.
Results of tests with (a) d ∈ { , . . . , } , D = 5, n = 1, and (cid:15) = 10 − ; (b) d = 5, D ∈ { , . . . , } , n = 1, and (cid:15) = 10 − ; (c) d = 5, D = 5, n ∈ { , . . . , } , and (cid:15) = 10 − ; (d) d = 5, D = 5, n = 1, and (cid:15) ∈ { , − , . . . , − } . is a lower bound on the complexity. The number of repetitions A depends onthe input in a complicated way which we do not attempt to characterize here.We benchmarked EvalCertification using random polynomials generatedby the command randpoly in Maple . All tests computed roots of a randompolynomial p ( x ) with integer coefficients between − and 10 and evaluatedrational functions where the numerator and denominator were polynomials ofdegree D . The average was taken over 50 random selections. Figure 1 showsthe results of the benchmarking tests, which were performed on Ubuntu 18.04running Maple d of p ( x ), the value of D , the number of functions n to evaluate, andthe size of the output error (cid:15) . One application of our approach to certified evaluations is to certifiably de-cide whether or not all coefficients of the Taylor series expansion centered atthe origin are non-negative for a given real rational function r ( x ). We focuson non-negativity since non-positivity is equivalent to non-negativity for − r ( x )and alternating in sign is equivalent to non-negativity for r ( − x ). The followingmethod uses certified evaluations to obtain information about the coefficientsin the tail of the Taylor series expansion reducing the problem to only needingto inspect finitely many coefficients. This approach assumes that the functiondoes not have a pole at the origin, its denominator has only simple roots, and itsdenominator has a real positive root that is strictly smallest in modulus amongstall its roots. This approach can be extended to more general settings, but willnot considered here due to space considerations.We will make use of the following standard theorem. Theorem 1.
Let p ( x ) , q ( x ) ∈ R [ x ] such that p ( x ) and q ( x ) have no commonroot, q (0) (cid:54) = 0 and deg( p ( x )) < deg( q ( x )) = d . If q ( x ) has only simple roots say α , . . . , α d ∈ C , then r ( x ) = p ( x ) /q ( x ) has a Taylor series expansion of theform r ( x ) = (cid:80) ∞ n =0 r n x n converging for all x ∈ C with | x | < min {| α | , . . . , | α d |} .Furthermore, for all n ≥ , r n = − d (cid:88) i =1 p ( α i ) α i q (cid:48) ( α i ) α − ni . (9)Theorem 1 follows from partial fraction decomposition of rational functions orusing linear recurrences. For completeness, we provide a proof in the Appendix.Using Theorem 1, we obtain the following result on the eventual behavior of thecoefficients of the Taylor series of certain rational functions. Theorem 2.
With the setup from Theorem 1, define C i = − p ( α i ) / ( α i q (cid:48) ( α i )) for i = 1 , . . . , d . If α ∈ R is such that | α | < min {| α | , . . . , | α d |} , then thereexists N after which exactly one of the following conditions on r n holds:1. If α > and C > , then r n > for all n > N .2. If α > and C < , then r n < for all n > N .3. If α < , then r n is alternating in sign for all n > N , i.e., ( − n · r n > for all n > N or ( − n · r n < for all n > N .Moreover, one may take N = log( K ) / log( M/m ) where K = (cid:80) di =2 | C i | / | C | , m = | α | , and M = min {| α | , . . . , | α d |} . A proof of Theorem 2 is provided in the Appendix. Theorems 1 and 2 yieldthe following.
Corollary 1.
Suppose that f ( x ) , q ( x ) ∈ R [ x ] have no common root, q (0) (cid:54) = 0 ,and q ( x ) has only simple roots, namely α , . . . , α d ∈ C , such that α ∈ R and | α | < min {| α | , . . . , | α d |} . Let g ( x ) , p ( x ) ∈ R [ x ] be the unique polyno-mials such that f ( x ) = q ( x ) · g ( x ) + p ( x ) with deg( p ( x )) < deg( q ( x )) . De-fine C i = − p ( α i ) / ( α i q (cid:48) ( α i )) for i = 1 , . . . , d . Then, f ( x ) /q ( x ) has a Tay-lor series expansion f ( x ) /q ( x ) = (cid:80) ∞ n =0 R n x n converging for all x ∈ C with | x | < min {| α | , . . . , | α d |} and there is a threshold N so that exactly one of thefollowing conditions on R n holds:1. If α > and C > , then R n > for all n > N .2. If α > and C < , then R n < for all n > N .3. If α < , then R n is alternating in sign for all n > N , i.e. ( − n · R n > for all n > N or ( − n · R n < for all n > N .One can take N = max { deg( f ( x )) − deg( q ( x )) + 1 , log( K ) / log( M/m ) } where K = (cid:80) di =2 | C i | / | C | , m = | α | , and M = min {| α | , . . . , | α d |} .Proof. Since f ( x ) /q ( x ) = g ( x ) + p ( x ) /q ( x ), applying Theorem 1 yields the firstpart. Since the Taylor series coefficients of f ( x ) /q ( x ) and p ( x ) /q ( x ) are same for n > deg( f ( x )) − deg( g ( x )), the second part immediately follows from Theorem 2. ertified evaluations 11 One key to utilizing Theorem 2 and Corollary 1 is to certify that q ( x ) sat-isfies the requisite assumptions. Validating that q ( x ) has only simple roots fol-lows from computing an irreducible factorization as in Section 3 and checkingif every factor has multiplicity 1. Section 6.1 describes a certified approach toverify the remaining conditions on q ( x ). Section 6.2 yields a complete algorithmfor certifiably deciding non-negativity of all Taylor series coefficients which isdemonstrated on two examples. Given a polynomial q ( x ) ∈ R [ x ] with only simple roots and q (0) (cid:54) = 0, the fol-lowing describes a method to certifiably determine if q ( x ) has a positive rootthat is strictly smallest in modulus amongst all its roots. This method usesthe ability to certifiably approximate all real points in zero-dimensional semi-algebraic sets. Computationally, this can be accomplished using the command RealRootIsolate in Maple .The first step is to certifiably determine if q ( x ) has a positive root via P = { p ∈ R : q ( p ) = 0 , p > } . If P = ∅ , then one returns that q ( x ) does not have a positive root. Otherwise,one proceeds to test the modulus condition for α = min P .The modulus condition needs to be tested against negative roots and non-realroots. For negative roots, consider N = { n ∈ R : q ( − n ) = 0 , n > } and B = { b ∈ R : q ( b ) = q ( − b ) = 0 , b > } . By using certified approximations of α and points in N and B of decreasingerror, one can certifiably determine which of the following holds: α < min N , α > min N , or α ∈ B ⊂ P . If α > min N or α ∈ B , then one returns that q ( x )does has not a positive root that is strictly smallest in modulus amongst all itsroots. Otherwise, one proceeds to the non-real roots by considering L = { ( r, a, b ) ∈ R : q ( r ) = 0 , r > , q ( a + ib ) = 0 , b > , a + b < r } and E = { ( r, a, b ) ∈ R : q ( r ) = 0 , r > , q ( a + ib ) = 0 , b > , a + b = r } . Note that q ( a + ib ) = 0 provides two real polynomial conditions on ( a, b ) ∈ R via the real and imaginary parts as in (8) so that L and E are clearly zero-dimensional semi-algebraic sets. Moreover, for the projection map π ( r, a, b ) = r , π ( L ∪ E ) ⊂ P . By using certified approximations of α and points in L and E ofdecreasing error, one can certifiably determine if α ∈ π ( L∪E ) or α / ∈ π ( L∪E ).If the former holds, then one returns that q ( x ) does has not a positive root thatis strictly smallest in modulus amongst all its roots. If the later holds, then onereturns that q ( x ) does indeed have a positive root that is strictly smallest inmodulus amongst all its roots. Suppose that f ( x ) , q ( x ) ∈ R [ x ] which satisfy the assumptions in Corollary 1. Thefollowing describes a method to certifiably determine if all of the coefficients R n of the Taylor series expansion for f ( x ) /q ( x ) centered at the origin are non-negative or provides an integer n such that R n < g ( x ) , p ( x ) ∈ R [ x ] withdeg( p ( x )) < deg( q ( x )) such that f ( x ) = q ( x ) · g ( x ) + p ( x ). Define h ( x ) = x · q (cid:48) ( x )and C ( x ) = − p ( x ) /h ( x ). Hence, d = deg( q ( x )) = deg( h ( x )) such that q ( x )and h ( x ) have no common roots. As in Corollary 1, let α , . . . , α d be the rootsof q ( x ) with α ∈ R > such that α < min { α , . . . , α d } . Let β , . . . , β d ∈ C (notnecessarily all distinct) be the roots of h ( x ).Certified evaluations at the roots of q ( x ) and h ( x ) with error bound (cid:15) k = 2 − k for k = 1 , , . . . can be used until the following termination conditions are met:1. α ∗ i and β ∗ j are such that α ∗ ∈ R , | α ∗ i − α i | < (cid:15) k , and | β ∗ j − β j | < (cid:15) k
2. the set { , α ∗ , . . . , α ∗ d } is 2 · (cid:15) k separated, i.e., | s − t | ≥ (2 (cid:15) k ) for all distinct s, t in this set,3. γ ∗ ≤ min {| α ∗ i − β ∗ j | : 1 ≤ i, j ≤ d } such that γ ∗ > · (cid:15) k + (cid:15) / (4 d ) k ,4. for m ∗ = α ∗ + (cid:15) k and M ∗ ≤ min {| α ∗ | , . . . , | α ∗ d |} − (cid:15) k , one has m ∗ < M ∗ ,5. C ∗ i such that | C ∗ i − C ( α i ) | < (cid:15) k ,6. L ∗ i such that L ∗ i ≥ | c d | − (cid:80) d − (cid:96) =0 | u ( (cid:96) ) ( α ∗ i ) | (cid:15) (cid:96)k /(cid:96) ! where c d is the leading coef-ficient of q ( x ) and u ( x ) = − p (cid:48) ( x ) h ( x ) + p ( x ) h (cid:48) ( x ), and7. either (a) C ∗ + L ∗ √ (cid:15) k < C ∗ + L ∗ √ (cid:15) k > (cid:15) = (cid:15) k be the value where thetermination conditions are met and calculate 0 < b ∗ ≤ min {| C ∗ ± L ∗ √ (cid:15) |} , K ∗ ≥ (1 /b ∗ ) (cid:80) di =2 ( | C ∗ i | + L ∗ i √ (cid:15) ) , A ∗ ≥ log( K ∗ ) / log( M ∗ /m ∗ ) ,N ∗ = max { deg( f ( x )) − deg( q ( x )) , (cid:100) A ∗ (cid:101)} . The inequalities in Items 3, 4, and 6 and these values above are meant to signifythe rounding direction in machine precision used to compute the values. Withthis, if C ∗ + L ∗ √ (cid:15) k <
0, then return n = N ∗ + 1 in which R n <
0. Otherwise,the non-negativity of all R n is equivalent to the non-negativity of R , . . . , R N ∗ which can be computed by explicit computation. If there exists n ∈ { , . . . , N ∗ } such that R n <
0, return n . Otherwise, return that all R n are non-negative.Returning to the proof of correctness, note that Items 1, 5, and 6 follow fromevaluation errors and δ = min {| α i − α j | , | α i | : 1 ≤ i < j ≤ d } > . By Item 1, | α ∗ i − α ∗ j | ≥ | α i − α j | − (cid:15) k and | α ∗ i | ≥ | α i | − (cid:15) k . Thus, δ ∗ = min {| α ∗ i − α ∗ j | , | α ∗ i | : 1 ≤ i < j ≤ d } ≥ δ − (cid:15) k − η ertified evaluations 13 where η > δ ∗ .Since (cid:15) k → η can be made arbitrarily small, eventually δ ∗ − (cid:15) k − η > (cid:15) k and Item 2 will be met.Since q ( x ) and h ( x ) have no roots in common, consider γ = min {| α i − β j | : 1 ≤ i, j ≤ d } > . We have γ ∗ ≥ γ − (cid:15) k − η where η > γ ∗ . Eventually, γ ∗ ≥ γ − (cid:15) k − η > (cid:15) k + ν + (cid:15) / (4 d ) k where ν > (cid:15) / (4 d ) k and Item 3will be met.Since ∆ = M − m >
0, we have M ∗ > M − (cid:15) k − η where η > M ∗ and also m ∗ = α ∗ + (cid:15) k < α + 2 (cid:15) k . Thus, M ∗ − m ∗ ≥ ∆ − (cid:15) k − η so that eventually Item 4 will be met.Since C ( x ) = p ( x ) /h ( x ), we have C (cid:48) ( x ) = u ( x ) /h ( x ). Assuming the previousitems have all been met, we have | β j − α ∗ i | ≥ | β ∗ j − α ∗ i | − | β j − β J | ∗ > (cid:15) k + (cid:15) / (4 d ) k − (cid:15) k > (cid:15) k . Hence, h ( x ) has no roots in B ( α ∗ i , (cid:15) k ) and C (cid:48) ( x ) is continuous on B ( α ∗ i , (cid:15) k ). Fix z ∈ B ( α ∗ i , (cid:15) k ) and let ζ be the straight line segment contour from α ∗ i to z in B ( α ∗ i , (cid:15) k ). Since C (cid:48) ( x ) exists on B ( α ∗ i , (cid:15) k ), C ( z ) − C ( α ∗ i ) = (cid:82) ζ C (cid:48) ( x ) dx . Thus, | C ( z ) − C ( α ∗ i ) | ≤ P · | z − α ∗ i | where P = max {| C (cid:48) ( x ) | : x ∈ B ( α ∗ i , (cid:15) k ) } .Therefore, we have P ≤ P /P where P = max {| u ( x ) | : x ∈ B ( α ∗ i , (cid:15) k ) } and P = min {| h ( x ) | : x ∈ B ( α ∗ i , (cid:15) k ) } . Since u ( z ) = (cid:80) d − (cid:96) =0 u ( (cid:96) ) ( α ∗ i )( z − α ∗ i ) (cid:96) /(cid:96) !,we have P ≤ (cid:80) d − (cid:96) =0 | g ( (cid:96) ) ( α ∗ i ) | (cid:15) (cid:96)k /(cid:96) !(1 + η ) where η > | u ( (cid:96) ) ( α ∗ i ) | . Since h ( z ) = c d d (cid:89) j =1 ( z − β j ) and | z − β j | ≥ | α ∗ i − β ∗ j | − | z − α ∗ i | − | β ∗ j − β j | > (cid:15) / (4 d ) k , we have | h ( z ) | ≥ | c d | (cid:15) / . Thus, P ≥ | c d | (cid:15) / and | C ( z ) − C ( α ∗ i ) | ≤ L ∗ i √ (cid:15) .Since B ( α ∗ i , (cid:15) k ) ⊂ B ( α i , (cid:15) k ) ⊂ B ( α i , | u ( (cid:96) ) ( α ∗ i ) | will be boundedabove by the corresponding maximum values of | u ( (cid:96) ) ( z ) | for z ∈ B ( α , P varies in each step k , P and L ∗ will be uniformly bounded abovefor all k . By the proof of Theorem 2, C (cid:54) = 0 so Item 7(a) will eventually be metif C < C > N ∗ which is reasonable close to the value of N in Corollary 1, one may continue todecrease (cid:15) k past the point where all termination conditions are first met. Thereason for this is to separate m ∗ and M ∗ as far all possible, i.e., to match theactual gap between m and M as closely as possible. This could be wise especiallywhen m ∗ is very close to M ∗ in which case log( M ∗ /m ∗ ) will be very close to 0so that N ∗ will be very large. Example 1.
To demonstrate the approach, consider the rational functions(1 − x − x + x ) − and (1 − x − x + x ) − . The implementation of this approach in
Maple certifies that both rational func-tions have Taylor series expansions centered at the origin where all of the co-efficients are non-negative. The value of N from Corollary 1 which could becertified by the method described above was N ∗ = 204 and N ∗ = 55, respec-tively. Thus, it was easy to utilize series in Maple to check the non-negativityof the Taylor series coefficients up to N ∗ combined with Corollary 1 for the tail. This manuscript developed techniques for certified evaluations of locally H¨oldercontinuous functions at roots of polynomials along with an implementation in
Maple . These techniques were demonstrated on several problems including certi-fied bounds on critical values and proving non-negativity of coefficients in Taylorseries expansions. Although this paper focused on roots of univariate polynomi-als, it is natural to extend to multivariate polynomial systems in the future.
Acknowledgments
JDH was supported in part by NSF CCF 1812746. CDS was supported in partby Simons Foundation grant 360486.
References
1. Boulier, F., Chen, C., Lemaire, F., Maza, M.M.: Real root isolation of regularchains. In: The Joint Conference of ASCM 2009 and MACIS 2009, COE Lect.Note, vol. 22, pp. 15–29. Kyushu Univ. Fac. Math., Fukuoka (2009)2. Gargantini, I., Henrici, P.: Circular arithmetic and the determination of polynomialzeros. Numer. Math. , 305–320 (1971/72)3. Kantorovich, L.V.: On Newton’s method for functional equations. Doklady Akad.Nauk SSSR (N.S.) , 1237–1240 (1948)4. Krawczyk, R.: Newton-Algorithmen zur Bestimmung von Nullstellen mit Fehler-schranken. Computing (Arch. Elektron. Rechnen) , 187–201 (1969)5. Moore, R.E.: A test for existence of solutions to nonlinear systems. SIAM J. Numer.Anal. (4), 611–615 (1977)6. Moore, R.E., Kearfott, R.B., Cloud, M.J.: Introduction to interval analysis. Societyfor Industrial and Applied Mathematics (SIAM), Philadelphia, PA (2009)7. Rioboo, R.: Real algebraic closure of an ordered field: Implementation in axiom.In: Papers from the International Symposium on Symbolic and Algebraic Compu-tation. p. 206–215. ISSAC ’92, Association for Computing Machinery, New York,NY, USA (1992)8. Rump, S.M.: Verification methods: rigorous results using floating-point arithmetic.Acta Numer. , 287–449 (2010)ertified evaluations 159. Smale, S.: Newton’s method estimates from data at one point. In: The Mergingof Disciplines: New Directions in Pure, Applied, and Computational Mathematics(Laramie, Wyo., 1985), pp. 185–196. Springer, New York (1986)10. Xia, B., Yang, L.: An algorithm for isolating the real solutions of semi-algebraicsystems. J. Symbolic Comput. (5), 461–477 (2002)11. Xia, B., Zhang, T.: Real solution isolation using interval arithmetic. Comput. Math.Appl. (6-7), 853–860 (2006)12. Yang, L., Hou, X., Xia, B.: A complete algorithm for automated discovering of aclass of inequality-type theorems. Sci. China Ser. F (1), 33–49 (2001) Appendix
Proof of Theorem 1.
Suppose that C (cid:54) = 0 such that q ( x ) = C · (cid:81) di =1 ( x − α i ). Thus,we know q (cid:48) ( x ) = C · (cid:80) di =1 (cid:81) j (cid:54) = i ( x − α j ) and q (cid:48) ( α i ) = C · (cid:81) j (cid:54) = i ( α i − α j ) (cid:54) = 0 forall i . Let p i ( x ) = q ( x ) / ( x − α i ) = C · (cid:81) j (cid:54) = i ( x − α j ). Hence, p i ( α i ) = q (cid:48) ( α i ) and p i ( α j ) = 0 if j (cid:54) = i . The polynomials p , . . . , p d are linearly independent since,if (cid:80) di =1 a i p i ( x ) = 0, then evaluating at x = α j yields a j · q (cid:48) ( α j ) = 0 whichimplies a j = 0. Thus, they must form a basis for the d -dimensional vector spaceof polynomials of degree at most d − p ( x ) has degree at most d −
1, there are unique constants a i so that (cid:80) di =1 a i p i ( x ) = p ( x ). Evaluating at x = α j yields a j q (cid:48) ( α j ) = p ( α j ) so that a j = p ( α j ) /q (cid:48) ( α j ). Therefore, for all x ∈ C \ { α , . . . , α d } , p ( x ) q ( x ) = d (cid:88) i =1 p ( α i ) q (cid:48) ( α i ) 1 x − α i = d (cid:88) i =1 − p ( α i ) α i q (cid:48) ( α i ) 11 − x/α i . (10)The terms in (10) have a Taylor series expansion centered at the origin thatconverge for all x with | x | < min {| α | , . . . , | α d |} such that, as (9) claims, p ( x ) q ( x ) = d (cid:88) i =1 − p ( α i ) α i q (cid:48) ( α i ) ∞ (cid:88) n =0 α − ni x n = ∞ (cid:88) n =0 (cid:32) − d (cid:88) i =1 p ( α i ) α i q (cid:48) ( α i ) α − ni (cid:33) x n . Proof of Theorem 2.
Clearly, one has r n = d n dz n ( p ( z ) /q ( z )) | z =0 . Since p ( x ) and q ( x )have real coefficients, r n is real for all n ≥
0. For i ∈ { , . . . , d } , let t in = C i α − ni so that (9) reduces to r n = (cid:80) di =1 t in . Moreover, α ∈ R \{ } implies C ∈ R \{ } .Clearly, if α <
0, then t n is alternating in sign.Consider the case when α >
0. First, note that t n and C always have thesame sign. The following derives a threshold N such that | r n − t n | < | t n | for all n > N . Given such an N , r n will have the same sign as t n and C for n > N and the theorem will be proved. To that end, since ( r n − t n ) /t n = (cid:80) di =2 t in /t n , | r n − t n || t n | ≤ d (cid:88) i =2 | C i || C | | α | n | α i | n ≤ K (cid:16) mM (cid:17) n for all n . Since, by assumption, m/M <