> <\body> <\hide-preamble> >> >> \\\>> |||\>|\>>> \; |<\doc-note> Grégoire Lecerf and Arnaud Minondo have been supported by the French NODE project. Joris van der Hoeven has been supported by an ERC-2023-ADG grant for the ODELIX project (number 101142171). <\wide-tabular> ||||||| ||LOGO_ERC-FLAG_FP.png>|61.62pt|28.51pt||>>>>> ||<\author-affiliation> Laboratoire d'informatique de l'École polytechnique (LIX, UMR 7161 CNRS) CNRS, École polytechnique, Institut Polytechnique de Paris Bâtiment Alan Turing, CS35003 1, rue Honoré d'Estienne d'Orves 91120 Palaiseau, France |>>||<\author-affiliation> Laboratoire d'informatique de l'École polytechnique (LIX, UMR 7161 CNRS) CNRS, École polytechnique, Institut Polytechnique de Paris Bâtiment Alan Turing, CS35003 1, rue Honoré d'Estienne d'Orves 91120 Palaiseau, France |>>||<\author-affiliation> Laboratoire d'informatique de l'École polytechnique (LIX, UMR 7161 CNRS) CNRS, École polytechnique, Institut Polytechnique de Paris Bâtiment Alan Turing, CS35003 1, rue Honoré d'Estienne d'Orves 91120 Palaiseau, France |<\author-email> minondo@lix.polytechnique.fr >>||> How to automatically determine reliable error bounds for a numerical computation? One traditional approach is to systematically replace floating point approximations by intervals or balls that are guaranteed to contain the exact numbers one is interested in. However, operations on intervals or balls are more expensive than operations on floating point numbers, so this approach involves a non-trivial overhead. In this paper, we present several approaches to remove this overhead, under the assumption that the function that we wish to evaluate is given as a straight-line program (SLP). We will first study the case when the arguments of our function lie in fixed balls. For polynomial SLPs, we next consider the \Pglobal\Q case where this restriction on the arguments is removed. We will also investigate the computation of bounds for first and higher order derivatives of . |> Interval arithmetic is a popular technique to calculate guaranteed error bounds for approximate results of numerical computations. The idea is to systematically replace floating point approximations by small intervals around the exact numbers that we are interested in. Basic arithmetic operations on floating point numbers are replaced accordingly with the corresponding operations on intervals. When computing with complex numbers or when working with multiple precision, it is more convenient to use balls instead of intervals. In this paper, we will always do so and this variant of interval arithmetic is called . Unfortunately, ball arithmetic suffers from a non-trivial overhead: floating point balls take twice the space of floating point numbers and basic arithmetic operations are between two and approximately ten times more expensive. For certain applications, it may therefore be preferable to avoid the systematic use of balls for individual operations. Instead, one may analyze the error for larger groups of operations. For instance, when naively multiplying two double precision n> matrices whose coefficients are all bounded by in norm(say), then it is guaranteed that the maximum error of an entry of the result is bounded by *2>, without having to do any individual operations on balls. The goal of this paper is to compute reliable error bounds in a systematic fashion, while avoiding the overhead of ball arithmetic. We will focus on the case when the function that we wish to evaluate is given by a (SLP). Such a program is essentially asequence of basic arithmetic instructions like additions, subtractions, multiplications, and possibly divisions. For instance, n> matrix multiplication can be computed using an SLP. The SLP framework is actually surprisingly general: at least conceptually, the trace of the execution of a more general program that involves loops or subroutines can often be regarded as an SLP. So consider a function \\> that can be computed using an SLP, where =\> or =\>. Given \> and \>>, let > denote the ball with center and radius . We will consider several problems: Given an approximate evaluation ,\,b|)>\f,\,a|)>=f> using floating point arithmetic, can we efficiently compute a bound for the error? Given f> and ,\,r\\>>, can we efficiently compute ,\,s\\>> such that ,r|)>\\\\,r|)>|)>\\,s|)>\\\\,s|)>>? Given balls ,r|)>,\,\,r|)>>, an index ,n|}>>, and ,\,k>\,m|}>>, can we efficiently compute a bound for f|\ a>*\*\ a>>>|\|>> on ,r|)>\\\\,r|)>>? > Here \Pefficiently\Q means that the cost of the bound computation should not exceed >>, ideally speaking. In particular, the cost should not depend on the length of the SLP, but we do allow ourselves to perform precomputations with such a higher cost. We may regard as a special case of by taking =\=r=0>. If =0>, then becomes essentially a special case of , since we may use |\|>+s> as the required bound. Of course, there is a trade-off between the cost of bound computations and the sharpness of the obtained bounds. We will allow our bounds to be less sharp than those obtained using traditional ball arithmetic. But we still wish them to be as tight as possible under the constraint that the cost of the bound computations should remain small. We start with the special case when the centers ,\,a> lie in some fixed balls ,\,B>. For this purpose we introduce a special variant of ball arithmetic, called : see section. A matryoshka is a ball whose center is itself a ball and its radius is a number in >>. Intuitively speaking, a matryoshka allows us to compute enclosures for \Pballs inside balls\Q. By evaluating at matryoshki with centers ,\,B> and zero radii, we will be able to answer question : see section. Using this to compute bounds for the gradient of on \\\B>, we will also be able to answer question in the case when ,r|)>\B,\,\,r|)>\B>. In section, we will actually describe a particularly efficient way to evaluate SLPs at matryoshki. For this, we will adapt transient ball arithmetic from. This variant of ball arithmetic has the advantage that, during the computations of error bounds for individual operations, no adjustments are necessary to take rounding errors into account. We originally developed this technique for SLPs that only use ring operations. In section, we will extend it to SLPs that may also involve divisions. For polynomial SLPs that do not involve any divisions, we will show in section that it is actually possible to release the condition that ,\,a> must be contained in fixed balls ,\,B>. The idea is to first reduce to the case when the components :\\\>> with ,n>> are homogeneous. Then *a|)>=\>*f> for some >, so we may always rescale such that it fits into the unit poly-ball >. We may then apply the theory from sections and. Reliable numeric homotopy continuation is a typical application for which it is important to efficiently evaluate polynomial SLPs at arbitraryballs. Our final section is devoted to question . The main idea is to compute bounds for the |\|>> on balls ,r+\|)>,\,\,r+\|)>> with the same centers but larger radii. We will then use Cauchy's formula to obtain bounds for the derivatives of without having to explicitly compute these derivatives. Bounds for the derivatives of are in particular useful when developing reliable counterparts of Runge\UKutta methods for the integration of systems of ordinary differential equations. We intend to provide more details about this application in an upcoming paper. Throughout this paper, we assume that we work with a fixed floating point format that conforms to the IEEE 754 standard. We write for the bit precision, the number of fractional bits of the mantissa plus one. We also denote the minimal and maximal allowed exponents by > and >. For 754 double precision numbers, this means that , =-1022> and =1023>. We denote the set of hardware floating point numbers by >. Given an >-algebra >, we will also denote the corresponding approximate version by >. For instance, if =\=\|]>>, then we have =\|]>>. The IEEE 754 standard imposes correct rounding of all basic arithmetic operations. In this paper we will systematically use the mode. We denote by>> the result of rounding \> according to this mode. The quantity >\>-x|\|>>> stands for the corresponding rounding error, which may be>. Given a single operation >\,\|}>>>, it will be convenient to write >y> for y|)>>>. For compound expressions >, we will also write |]>> for the full evaluation of> using the rounding mode>. For instance, *b|]>=x>\>y>+>>\>a>|)>\>b>>. We denote by |\>>> any upper bound function for >> that is easy to compute. In absence of underflow, one may take |\>>=>|\|>*2>. If we want to allow for underflows during computations, then we can take |\>>=>|\|>*2+2-p+1>> instead, where -p+1>> is the smallest positive subnormal number in >. If \>, then we may still take |\>>y|)>=>y|\|>*\>> since no underflow occurs in that special case. See for more details about these facts. For bounds on rounding errors, the following lemma will be useful. <\lemma> Let \> and \0> be such that \\>. Then |)>\1+*\>. <\proof> Let \|i+1>> for all \>. Since \1>, we have <\eqnarray*> ||)>>||+|2>*\+3>*\>>||>|+*\|2>+*\|6>*\*\>>>||>|++|6>|)>*\>>||>|*\.>>>> Let > be an >-algebra and let |\|>> be a norm on >. We will typically take =\> or =\>, but more general normed algebras are also allowed. Given \> and \>, let >\\,\r|}>> be the closed ball with center and radius. We denote by ,\|)>> the set of all such balls. The aim of ball arithmetic is to provide a systematic way to bound the errors of numerical computations. The idea is to systematically replace numerical approximations of elements in > by balls that are guaranteed to contain the true mathematical values. It will be useful to introduce a separate notation > for this type of semantics: given aball \\,\|)>> and a number, we say that > if \x>, <\equation*> \\x\x\\. We also introduce poly-balls, which are vectors of balls, for situations where it is required to reason \Pcoordinate wise\Q. For all ,\,a|)>\\> and all ,\,r|)>\\>, we denote the poly-ball \,r|)>,\,\,r|)>|)>\\,\|)>>>. We extend this enclosure relation to poly-balls \\,\|)>> and \> as follows: <\equation> \\x\\\x\\\\\x. In the sequel, we will use a bold font for ball enclosures and a normal font for values at actual points. Note that the smallest enclosure of \> is the \Pexact\Q ball >, and we will regard > as being embedded into ,\|)>> in this way. Let us now recall how to perform arithmetic operations in a way that is compatible with the \Penclosure semantics\Q. Given a function \\\>, a of is a function \,\|)>\\,\|)>> that satisfies the <\equation*> \\x\f|)>\f for all \\,\|)>> and \>. (Note that we voluntarily used the same name for the function and its lift.) For instance, the basic arithmetic operations admit the following ball lifts: <\eqnarray*> \\>|>|b,r+s|)>>>|\\>|>|+r|)>*s+*r|)>.>>>> Other operations can be lifted in a similar way (in section below, we will in particular study division). We can also extend the norm from > to ,\|)>> via <\eqnarray*> |\|>\\,\|)>>|>|>>|>|>|+r.>>>> The extended norm is sub-additive, sub-multiplicative, and positive definite. \ For any ball =\>, we have +|)>=\+\=>\0>, so ,\|)>>> is clearly not an additive group in the mathematical sense. Nonetheless, we may consider ,\|)>>> to be a normed >-algebra from a computer science perspective, because all the relevant operations >, and |\|>> are \Pimplemented\Q. Allowing ourselves this shortcut, we may apply the theory from the previous subsection and formally obtain a ball arithmetic on ,\|)>,\|)>>. It turns out that this construction actually makes sense from a mathematical perspective, provided that we appropriately adapt the semantics for the notion of \Penclosure\Q. Let ,\,\|)>\\,\|)>,\|)>>. An element \\,r|)>> of ,\,\|)>> will be called a . Given a matryoshka =\>, we call its , its , and its . The appropriate enclosure relation for matryoshki is defined as follows: given a matryoshka =\> and aball =\>>, we define <\equation*> \\\\\\as\r. Conceptually, the \Psmall embedded ball\Q > is contained in the matryoshka, which corresponds to the \Pbig ball\Q >: see Figure. The enclosure relation naturally extends to vectors as in().|gr-frame|>|gr-geometry||gr-grid||gr-grid-old||1>|gr-edit-grid-aspect|||>|gr-edit-grid||gr-edit-grid-old||1>|gr-dash-style|10|gr-transformation||||>||>||||>||>>|||>|||>>||>>|>||>|>|>|>>> Representation of an embedded ball inside a matryoshka. > We use the abbreviation ,\,\|)>\\,\|)>,\|)>> for the set of matryoshki. Afunction ,\|)>\\,\|)>> is said to to ,\,\|)>\\,\,\|)>> if it satisfies the inclusion principle: <\equation*> \\\\f|)>\f|)> for all \\,\,\|)>> and \\,\|)>>. Exactly the same formulas() and() can be used in order to lift the basic arithmetic operations: <\eqnarray*> ,r|)>\\,s|)>>|>|\\,r+s|)>>>|,r|)>\\,s|)>>|>|\\,|\|>+r|)>*s+|\|>*r|)>,>>>> for all ,r|)>,\,s|)>\\,\,\|)>>. (Note that it was this time convenient to use >, > as a notation for the centers of our matryoshki.) We may also extend the norm on ,\|)>> to matryoshki via <\eqnarray*> |\|>\\,\,\|)>>|>|>>|,r|)>>|>||\|>+r>>>> which again remains sub-additive, sub-multiplicative, and positive definite. <\remark> In principle, we could repeat the process and recursively consider matryoshki as centers of even larger matryoshki. While attractive from a folkloric perspective, it turns out that the basic matryoshki are more efficient for the applications we know of. We will write ,\|)>> and ,\,\|)>> for the approximate versions of ,\|)>> and ,\,\|)>> when working with machine floating point numbers in > instead of exact real numbers in >. In that case, the formulas from section and need to be adjusted in order to take into account rounding errors. For instance, if =\>, then one may replace the formulas(), (), and() by <\equation> ||\>\>|>|>b,>|\>>b|)>|)>*2|)>|]>|)>>>|\>\>|>|>b,>+r|)>*s+*r+|\>>b|)>|)>*2|)>|]>|)>>>||\|>>>|>|>+r|)>*2|)>|]>.>>>>> See, ,. We will call this . Unfortunately, these formulas are far more complicated than(),(), and(), so bounding the rounding errors in this way gives rise to a significant additional computational overhead. An alternative approach is to continue to use the non-adjusted formulas <\equation> ||\>\>|>|>b,>|)>>>|\>\>|>|>b,>+r|)>*s+*r|]>|)>>>||\|>>>|>|>+r|]>.>>>>> This type of arithmetic was called in. This arithmetic is not a ball lift, but we will see below how to take advantage of it by inflating the radii of the input balls of a given SLP. Of course, we need to carefully examine the amount of inflation that is required to ensure the correctness of our final bounds. This will be the purpose of the next section in the case where we evaluate an entire SLP using transient ball arithmetic. We denote by >\\\*2|)>> a quantity that satisfies >\> and <\equation> |>b-a\b|\|>>|>|>b|\|>*\>>>|>-|\|>>|>|>*\>,>>>>> for all \> and >\|}>>, in absence of underflows and overflows. One may for instance take >\2> and >\4\2>, whenever 16>; seeA of the preprint version>. We recall the following lemma, which provides a useful error estimate for the radius of a transient ball product. <\lemma> 3>> For all \> and \> such that the computation of <\equation*> R=+r|)>*s+*r|]> involves no underflows or overflows, we have +r|)>*s+*r\R*>|)>>. In the presence of overflows or underflows, additional adjustments are required. Overflows actually cause no problems, because the IEEE 754 standard rounds every number beyond the largest representable floating point number to infinity, so the radii of balls automatically become infinite whenever an overflow occurs. In order to protect ourselves against underflows, the requirements() need to be changed into <\equation> ||>b|)>-b|)>|\|>>|>|>b|\|>*\>>>|>b-a\b|\|>>|>|>b|\|>*\>+\>>>|>-|\|>>|>|>*\>+\>.>>>>> Here we recall that additions and subtractions never lead to underflows. If =\>, then we may take >\2+1>>. If =\>, then a safe choice is >\8\2+1|)>/2>> (in fact, >\5\2+1>> is sufficient for multiplication). We refer to3.3> for details. Taking ,\|)>> instead of > for our center space, the formulas() naturally give rise to transient matryoshka arithmetic: <\eqnarray*> ||,r|)>\>\,s|)>>|>|\>\,>|)>>>|,r|)>\>\,s|)>>|>|\>\,>|\|>+r|)>*s+|\|>*r|]>|)>>>|,r|)>|\|>>>|>|>|\|>+r|]>.>>>> Here we understand the operations >>, >>, and |\|>>> on the centers are done use transient ball arithmetic. We will need an analogue of Lemma for matryoshki. Given balls =\>> and =\> in ,\|)>>, we define -\\\-\,|)>>. In what follows, let ,\|)>>\\\*2|)>> be such that ,\|)>>\> and <\equation> |\>\-\\\|\|>>|>|\>\|\|>*\,\|)>>>>||\|>>-|\|>|\|>>|>||\|>>*\,\|)>>,>>>>> for all ,\\\,\|)>> and >\|}>>, in absence of underflows and overflows. Let \max>,\>|)>>. Given \\\\,\|)>> and \\\\,\|)>>, Lemma gives +r|)>*s+*r\+r|)>*s+*r|]>*|)>>; similarly we have <\equation*> +r|)>*s+*r|]>*|)>\+r|)>*s+*r. This implies <\eqnarray*> \>\-\\\|\|>>|>|\>\|\|>*max|)>-1,1-|)>|)>>>||>|\>\|\|>*|)>.>>>> It thus suffices to take ,\|)>>\5*max>,\>|)>> for() to hold in the case of multiplication. Simpler similar computations show that this is also sufficient for the other operations. We can now state the analogue of Lemma. <\lemma> For all ,\\\,\|)>> and \> such that the computation of <\equation*> R=|\|>+r|)>*s+|\|>*r|]> involves no underflows or overflows, we have |\|>+r|)>*s+|\|>*r\R*,\|)>>|)>>. <\proof> The same as the proof of1>, . A > over a ring > is a sequence ,\,\> of instructions of the form <\eqnarray*> >|>|\C>>|>|>|\Y\Z,>>>> where ,Y,Z> are variables in a finite ordered set >, > constants in >, and >\|}>>>. Variables that appear for the first time in the sequence in the right-hand side of an instruction are called . Adistinguished subset of the set of all variables occurring at the left-hand side of instructions is called the set of . Variables which are not used for input or output are called and determine the amount of auxiliary memory needed to evaluate the program. The length >=l> of the sequence is called the of >. Let ,\,I> be the input variables of > and ,\,O> the output variables, listed in increasing order. Then we associate an >:\\\> to > as follows: given ,\,a|)>\\>, we assign > to > for ,m>, then evaluate the instructions of > in sequence, and finally read off the values of ,\,O>, which determine >,\,a|)>>. To each instruction >, one may associate the > as follows. Let =1>, and assume that ,\,q> have been defined for some ,l|}>>. Then we take =max >,\,q>|)>+1>, where \\\i> are those indices k> such that > is of the form \Y\Z>> with \,Z|}>>, >\|}>>>, and \,\,X|}>>. If no such indices exist, then we set =1>. Similarly, for each input variable > we define >=>,\,q>|)>>+1>, where \\\i> are those indices such that > is of the form \Y\Z> with \,Z|}>>>, >\|}>>>, and \,\,X|}>>. We also define >=>,\,q>,q>,\,q>|)>>>, where ,\,I>> are the input variables of > and ,\,i> all indices such that > is of the form\C>>. <\example> Let us consider =\5,x\a\a,x\x\x,x\x+a|)>>, of length. The input variables are > and >, and we distinguish > as the sole output variable. This SLP thus computes the function *a+a>. The associated computation graph, together with remaining path lengths are as pictured: |gr-frame|>|gr-geometry||gr-grid||gr-grid-old||1>|gr-edit-grid-aspect|||>|gr-edit-grid||gr-edit-grid-old||1>|gr-grid-aspect|||>|gr-grid-aspect-props|||>|||||||||||>||>||>|>|>|>|>|>|>>|>|>>||>>||>|=3>|>|=3>|>|=2>|>|=1>|>|>>> Our goal is to evaluate SLPs using transient ball arithmetic from section instead of rounded ball arithmetic. Transient arithmetic is faster, but some adjustments are required in order to guarantee the correctness of the end-results. In order to describe and analyze the correctness of these adjustments, it is useful to introduce one more variant of ball arithmetic. In , the operations on centers are approximate, but correctly rounded, whereas operations on radii are exact and certified: <\eqnarray*> ||>|)>\>\>|)>>|>|>b|]>,r>+s>+|\>>b|)>|)>>>|>|)>\>\>|)>>|>|>b|]>,+r>|)>*s>+*r>+|\>>b|)>|)>,>>>> where \> and >,s>\\>. The extra terms |\>>b|)>> in the radius ensure that these definitions satisfy the inclusion principle. In order to measure how far transient arithmetic can deviate from this idealized semi-exact arithmetic, let \+\+> be the -th harmonic number and let \H-H> for all k>. The following theorem shows how to certify transient ball evaluations: <\theorem> 5>> Let > be an SLP of length as above and let \0> be an arbitrary parameter such that \|)>>>>, where =\>> is as in. Consider two evaluations of> for different types of ball arithmetic. For the first evaluation, we use semi-exact ball arithmetic with |\>>=>|\|>*\>. For the second evaluation, we use transient ball arithmetic with the additional property that any input or constant ball >|)>>> is replaced by a larger ball > with <\equation*> r\max*|)>*q>>-1|)>,|)>*r>|)>, where <\equation*> \max|\>*\|)>,>||\H>>|)>>>*|1+\>*|)>>>|1+\>|)>.>>>>> Assume that no underflow or overflow occurs during the second evaluation. For any corresponding outputs >|)>> and> for the first and second evaluations, we then have>\t>. Given a fixed > and >=O|)>>, we note that we may take > and > such that |)>*q>>-1=O*q>*log q>|)>>. The value of the parameter > may be optimized as a function of the SLP and the input. Without entering details, > should be taken large when the input radii are small, and small when these radii are large. The latter case occurs typically within subdivision algorithms. For our main application when we want to compute error bounds, the input radii are typically small or even zero. For our implementations in and , we found it useful to further simplify the bounds from Theorem, assuming that > is \Psufficiently small\Q. <\corollary> With the notation from Theorem, assume that >|)>\\>. In order to apply Theorem, it is sufficient to take \\\>+1|)>*\> and <\equation*> \\max>+1|)>*|)>*|)>|\-\>|)>. <\proof> The harmonic series satisfies the well-known inequality >\H>>\ln q>+1>, ofwhich we only use the right-hand part. Lemma also yields |)>>>\1+\>. Given that \\>, it follows that <\equation*> |)>>>|1+\>|)>=|1+\-|)>>>>\|\-\>. Therefore, <\equation*> >>|)>>>*|)>>>|1+\>|)>>|>|>|)>+1|)>*|)>*|)>|\-\>>>>>>. The semi-exact ball arithmetic from the previous subsection naturally adapts to matryoshki as follows: given ,\\\,\|)>> and >,s>\\>, we define <\eqnarray*> |,r>|)>\>\,s>|)>>|>|\>\,r>+s>+|\>>\\|)>|)>,>>|,r>|)>\>\,s>|)>>|>|\>\,|\|>+r>|)>*s>+|\|>*r>+|\>>\\|)>|)>,>>>> where |\>>> is extended to balls <\eqnarray*> |\>>|)>>|>||\>>|\|>|)>.>>>> Recall that |\>>|\|>|)>> denotes an upper bound on the rounding errors involved in the computation of the norm of >. Operations on centers use the semi-exact arithmetic from the previous subsection. This arithmetic is therefore a matryoshka lift of the elementary operations. The following theorem specifies the amount of inflation that is required in order to certify SLP evaluations using \Ptransient matryoshka arithmetic\Q. <\theorem> Let =\,\>>, and let ,\,\> be as in Theorem. Consider two evaluations of> for two different types of arithmetic. For the first evaluation, we use semi-exact matryoshka arithmetic with |\>>=>|\|>*\>. For the second evaluation, we use transient matryoshka arithmetic with the additional property that any input or constant ball >,r>|)>> is replaced by a larger ball =\,r|)>> with <\equation*> ||>|*|)>*q>>-1|)>,|)>*R>|)>,>>||>||\|>*|)>*q>>-1|)>,|)>*r>|)>.>>>>> Assume that no underflow or overflow occurs during the second evaluation. For any corresponding outputs >,t>|)>> and> for the first and second evaluations, we then have>\T> and>\t>. <\proof> The theorem essentially follows by applying Theorem twice: once for the big and once for the small radii. This is clear for the big radii >> and since the formulas for the big radii of matryoshki correspond with the formula for the radii of the corresponding ball arithmetic. For the small radii, we use essentially the same proof as in. For convenience of the reader, we reproduce it here with the required minor adaptations. Let ,T,t|)>=\,t|)>> be the value of the variable > after the evaluation of ,\,\>>using transient matryoshki arithmetic. It will be convenient to systematically use star superscripts for the corresponding value ,T>,t>|)>=\>,t>|)>> when using semi-exact matryoshka arithmetic. Let us show by induction on that: <\equation> t\|\|>*|)>*q>-1|)>.> If > is of the form \\,t|)>>, then we are done since \q>>. Otherwise, > is of the form \Y\Z>. Writing =\,r|)>>, we claim that <\equation*> r\|\|>*|)>*+1|)>>-1|)>. This holds by assumption if > is an input variable. Otherwise, let k> be the largest index such that =X>. Then \q+1> by construction of >, whence our claim follows from the induction hypothesis. Similarly, writing =\,s|)>>, we have |\|>*|)>*+1|)>>-1|)>>>. Having shown our claim, let us first consider the case where >\>. In this case we obtain <\eqnarray*> ||>||\|>+|\|>|)>*|)>*+1|)>>-1|)>>>||>|\\|\|>*|)>*+1|)>>-1|)>.>>>> Combined with and the inequalities <\equation> ||>|>||)>>>||)>-1|)>*|)>>|>||)>-1,>>>>> we deduce: <\eqnarray*> |>||>>>||>|*|)>>>||>|\\|\|>*|)>*+1|)>>-1|)>*|)>>>||>|\>\|\|>*|)>*+1|)>>-1|)>*|)>*|)>>>||>|\>\|\|>*|)>*+1|)>>-1|)>*|)>>>||>|\>\|\|>*|)>*q>-1|)>.>>>> In case when >=>>, we obtain <\eqnarray*> |||\|>*s+|\|>*r+r*s>|>||\|>*s+|\|>*r>>||>|\\|\|>*|)>*+1|)>>-1|)>.>>>> Combined with Lemma, the inequalities, , and |)>\|)>> successively imply <\eqnarray*> ||>|||\|>+r|)>*s+|\|>*r|]>>>||>||\|>*s+|\|>*r+r*s|)>*|)>>>||>|\\|\|>*|)>*+1|)>>-1|)>*|)>>>||>|\>\|\|>*|)>*+1|)>>-1|)>*|)>*|)>>>||>|\>\|\|>*|)>*q>-1|)>,>>>> which achieves the induction. Then, for all ,l|}>>, we define <\equation*> \\>+|\*|)>>*H,q>>|)>*|)>>-q|)>> so that |)>\\>. By the inequality that we imposed on>, we have\1>>. Using asecond induction over , let us next prove that <\equation> t>\\*t.> Assume that this inequality holds up until order . If > is of the form \\,t|)>> then we are done by the fact that \|)>>. If > is of the form \Y\Z>, then let k> be the largest integers so that =Y> and =Z>, so we have ,q|)>\q+1>, <\equation*> max,\|)>\>+|\*|)>>*H+1,q>>|)>*|)>>-+1|)>|)>>, and <\equation*> \\|\|>*\|t>\||)>*q>-1|)>>\*q>\|q*\*|)>>, using. In particular, we obtain ,\|)>+\|)>*|)>\\>. We denote by >,r>|)>> (resp. >,s>|)>>) the matryoshka corresponding to the value of > (resp. >) in the semi-exact evaluation. Their corresponding values for the transient evaluation are denoted without the star superscript. Thanks to Theorem, we know that >\T>, whence >|\|>\|\|>> since >=C>. With and as above, using and, it follows that <\eqnarray*> |>>||>+s>+>|\|>*\>>||>|>+s>+|\|>*\>>||>|*max,\|)>+\*t|)>>>||>|,\|)>+\|)>*t*|)>>>||>|*t.>>>> If > is of the form \Y\Z>, then, thanks to Lemmas, , and inequality a similar calculation yields <\eqnarray*> ||>>||>|\|>*s>+>|\|>+s>|)>*r>+>|\|>*\>>||>||\|>*s>+|\|>+s|)>*r>+\*t>>||>||\|>*s+|\|>+s|)>*r|)>*max,\|)>+\*t>>||>|,\|)>+\|)>*t*|)>>>||>|*t,>>>> which completes the second induction and the proof of this theorem. <\remark> InIII.C>, it is discussed how to manage underflows and overflows in Theorem. A similar discussion applies to Theorem. Let us now give two direct applications of matryoshka arithmetic. Assume that we are given an SLP > that computes a function >:\\\> and denote by >:\\\> the function obtained by evaluating > using approximate floating point arithmetic over>. Matryoshki yield an efficient way to statically bound the error >-f> for inside some fixed poly-ball, as follows: <\proposition> Let =,\,\|)>\\,\|)>> be a fixed poly-ball and ,0|)>\,0|)>>,\,\,0|)>|)>\\,\,\|)>>. Let :,\,\|)>>\,\,\|)>>> be apolymatryoshka lift of and ,E|)>\,E|)>,\,\,E|)>|)>\,0|)>|)>>>. Then for all \> with \\,\,a\\> and ,n>, we have <\equation*> ,i>-f|\|>\E. <\proof> Let :\,\|)>\\,\|)>> be a ball lift for which > satisfies the matryoshka lift property. Since \a> by assumption, we have <\equation*> \,0|)>|)>=\,E|)>\\|)>=\>,r|)>, for some \>. The ball lift property also yields <\equation*> \>,r|)>\f, so ,i>-f|\|>\r\E> for ,n>. As our second application, we wish to construct a ball lift of that is almost as efficient to compute as itself, provided that the input balls are included into fixed large balls ,\,\\\,\|)>> as above. For this, we first compute bounds \\>> for the Jacobian matrix of : <\equation> sup\a> f|\ x>|\|>\B. For this, it suffices to evaluate a ball lift of the Jacobian of at >, which yields a matrix \\,\|)>m>>, after which we take \|\|>>. We recall that an SLP for the computation of the Jacobian can be constructed using the algorithm by Baur and Strassen. <\proposition> Let |log k|\>> for all integers 1> and \\>>. With the above notation, assume that \\>, and let be as in Proposition. For every =>\,\|)>>> with \\,\,\\\>, let <\equation*> \>|)>\\>,>**\|)>+3*\|]>|)>. Then >> is a ball lift of . (Note that this lift is only defined on the set of \\,\|)>> such that \\,\,\\\>.) <\proof> Let \\,\,b\\>. Then, for ,n>, we have <\equation*> -f|\|>\B*r+\+B*r, by() and the mean value theorem. In other words, <\equation*> \,B*r|)>\f. Applying Proposition, it follows that <\equation*> \>,E+B*r|)>\f. By computing the sum *r+\+B*r> via a balanced binary tree and using, weobtain <\eqnarray*> +B*r+\+B*r>|>|>+B*r+\+B*r|)>|]>*|)>+\>|)>*|)>>>||>|>+B*r+\+B*r|)>|]>*|)>+2*\>>||>|>+B*r+\+B*r|]>**\|)>+3*\|)>*|)>>>||>|>+B*r+\+B*r|)>**\|)>+3*\|]>,>>>> where we also used Lemma and the fact that the computation of is exact. It is classical to extend ball arithmetic to other operations like division, exponentiation, logarithm, trigonometric functions, etc. In this section, we will focus on the particular case of division. It will actually suffice to study reciprocals \x>, since >. Divisions and reciprocals raise the new difficulty of divisions by zero. We recall that the IEEE 754 provides a special not-a-number element NaN in > that corresponds to undefined results. All arithmetic operations on NaN return NaN, so we may check whether some error occurred during a computation, simply by checking whether one of the return values is NaN. For exact arithmetic, we will assume that > is extended in a similar way with an element such that \NaN>. In the remainder of this section, assume that =\> or =\>. In order to extend reciprocals to balls, it is convenient to introduce the function <\equation*> |\>\,>||x\0>>||| x\0>>>>> Then the exact reciprocal of a ball \\,\|)>> can be defined as follows: <\eqnarray*> |)>>|>|,r*|\>-r|)>*|)>|)>,>>>> where we note that <\equation*> r*\-r|)>*|)>=-r|)>*>=-r>->. If \\,\|)>> then we may also use the following formula <\eqnarray*> ||>|)>>|>|>,>|\>-r|)>*|)>|]>|)>>>>> for transient ball arithmetic and <\eqnarray*> |>>|)>|)>>|>|>,r*|\>-r|)>*|)>+|\>>|)>|)>.>>>> for semi-exact ball arithmetic. We assume that the rounded counterpart of the reciprocal verifies the condition <\equation> >-\|\|>\>|\|>*\>, for all invertible in >. For IEEE 754 arithmetic, this condition is naturally satisfied for>\2> when =\>, and in absence of overflows and underflows. If =\> and \>, then we compute reciprocals using the formula <\eqnarray*> >*b|)>>|>|>*b|)>*\+b|)>|]>.>>>> For this definition, we have <\eqnarray*> >*b|)>-\*b|)>|\|>>|>|*b|\|>*+b|)>|]>-\+b|)>|\|>+>*b|)>|\|>*\>.>>>> Now <\equation*> ||+b>|>|>a+b\>b|)>*>|)>>|>|>+b|]>*>|)>>>|+b>|>|>a+b\>b|)>*>|)>>|>|>+b|]>*>|)>,>>>>> whence <\equation*> +b|)>>|>|>+b|]>|)>*>|)>>|>|>+b|)>|]>*>|)>*>|)>>>|+b|)>>|>|>+b|]>|)>*>|)>>|>|>+b|)>|]>*>|)>*>|)>.>>>>> For >\>, it follows that <\eqnarray*> *b|\|>*+b|)>|]>-\+b|)>|\|>>|>|*b|\|>*+b|)>|]>*>+6*\>|)>|)>>>||>|*b|)>*\+b|)>|]>|\|>*>+10*\>|)>>>|||>*b|)>|\|>*>+10*\>|)>>>||>|>*b|)>|\|>*>|)>.>>>> Consequently, <\eqnarray*> >*b|)>-\*b|)>|\|>>|>|>*b|)>|\|>*>|)>.>>>> This shows that we may take >\5*\>>. Let us denote \max>,\>|)>> and assume that \>. We are now in a position to prove a counterpart of Lemma for division. <\lemma> Let \\,\|)>> and \\>> be such that -r>\\>. Assume that the computation of <\eqnarray*> |>|-r|)>*|)>|]>>>>> involves no overflows and no underflows. Then -r|)>*|)>\R*|)>+7>>. <\proof> We have <\eqnarray*> >-r>|>|*|)>-r>>|||-r+r*\|)>*|)>>>||>|-r|)>**\|)>*|)>,>>>> whence <\eqnarray*> -r|)>*>|>|>-r|)>**|)>**\|)>>>||>|>-r|)>*>*|)>**\|)>>>||>|>->r|)>*>*|)>**\|)>>>||>|>->r|)>\>>*|)>**\|)>,>>>> whence <\eqnarray*> -r|)>*|)>>|>|>->r|)>\>>|)>*|)>**\|)>>>||>|>>->r|)>\>>|)>*|)>**\|)>*|)>>>||>|>-r|)>*|)>|]>*|)>**\|)>*|)>.>>>> We conclude by observing that |)>**\|)>*|)>\|)>+7>>, since \>. We adapt the definition of the remaining path length to our extended notion of SLPs with reciprocals. We now take =max >,\,q>|)>+1>, where \\\i> are those indices k> such that > is of the form \Y\Z>> with \,Z|}>>, >\|}>>>, and \,\,X|}>>, or > is of the form \\|)>> and \,\,X|}>>. Lemma allows us to extend Theorem as follows. <\theorem> With the setup from Theorem, assume that \|)>+7|)>*q>>> and that>may also contain reciprocals. During the second evaluation, assume that no underflow or overflow occurs and that -r>\\>> for every computation of a reciprocal of >>, where \\>>. Assume in addition that *q>|)>\\>> and that the following two inequalities hold: <\equation*> \\H>>|)>+7|)>*q>>*|1+\>*|)>+7|)>*q>>|1+\>|)>,\\max+9|2>,|\>*\|)>. Then >\t> for any corresponding outputs >|)>> and> for the two evaluations. <\proof> We adapt the proof of5>, by showing that the two inductions still go through in the case of reciprocals. For the first induction, consider an instruction \\|)>>, where =\> and -r>\\>. Recall that <\equation*> r\*|)>*+1|)>>-1|)>. Using \4>, *q>|)>\\>, and Lemma, we note that <\eqnarray*> |)>-2>**\|)>>|>|-1|)>*\|)>**\|)>>>|||-\*-1|)>*\>>||>|.>>>> Combined with Lemmas and, we deduce that <\eqnarray*> >||-r|)>*|)>|]>>>||>|*-r|)>>*|)>+7|)>>>>||||\|>*/r-1>*|)>+7|)>>>>||>||\|>*|)>*+1|)>>-1|)>|2-|)>*+1|)>>>*|)>+7|)>>>>||>|>|\|>*|)>*|)>*+1|)>>-1|)>|2-*\|)>>*|)>+7|)>>>>||>|>|\|>*|)>*|)>*q>-1|)>*|)>>|2-*\|)>>*|)>+7|)>>>>||>|>|\|>*|)>*q>-1|)>*|)>-+9|)>>.>>||>|>|\|>*|)>*q>-1|)>.>>>> For the second induction, we redefine <\equation*> \\>+|\*|)>>*H,q>>|)>*|)>+7|)>*>-q|)>>. Assume \\|)>>, let >|)>> (resp. >) be the value of > in the semi-exact evaluation (resp. transient evaluation), and let be the biggest index so that =X>. Then <\eqnarray*> >>||>*\-r>|)>*|)>+|\>>|)>>>||>|>*\-r>|)>*|)>+\*t>>||>|*r*\-r|)>*|)>+\*t>>||>|*|)>+7>*t+\*t>>||>|*t.>>>> The two inductions still hold for the other operations since we increased > and >. If the conditions of the theorem are all satisfied for a given input poly-ball >, then we note that the numeric evaluation of the SLP at any point > with \> never results in adivision by zero. Conversely however, even when all these numerical evaluations are legit, it may happen that the condition that -r>\\> is violated for the computation of some reciprocal |)>>. This can either be due to an overly optimistic choice of > or to the classical phenomenon of overestimation. Intuitively speaking, a low choice of > means that balls that need to be inverted should neatly steer away from zero, even when inflating the radius by a small constant factor. Infact, this is often a reasonable requirement, in which case one may simply take \1>>. A somewhat larger choice like \10> remains appropriate and has the advantage of increasing the resilience of reciprocal computations (the condition -r>\\> being easier to satisfy). Large values of > deteriorate the quality of the obtained bounds for adubious further gain in resilience. The formulas(), (), and() for reciprocals can again be specialized to matryoshki modulo the precaution that we should interpret > as a lower bound for the norm. For \\> and >\\>|)>>, this leads to the definitions <\eqnarray*> |||,r|)>|)>>|>||)>,r*|\>|\|\>-r|)>*|\|\>|)>|)>>>|>,r|)>|)>>|>|>|)>,>|\>|\|\>-r|)>*|\|\>|)>|]>|)>>>|>>,r>|)>|)>>|>|>>|)>,r>*|\>|\>|\>-r>|)>*|\>|\>|)>+|\>>>|)>|)>|)>,>>>> where we use the lower bound notation |\|\>\-R>. In addition to(), we will assume that the rounded counterpart of the reciprocals also verifies the condition <\eqnarray*> |)>-\>|)>|\|>>|>|>|)>|\|>*\,\|)>>>>>> for all invertible balls > in ,\|)>>. For what follows, let \\,\|)>>>. <\lemma> Let ,\\\>> be such that \1-*\|)>*|)>\>. Let ,r|)>=>\,\,\|)>>> be such that |\|\>>\\> and |\|\>-r>\\>. If the computation of <\eqnarray*> |>||\|\>-r|)>*|\|\>|)>|]>>>>> involves no overflows and no underflows, then |\|\>-r|)>*|\|\>|)>\R*|)>+3|)>*+2|)>+5>>. <\proof> From <\eqnarray*> |\|\>>>|>|>->R>>||>|*|)>-R|)>*|)>>>|||-R+R*\|)>*|)>>>||>|-R|)>**\|)>*|)>>>||||\|\>*|)>>>>> and <\eqnarray*> |\|\>>-r>|>||\|\>*|)>-r>>||||\|\>-r+r*\|)>*|)>>>||>||\|\>-r|)>**\|)>*|)>,>>>> we obtain <\eqnarray*> |\|\>-r|)>*|\|\>>|>||\|\>>-r|)>*|\|\>>**\|)>*|)>>>||>||\|\>>->r|)>*|\|\>>**\|)>*|)>*|)>>>||>||\|\>>->r|)>\>|\|\>>**\|)>*|)>*|)>.>>>> It follows that <\eqnarray*> |\|\>-r|)>*|\|\>|)>>|>||\|\>>->r|)>\>|\|\>>|)>**\|)>*|)>*|)>>>||>|>|\|\>>->r|)>\>|\|\>>|)>**\|)>*|)>*|)>*|)>>>||>|>|\|\>-r|)>*|\|\>|)>|]>**\|)>*|)>*|)>*|)>.>>>> Using \> and \>, we have <\equation*> |)>=*\|)>*|)>\|)>+3> and then <\equation*> *\|)>*|)>*|)>*|)>\|)>+2|)>>*|)>\|)>+3|)>*+2|)>+5>. We are now ready to extend Theorem to SLPs with reciprocals. <\theorem> Let ,\\\>> be such that \1-*\|)>*|)>\> and let +3|)>*+2|)>>+5>. With the setup from Theorem, assume that \|)>>>>> and that > may also contain reciprocals. During the second evaluation, assume that no underflow or overflow occurs and that |\|\>>\\>> and |\|\>-r>\\> hold for every computation of a reciprocal of ,r|)>=\>>. Assume in addition that *q>|)>\\>> and that the following two inequalities hold: <\equation*> ||>|>|>>|)>>>*|1+\>*|)>>>|1+\>|)>>>|\>|>|,|\>*\|)>.>>>>> Then >\t> and >\T> for any corresponding outputs >,t>|)>> and> for the first and second evaluations. <\proof> We adapt the proof of Theorem by showing that the two inductions still go through in the case of reciprocals. For the first induction, consider an instruction \\|)>>, where =\=\,r|)>>, with |\|\>>\\> and |\|\>-r>\\>. Recall that <\eqnarray*> |>||\|>*|)>*+1|)>>-1|)>.>>>> Using Lemma, the identity |)>|\|>=|\|\>>>, as well as and, we obtain <\eqnarray*> >|||\|\>-r|)>*|\|\>|)>|]>>>||>||\|\>-r|)>*|\|\>>*|)>>>||||)>|\|>*|\|\>-r>*|)>>>||||)>|\|>*|\|\>/r-1>*|)>>>||>||)>|\|>*|\|>/r-1>*|)>>>||>||)>|\|>*|)>*+1|)>>-1|2-|)>*+1|)>>-1>*|)>>>||>|>|)>|\|>*|)>*|)>*+1|)>>-1|1-\*\>>*|)>>>||>|>|)>|\|>*|)>*+1|)>>-1|)>*|)>->>>||>|>|)>|\|>*|)>*q>-1|)>*|)>>*|)>->>>||>|>|)>|\|>*|)>*q>-1|)>.>>>> For the second induction, we redefine <\equation*> \\>+|\*|)>>*H,q>>|)>*|)>>-q|)>>. Assume \\|)>>, let >,r>|)>> (resp. ,r|)>>) be the value of> in the semi-exact evaluation (resp. transient evaluation), and let be the biggest index so that =X>. Then Lemma implies <\eqnarray*> >>||>*\|\|\>-r>|)>*|\|\>|)>+|\>>>|)>|)>>>||>|>*\|\|\>-r>|)>*|\|\>|)>+\*t>>||>|*r*\|\|\>-r|)>*|\|\>|)>+\*t>>||>|*|)>*t+\*t>>||>|*t.>>>> The two inductions still hold for the other operations since we increased > and >. Similar comments as those made after Theorem also apply to the above theorem. In particular, we recommend taking ,\\>. In the previous sections we have focussed on more efficient algorithms for the reliable evaluation of functions on fixed bounded poly-balls. This technique applies to very general SLPs that may involve operations like division, which are only locally defined. Nonetheless, polynomial SLPs, which only involve the ring operations > are important for many applications such as numerical homotopy continuation. In this special case, we will show how remove the locality restriction and allow for the global evaluation of such SLPs in a reliable and efficient way. Consider a polynomial map <\eqnarray*> >|>|>>|,\,x|)>>|>|=,\,x|)>,\,P,\,x|)>|)>>>>> for polynomials ,\,P\\,\,x|]>>. Given ,n|}>>, the polynomial > is not homogeneous, in general, but there exists a homogeneous polynomial \\,\,x|]>> such that ,\,x|)>=,\,x,1|)>>>. This polynomial is unique up to multiplication by powers of >. The polynomials ,\,P> give rise to a polynomial map <\eqnarray*> :\>|>|>>|,\,x|)>>|>|,\,x|)>,\,P,\,x|)>|)>,>>>> which we call a of >. Assume now that the map >> can be computed using an SLP> of length without divisions. Let us show how to construct an SLP > such that \E>> is a homogenization of. Consider the formal evaluation of > over ,\,x|]>>, by applying >> to the input arguments ,\,x>. Then the output value > and the input arguments> and > of an instruction > of the form \Y\Z> (>\|}>>) are polynomials in ,\,x|]>> that we denote by >> >>, and >>, respectively (for instructions \C>, we set >\C>). We can recursively compute upper bounds >,k>>, >,k>>, >,k>> (abbreviated as>>>, >>>, >>>) for the total degrees of these polynomials: <\itemize> If > is of the form \C>, then we take >>\0>. If > is of the form \Y\Z>, then >>\max>>,d>>|)>>. If > is of the form \Y\Z>, then >>\d>>+d>>>. Here >>\1> ( >>\1>) whenever > ( >) is an input variable. Now consider the program |~>> obtained by rewriting each instruction > as follows: <\itemize> If > is of the form \C> or \Y\Z>, then > is rewritten into itself. If > is of the form \Y\Z> with >>=d>>>, then > is rewritten into itself. If > is of the form \Y\Z> with >>\d>>>, then > is rewritten into two instructions \Y\x>>-d>>>> and \A\Z>, for some new auxiliary variable >. If > is of the form \Y\Z> with >>\d>>>, then > is rewritten into two instructions \Z\x>>-d>>>> and \Y\A>, for some new auxiliary variable >. By induction, one verifies that the image of > under this rewriting computes the unique homogenization ,k>> of ,k>> of degree >>> for ,l>. We obtain > by prepending|~>> with instructions that compute all powers > that occur in |~>>. This can for instance be done using binary powering: \x\x> and \x\x> for all 1>. <\example> Applying the homogenization procedure to the left-hand SLP below, we obtain the right-hand one: <\equation*> ||\x\x>>|||\x\x>>|||\x\x>>|||\x\x>>|\x\x>||\x\x>>|\V+x>||\x\x>>|||homogenize>>>|\V+A>>|\V\V>||\V\V>>|\V+x>||\x\x>>|||\V+A>>|\V\V>||\V\V>>|\V+x>||\x\x>>|||\V+A>>>>> In this example, >>=d>>=2>, >>=d>>=4>, and >>=d>>=8>. In the last instruction, we thus have >>=8\d>>=1>, whence the exponent >>-d>>> in \x\x>. <\remark> In the worst case, computing the powers of > using binary powering may give rise to a logarithmic overhead. For instance, there is a straightforward SLP> of length proportional to that computes the polynomials >+1> for ,n>. But when computing using binary powering in order to compute its homogenization, the length of> is proportional to . An alternative way to compute the required powers is to systematically compute >>>> and |)>>>>> for all . In that case >>-d>>>> and >>-d>>>> can always be obtained using a simple multiplication and the length of > remains bounded by > for an SLP > of length . Of course, this requires division to be part of the SLP signature, or > to be passed an argument, in addition to >. Consider an output variable > with ,n|}>> and denote by ,i>:\\\> and ,i>:\\\>> the -th component of >> and >>, respectively. If ,l|}>> is largest with =O>, then we define \d>>>, and we note that ,i>> is the unique homogenization of ,i>> of degree >. Let ,\,d|)>> and define \*t>,\,A*t>|)>> for any ,\,A\\>. For any \> and \>, it then follows that <\equation> E>=E>*t. We will denote by \\,\,x|]>> the polynomials with >=P> for all \>. With the notation from the previous subsection, assume that =\> with =\> or =\>. Suppose that we wish to evaluate at a point \> and let <\equation> t\max|\|>,\,|\|>,1|)>,u\,x\,\,u*x,u|)>. Of course, we may simply evaluate > at , which yields =E>>. But thanks to(), we also have the projective evaluation method <\equation*> P=E>|)>*t. The advantage here is that <\equation> |\|>>\max |\|>,\,|\|>|)>\1, so we only evaluate > at points in the poly-ball >. This allows us to apply the material from the previous sections concerning the evaluations of SLPs at points inside a fixed ball. Note that we may easily add the function to SLPs since its computation is always exact. If =\> then we apply to the real and imaginary parts. For instance, let ,r|)>,\,\,r|)>|)>\\,\|)>> be the evaluation of > at>> using ball arithmetic and set \|\|>+r> for ,n>. Then we have <\equation> |\|>\M*t> for all \> and ,n>. Consequently, we may compute upper bounds for |\|>,\,|\|>>> in time +\+log d|)>>, which is typically much smaller than the length of >. We may use such global bounds for the construction of more efficient ball lifts, although the resulting error bounds may be less sharp. For this, we first evaluate the Jacobian matrix of >> at >, which yields bounds <\equation> E,i>|\ x>|\|>\M, for all ,n>, ,m+1>, and \>. For any ball =,r|)>,\,\,r|)>|)>> with ,r|)>\\> for ,m+1>, we then have <\equation> P|)>\\,i>,M*r+\+M*r|)>. This allows us to compute an enclosure of |)>> using a single evaluation of > at plus > extra operations (the quantity > can be further reduced to > by replacing the > by ,\,M|)>> or ,\,M|)>>). Similarly, if =,r|)>>,\,\,r|)>|)>> is any ball, then we define <\equation> t\max|\|>+r,\,|\|>+r,1|)>,u\,\\,\,u*\,\|)>. For ,n>, the relations() and() now yield <\equation> P|)>\\,i>,*r+\+M*r|)>*t-1>|)>. (If =0>, then ,i>> is a constant, and we may use zero as the radius.) This allows us to compute an enclosure of |)>> using one evaluation of > at plus +\+log d|)>> extra operations in >. Note that the homogenization > was only required in order to compute the constants >, but not for evaluating the right-hand side of() at a specific >. In the previous subsection, we assumed infinitely precise arithmetic over > or >. Let us now show how to adapt this material to the case when =\>, where =\> or=\>>. As usual, assume \max >,\>|)>\>, and let us first show how to deal with rounding errors in the absence of overflows and underflows. We first replace() by <\eqnarray*> >>|>|>|\|>,\,|\|>,1|)>*|)>|]>>>|>>|>||)>/t>|]>>>|>>|>|>*x,\,u>*x,u>|)>|]>>>>> and verify that <\eqnarray*> >>|>|>|\|>,\,|\|>,1|)>|]>*|)>*|)>\t*|)>*|)>\t>>|>>|>|>|\|>,\,|\|>,1|)>|]>*|)>*|)>\t*|)>*|)>\t*|)>>>|>>|>||)>/>*|)>|)>\|)>*u/|)>\u*|)>>>|>>|>||)>/>*|)>|)>\|)>*u/|)>\u*|)>.>>|>|)>|\|>>|>|>*x|\|>*|)>\u*|)>*t*|)>=1,,m|)>.>>>> We have <\eqnarray*> >-x|\|>>>|>|>|u>*x-x|\|>>+>->|u>*x|\|>>>>||>|>|u>-1|\|>*|\|>>+>|\|>>*\>>||>|>|u>-1|\|>+\>>||>|.>>>> Let ,n|}>>. Evaluating > at ,|\>,\|)>>, for the matryoshka \\|)>>, by Proposition for instance, we may compute a \\> with <\equation*> >|)>-P|)>|\|>\\ for all \>, as well as <\equation*> >,i>|]>-P|\|>\\ for all \> with >\1>. It follows that <\eqnarray*> >,i>>|)>|]>*t>>-P|\|>>|>|>|)>*t>>-P|)>*t>|\|>+\*t>>>>||>|*>|t>|)>>-1|)>|\|>+2*\*t>>>>||>||\|>*|)>>-1|)>+2*\*t>>.>>>> Setting <\equation*> P>\>,i>>|)>*t>>|]>, we then get <\eqnarray*> >-P|\|>>|>|>,i>>|)>|]>*t>>-P|\|>+>,i>>|)>|]>*t>>|\|>|)>+1>-1|)>>>||>|>,i>>|)>|]>*t>>-P|\|>*|)>+1>+|\|>|)>+1>-1|)>>>||>||\|>*|)>>-1|)>+2*\*t>>|)>*|)>+1>+|\|>|)>+1>-1|)>>>||>||\|>*|)>+1>-1|)>+2*\*t>>*|)>+1>.>>>> Note that for all 2> and \a/2> we have <\equation*> a*log|)>\log*\|)>, whence |)>-1\a*\>. Consequently, provided that \|)>/2>, the bound simplifies into <\eqnarray*> >-P|\|>>|>||)>*\*>|\|>+3*\*t>>.>>>> This yields an easy way to compute error bound <\eqnarray*> >-P|\|>>|>|>|)>*\*>|\|>+3*\*t>>|)>*|)>+3>|]>>>>> \ for the approximate numeric evaluation of > at any point \>. Let us now turn to the bound(). Using traditional ball arithmetic (formulas) over ,\|)>>>, we may still compute \\> with |\|>\M> for all \> with >\1>>. Then we simply have <\equation*> |\|>\M*t>> for all \>, since t>>. In a similar way, bounds > that satisfy() can be computed using traditional ball arithmetic. For any ball =,r|)>>,\,\,r|)>|)>>\\,\|)>> and with we notation from(), the enclosure() still holds. Combining <\equation*> P|)>\\,i>,*r+\+M*r|)>*max|\|>>,\,|\|>>,1|)>-1>|)>. with(), this yields <\equation*> P|)>\\>,i>|]>,*r+\+M*r|)>*t>-1>+|)>*\*>|\|>+3*\*t>>|)>. Hence <\eqnarray*> |)>>|>|>,i>|]>,\|)>,>>>> where <\eqnarray*> >|>|>*r+\+M*r|)>*t>-1>+|)>*\*>|\|>+3*\*t>>|)>*|)>+10>|]>.>>>> This yields a ball lift for that is almost as efficient to compute as a mere numeric evaluation of >>. Let us finally analyze how to deal with underflows and overflows. Let \\>\\>> be such that() holds. As long as the computation of >> does not underflow, the inequality >-x|\|>>\12*\>> remains valid, even in the case of overflows(provided that \2>). Consequently, the relation() still holds, so it suffices to replace()by <\eqnarray*> >-P|\|>>|>|>|)>*\*>|\|>+3*\*t>>|)>*|)>+3>+\|]>.>>>> This indeed counters the possible underflow for the multiplication of |)>*\> with >|\|>>. Note that the relation trivially holds in the case when the right-hand side overflows, which happens in particular when the computation of >> overflows. Similarly, it suffices to replace() by <\eqnarray*> >|>|>*r+\+M*r|)>*t>-1>+|)>*\*>|\|>+3*\*t>>|)>*|)>+10>+m*\|]>.>>>> If the computation of >> does not overflow, but the computation of >> underflows, then the IEEE 754 norm implies that we must have >\0.35*\>, where > is the largest positive floating point number in > with \\>. Consequently, the computation *t>>> overflows whenever \0>, provided that we compute this product as > times >>>. The adjusted bounds are therefore valid in general. For several applications, it is useful to not only compute bounds for the function itself, but also for some of its iterated derivatives. Let us first assume that \\> is an analytic function defined on an open subset > of>. Given a derivation order \> and aball =\\\>, our goal is to compute a bound for \z> >|\|>>. If is actually an explicit polynomial <\equation*> f=f*z+\+f and , then a first option is to use the crude bound <\eqnarray*> \z> >|\|>>|>|!|i!>*|\|>*R.>>>> Of course, the case where 0> may be reduced to the case where via the change of variables z+a>. Assume now that is given through an SLP, which possibly involves other operations as >, such as division. In that case, the above crude bound typically becomes suboptimal, and does not even apply if is no longer a polynomial. A general alternative method to compute iterated derivatives is to use Cauchy's formula <\equation> f>=*\>*>|>*\ u, where \\> is some circle around \>. Taking \\> to be the circle with center and radius for some 0>, this yields the bound <\eqnarray*> \z> >|\|>>|>||f|)>|\>|r>,>>>> where |\|\>\+s> denotes an upper bound for > where \>. It is an interesting question how to choose . One practical approach is to start with any 0> such that \\>. Next we may compute both |f|)>|\>> and |f|)>|)>|\>> for some other \r> and check which value provided a better bound. If, say :=r/2>, yields a better bound, then we may next try \r/4>. If the bound for is better than the one for , we may continue with \r/|)>> and so on until we find a bound that suits us. For simple explicit functions , it is also instructive to investigate what is the optimal choice for . For instance, if and > with k>, then we have |f|)>|\>\R>, so the bound() reduces to <\eqnarray*> \z> >|\|>>|>||r>.>>>> The right-hand side is minimal when <\equation> r=*R. In general, we may consider the power series expansion at <\equation*> f=f+f*z+f*z+\. For such a power series expansion and a fixed 0> with \\>, let \> be largest such that |\|>*> is maximal. We regard > as the \Pnumerical degree\Q of on the disk >). Then the relation() suggests that we should have k*R/-k|)>> for the optimal value of . For higher dimensional generalizations, consider an analytic map \\> for some open set \\>. Let |\|>> be the standard Euclidean norm on > for any \>. Given \> and >|)>>, we denote by <\equation*> \\\,R|)>\\\\,R|)> the polydisk associated to the poly-ball >. For \> and assuming that \\>>, we denote the operator norm of the -th derivative f> of on > by <\equation*> |D f|\<\|\|\>>>\sup\>|\|>=\=|\|>=1> f|)>,\,u|)>|\|>|)>. Here we recall that <\eqnarray*> f|)>,\,h|)>>||=1>>=1> f|\ z>*\*\ z>>*h>*\*h>.>>>> Given >|)>> and \\\\,R+r|)>\\\\,r|)>> with \\>, we have the following -dimensional generalization of(): <\eqnarray*> +\+k> f|\z>*\*\z>>>|>|!*\*k!|*\|)>>*>|-z|)>+1>*\*-z|)>+1>>*\u*\*\ u,>>>> for any \>. For any \z>, it follows that <\eqnarray*> +\+k> f|\z>*\*\z>>|\|>>|>|!*\*k!*|f|)>|\>|r>*\*r>>,>>>> where |\|\>> stands for an upper bound for > with \z>. Translated into operator norms, this yields <\eqnarray*> |D f|\<\|\|\>>>>|>||f|)>|\>*+\+k=k>,\,k>*!*\*k!|r>*\*r>>>>||||f|)>|\>*+\+k=k>>*\*r>>.>>>> A practical approach to find an which approximately minimizes the right-hand side is similar to what we did in the univariate case: starting from a given , we vary it in adichotomic way until we reach an acceptably good approximation. This time, we need to vary the individual components ,\,r> of in turn, which makes the approximation process roughly times more expensive. Multivariate analogues of() are more complicated and we leave this issue for future work. <\bibliography|bib|tm-plain|> <\bib-list|21> A.Ahlbäck, J.vander HoevenG.Lecerf. JIL: a high performance library for straight-line programs. , 2025. G.AlefeldJ.Herzberger. . Academic Press, New York, 1983. W.BaurV.Strassen. The complexity of partial derivatives. , 22:317\U330, 1983. C.BeltránA.Leykin. Robust certified numerical homotopy tracking. , 13(2):253\U295, 2013. P.Bürgisser, M.ClausenM.A.Shokrollahi. , 315. Springer-Verlag, 1997. T.DuffK.Lee. Certified homotopy tracking using the Krawczyk method. , ISSAC '24, 274\U282. New York, NY, USA, 2024. ACM. A.GuillemotP.Lairez. Validated numerics for algebraic path tracking. , 36\U45. ACM, 2024. J.vander Hoeven. Ball arithmetic. A.Beckmann, C.GaÿnerB.Löwe, , 6Preprint-Reihe Mathematik, 179\U208. Ernst-Moritz-Arndt-Universität Greifswald, February 2010. International Workshop. J.vander Hoeven. Reliable homotopy continuation. , HAL, 2011. . J.vander HoevenG.Lecerf. Evaluating straight-line programs over balls. , 142\U149. 2016. Extended preprint version with Appendix at . J.vander HoevenG.Lecerf. Towards a library for straight-line programs. , HAL, 2025. . J.vander Hoeven, G.Lecerf, B.Mourrain etal. Mathemagix. 2002. . L.Jaulin, M.Kieffer, O.DidritE.Walter. Applied interval analysis. , 2001. F.Johansson. Arb: a C library for ball arithmetic. , 47(3/4):166\U169, 2014. R.B.Kearfott. An interval step control for continuation methods. , 31(3):892\U914, 1994. U.W.Kulisch. Computer arithmetic and validity. Theory, implementations and applications. , (33), 2008. R.E.Moore. . Prentice Hall, Englewood Cliff, 1966. R.E.Moore, R.B.KearfottM.J.Cloud. . SIAM Press, 2009. A.Neumaier. . Cambridge University Press, Cambridge, 1990. S.M.Rump. Fast and parallel interval arithmetic. , 39(3):534\U554, 1999. S.M.Rump. Verification methods: rigorous results using floating-point arithmetic. , 19:287\U449, 2010. <\initial> <\collection> <\attachments> <\collection> <\associate|bib-bibliography> <\db-entry|+q3uxtAeUIg|book|Alefeld1983> <|db-entry> J. > <\db-entry|+q3uxtAeUIi|article|Jaulin2001> <|db-entry> M. O. E. > <\db-entry|+q3uxtAeUIk|article|Kulisch2008> <|db-entry> > <\db-entry|+q3uxtAeUIn|book|Moore1966Book> <|db-entry> > <\db-entry|+q3uxtAeUIq|book|Moore2009Book> <|db-entry> R.B. M. J. > <\db-entry|+29X0JvK0Rz0|book|Neu90> <|db-entry> > <\db-entry|+q3uxtAeUJF|article|Rump1999parallel> <|db-entry> > <\db-entry|+q3uxtAeUIv|article|Rump2010> <|db-entry> > <\db-entry|+29X0JvK0S84|inproceedings|vdH:ball:greifswald> <|db-entry> > C. B. > <\db-entry|+29X0JvK0Rtv|article|Joh14> <|db-entry> > <\db-entry|+tRqnglh3Bg4WFp|book|BuClSh97> <|db-entry> M. M. A. > <\db-entry|+ZEhiZRx1lFqmsD3|techreport|vdH:slp> <|db-entry> G. > > <\db-entry|+2V95Ag9w7Pe|inproceedings|vdH:dagball-plus> <|db-entry> G. > > <\db-entry|+29X0JvK0RjT|article|BeltranLeykin2013> <|db-entry> A. > <\db-entry|+33Z4a4KIs9|inproceedings|DuffLee2024> <|db-entry> K. > <\db-entry|+2V95Ag9w7Pg|inproceedings|GL24> <|db-entry> P. > <\db-entry|+29X0JvK0S8K|techreport|vdH:homotopy> <|db-entry> > > <\db-entry|+29X0JvK0Rvc|article|KX94> <|db-entry> > <\db-entry|+29X0JvK0SBS|misc|vdH:mmx> <|db-entry> G. B. > > <\db-entry|+UCHwtGN1sHSR59v|misc|JIL> <|db-entry> J. van der G. > > <\db-entry|+2V95Ag9w7Ph|article|BS83> <|db-entry> V. > <\references> <\collection> > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > <\auxiliary> <\collection> <\associate|bib> Alefeld1983 Jaulin2001 Kulisch2008 Moore1966Book Moore2009Book Neu90 Rump1999parallel Rump2010 vdH:ball:greifswald Joh14 BuClSh97 vdH:slp vdH:dagball-plus BeltranLeykin2013 DuffLee2024 GL24 vdH:homotopy KX94 vdH:dagball-plus vdH:dagball-plus vdH:dagball-plus vdH:dagball-plus vdH:dagball-plus vdH:dagball-plus vdH:dagball-plus vdH:dagball-plus vdH:mmx JIL vdH:dagball-plus vdH:dagball-plus BS83 vdH:dagball-plus <\associate|figure> |2.1>|> Representation of an embedded ball inside a matryoshka. |> <\associate|toc> |math-font-series||font-shape||1.Introduction> |.>>>>|> |math-font-series||font-shape||2.Different types of ball arithmetic> |.>>>>|> |2.1.IEEE floating point arithmetic and notation |.>>>>|> > |2.2.Ball arithmetic |.>>>>|> > |2.3.Matryoshka arithmetic |.>>>>|> > |2.4.Rounded and transient ball arithmetic |.>>>>|> > |2.5.Transient matryoshka arithmetic |.>>>>|> > |math-font-series||font-shape||3.Evaluating straight line programs> |.>>>>|> |3.1.Straight-line programs |.>>>>|> > |3.2.Transient ball evaluations |.>>>>|> > |3.3.From balls to matryoshki |.>>>>|> > |3.4.Applications |.>>>>|> > |math-font-series||font-shape||4.Division> |.>>>>|> |4.1.Reciprocals of balls |.>>>>|> > |4.2.Transient evaluation |.>>>>|> > |4.3.Reciprocals of matryoshki |.>>>>|> > |math-font-series||font-shape||5.Global projective bounds> |.>>>>|> |5.1.Homogenization of SLPs |.>>>>|> > |5.2.Global bounds through homogenization |.>>>>|> > |5.3.Managing rounding errors |.>>>>|> > |math-font-series||font-shape||6.Bounding derivatives> |.>>>>|> |6.1.The univariate case |.>>>>|> > |6.2.Bounds for derivatives of multivariate functions |.>>>>|> > |math-font-series||font-shape||Bibliography> |.>>>>|> <\links> <\collection> > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > >