Certifying floating-point implementations using Gappa

We have already mentioned that floating-point (FP) numbers do not possess basic properties of real numbers. The following FP code sequence, due to Dekker [8], illustrates this:

This sequence consists only of three operations. The first one computes the FP sum of the two numbers $\mathtt{a}$ and $\mathtt{b}$. The second one would always return $\mathtt{b}$ and the third one $0$, if this FP sum was exact. As the numbers here are FP numbers, however, the sum is often inexact, because of the rounding. In IEEE-754 arithmetic with round-to-nearest, under certain conditions, this algorithm computes in $\mathtt{r}$ the error committed by this first rounding. In other words, it ensures that $\mathtt{r}+\mathtt{s}=\mathtt{a}+\mathtt{b}$ in addition to the fact that $\mathtt{s}$ is the FP number closest to $\mathtt{a}+\mathtt{b}$. The Fast2Sum algorithm provides us with an exact representation of the sum of two FP numbers as a pair of FP numbers, a very useful operation.

This example illustrates an important point, which pervades all of this article: FP numbers may be an approximation of the reals that fails to ensure basic properties of the reals, but they are also a very well-defined set of rational numbers, which have other well-defined properties, upon which it is possible to build mathematical proofs such as the proof of the Fast2Sum algorithm.

Let us come back to the condition under which the Fast2Sum algorithm works. This condition is that the exponent of a is larger than or equal to that of b, which will be true for instance if $|\mathtt{a}|\geq|\mathtt{b}|$. If one is to use this algorithm, one has first to prove that this condition is true. Note that alternatives to the Fast2Sum exist for the case when one is unable to prove this condition. The version by Knuth [9] requires $6$ operations instead of $3$. Here, being able to prove the condition, which is a property on values of the code, will result in better performance.

The proof of the properties of the Fast2Sum sequence (three FP operations) requires several pages [8], and is indeed currently out of reach of the Gappa tool, basically because it can not be reduced to manipulating ranges and errors. This is not a problem, since this algorithm has already been proven using machine checkers [10]. We consider it as a generic building-block of larger floating-point programs, and the focus of our approach is to automate the proof of such larger programs. In the case of the Fast2Sum, this means proving the condition.

Let us now introduce the class of larger FP programs that we have targeted in this work: implementations of elementary functions

Current floating-point implementations of elementary functions [11, 12, 13, 14, 15] have several features that make their proof challenging:

• •

  The code size is too large to be safely proven by hand. In the first versions of the CRlibm project, the complete paper proof of a single function required tens of pages. It is difficult to trust such proofs.

• •

  The code is optimized for performance, making extensive use of floating-point tricks such as the Fast2Sum above. As a consequence, classical tools of real analysis cannot be straightforwardly applied. Very often, considering the same operations on real numbers would be simply meaningless.

• •

  The code is bound to evolve for optimization purpose, because better algorithms may be found, but also because the processor technology evolves. Such changes will impose to rewrite the proof, which is both tedious and error-prone.

• •

  Much of the knowledge required to prove error bounds on the code is implicit or hidden, be it behind the semantics of the programming language (which defines implicit parenthesing, for examples), or in the various approximations made. Therefore, the mere translation of a piece of code into a set of mathematical variables that represent the values manipulated by this code is tedious and error-prone if done by hand.

Fortunately, FP elementary function implementations also have pleasant features that make their proof tractable:

• •

  There is a clear and simple definition of the mathematical object that the floating-point code is supposed to approximate. This will not always be the case of e.g. scientific simulation code.

• •

  The code size is small enough to be tractable, typically less than a hundred floating-point operations.

• •

  The control flow is simple, consisting mostly of straight-line code with a few tests but no loops.

The following elaborates on these features.

The evaluation of any mathematical function entails two main sources of errors.

• •

  Approximation errors (also called methodical errors), such as the error of approximating a function with a polynomial. One may have a mathematical bound for them (given by a Taylor formula for instance), or one may have to compute such a bound using numerics [18], for example if the polynomial has been computed using Remez algorithm.

• •

  Rounding errors, produced by most floating-point operations of the code.

The distinction between both types of errors is sometimes arbitrary: for example, the error due to rounding the polynomial coefficients to floating point numbers is usually included in the approximation error of the polynomial. The same holds for the rounding of table values, which is accounted for more accurately as approximation error than as rounding error. This point is mentioned here because a lack of accuracy in the definition of the various errors involved in a given code may lead to one of them being forgotten.

Efficient code is especially difficult to analyze and prove because of all the techniques and tricks used by expert programmers.

For instance, many floating-point operations are exact, and the experienced developer of floating-point code will try to use them. Examples include multiplication by a power of two, subtraction of numbers of similar magnitude thanks to Sterbenz’ Lemma [19], exact addition and exact multiplication algorithms (returning a double-double), multiplication of a small integer by a floating-point number whose mantissa ends with enough zeroes, etc.

The expert programmer will also do his best to avoid computing more accurately than strictly needed. He will remove from the computation some operations that are expected not to improve the accuracy of the result by much. This can be expressed as an additional approximation. However, it soon becomes difficult to know what is an approximation to what, especially as the computations are re-parenthesized to maximize floating-point accuracy.

To illustrate the resulting code obfuscation, let us introduce the piece of code that will serve as a running example along this article.

Listing 2 is an extract of the code of a sine function in the CRlibm library. These three lines compute the value of an odd polynomial $p(y)=y+s_{3}\times y^{3}+s_{5}\times y^{5}+s_{7}\times y^{7}$ close to the Taylor approximation of the sine (its degree-1 coefficient is equal to $1$). In our algorithm, the reduced argument $y$ is ideally obtained by subtracting to the FP input $x$ an integer multiple of $\pi/256$. As a consequence $y\in[-\pi/512,\pi/512]\ \subset\ [-2^{-7},2^{-7}]$.

However, as $y$ is irrational, the implementation of this range reduction has to return a number more accurate than a double, otherwise there is no hope of achieving an accuracy of the sine that matches floating-point double precision. In our implementation, range reduction therefore returns a double-double $\mathtt{yh}+\mathtt{yl}$.

To minimise the number of operations, Horner rule is used for the polynomial evaluation: $p(y)=y+y^{3}\times(s_{3}+y^{2}\times(s_{5}+y^{2}\times s_{7}))$. For a double-double input $y=\mathtt{yh}+\mathtt{yl}$, the expression to compute is thus $(\mathtt{yh}+\mathtt{yl})+(\mathtt{yh}+\mathtt{yl})^{3}\times(s_{3}+(\mathtt{yh}+\mathtt{yl})^{2}\times(s_{5}+(\mathtt{yh}+\mathtt{yl})^{2}\times s_{7})).$

The actual code uses an approximation to this expression: the computation will be accurate enough if all the Horner steps except the last one are computed in double-precision. Thus, $y_{l}$ will be neglected for these iterations, and coefficients $s_{3}$ to $s_{7}$ will be stored as double-precision numbers noted $\mathtt{s3}$, $\mathtt{s5}$ and $\mathtt{s7}$. The previous expression becomes:

However, if this expression is computed as parenthesized above, it will not be very accurate. Specifically, the floating-point addition $\mathtt{yh}+\mathtt{yl}$ will (by definition of a double-double) return $\mathtt{yh}$, so the information held by $\mathtt{yl}$ will be completely lost. Fortunately, the rest of the Horner evaluation also has much smaller magnitude than $\mathtt{yh}$ (this is deduced from $y\in[-2^{-7},2^{-7}]$, therefore $y^{3}\in[-2^{-21},2^{-21}]$). The following parenthesing is therefore much more accurate:

In this last version of the expression, only the leftmost addition must be accurate: we will use a Fast2Sum (which as we saw is an exact addition of two doubles returning a double-double). The other operations use the native — and therefore fast — double-precision arithmetic. We obtain the code of Listing 2.

To sum up, this code implements the evaluation of a polynomial with many layers of approximation. For instance, variable yh2 approximates $y^{2}$ through the following layers:

• •

  $y$ was approximated by $\mathtt{yh}+\mathtt{yl}$ with the relative accuracy $\epsilon_{\mathrm{argred}}$

• •

  $\mathtt{yh}+\mathtt{yl}$ is approximated by $\mathtt{yh}$ in most of the computation,

• •

  $\mathtt{yh}^{2}$ is approximated by yh2=yh*yh, with a floating-point rounding error.

In addition, the polynomial is an approximation to the sine function, with a relative error bound of $\epsilon_{\mathrm{approx}}$ which is supposed known (how it was obtained it is out of the scope of this paper [18]).

Thus, the difficulty of evaluating a tight bound on an elementary function implementation is to combine all these errors without forgetting any of them, and without using overly pessimistic bounds when combining several sources of errors. The typical trade-off here will be that a tight bound requires considerable more work than a loose bound (and its proof might inspire considerably less confidence). Some readers may get an idea of this trade-off by relating each intermediate value with its error to confidence intervals, and propagating these errors using interval arithmetic. In many cases, a tighter error will be obtained by splitting confidence intervals into several cases, and treating them separately, at the expense of an explosion of the number of cases. This is one of the tasks that Gappa will helpfully automate.

We have not yet explained why a tight error bound is required in order to obtain a correctly rounded implementation. This question is surveyed in [5]. To sum it up, an error bound is needed to guarantee correct rounding, and the tighter the bound, the more efficient the implementation. A related problem is that of proving the behaviour of interval elementary functions [20]. In this case, a bound is required to ensure that the interval returned by the function contains the image of the input interval. A loose bound here means returning a larger interval than possible, and hence useless interval bloat. In both case, the tighter the bound, the better the implementation.

As a summary, proofs written for version of the CRlibm project up to versions 0.8 [21] are typically composed of several pages of paper proof and several pages of supporting Maple for a few lines of code. This provides an excellent documentation and helps maintaining the code, but experience has consistently shown that such proofs are extremely error-prone. Implementing the error computation in Maple was a first step towards the automation of this process, but if it helps avoiding computation mistakes, it does not prevent methodological mistakes. Gappa was designed, among other objectives, in order to fill this void.

There has been other attempts of assisted proofs of elementary functions or similar floating-point code. The pure formal proof approach of Harrison [6, 7, 22] goes deeper than the Gappa approach, as it accounts for approximation errors. However it is accessible only to experts of formal proofs, and fragile in case of a change to the code. The approach of Krämer et al [23, 24] relies on operator overloading and does not provide a formal proof.

Section 4 will give examples of Gappa’s syntax and show that Gappa can be applied to mathematical expressions much more complex than just $x+1$, and in particular to floating-point approximation of elementary functions. This requires describing floating-point arithmetic expressions within Gappa.

Gappa only manipulates expressions on real numbers. In the property $x+1\in[2,3]$, $x$ is just a universally-quantified real number and the operator $+$ is the usual addition on real numbers $\mathbb{R}$. Floating-point arithmetic is expressed through the use of “rounding operators”: functions from $\mathbb{R}$ to $\mathbb{R}$ that associate to a real number $x$ its rounded value $\circ(x)$ in a specific format. These operators are sufficient to express properties of code relying on most floating-point or fixed-point arithmetics.

Verifying that a computed value is close to an ideal value can now be done by computing an enclosure of the error between these two values. For example, the property “$x\in[1,2]\Rightarrow\circ(\circ(2\times x)-1)-(2\times x-1)\in[?,?]$” expresses the absolute error caused during the floating-point computation of the following numerical code:

Infinities and NaNs (Not-a-Numbers) are not an implicit part of this formalism: The rounding operators return a real value and there is no upper bound on the magnitude of the floating-point numbers. This means that NaNs and overflows will not be generated nor propagated as they would in IEEE-754 arithmetic. However, one may use Gappa to prove very useful properties, for instance that overflows, or NaNs due to some division by $0$, cannot occur in a given code: This can be expressed in terms of intervals. What one cannot prove are properties depending on the correct propagation of infinities and NaNs in the code.

Thanks to the inclusion property of interval arithmetic, if $x$ is an element of $[0,3]$ and $y$ an element of $[1,2]$, then $x+y$ is an element of the interval sum $[0,3]+[1,2]=[1,5]$. This technique based on interval evaluation can be applied to any expression on real numbers. That is how Gappa computes the enclosures requested by the user.

Interval arithmetic is not restricted to this role though. Indeed the interval sum $[0,3]+[1,2]=[1,5]$ do not only give bounds on $x+y$, it can also be seen as a proof of $x+y\in[1,5]$. Such a computation can be formally included as an hypothesis of the theorem on the enclosure of the sum of two real numbers. This method is known as computational reflexivity [27] and allows for the proofs to be machine-checkable. That is how the formal proofs generated by Gappa can be checked independently without requiring any human interaction with a proof assistant.

Such “computable” theorems are available for the Coq [28] and HOL Light [29] proof assistants. Previous work [30] on using interval arithmetic for proving numerical theorems has shown that a similar approach can be applied for the PVS [31] proof assistant. As long as a proof checker is able to do basic computations on integers, the theorems Gappa relies on could be provided. As a consequence, the output of Gappa can be targeted to a wide range of formal certification frameworks, if needed.

Enclosures are not the only useful predicates. As intervals are connected subsets of the real numbers, they are not powerful enough to prove some properties on discrete sets like floating-point numbers or integers. So Gappa handles other classes of predicates for an expression $x$:

As with intervals, Gappa can compute with these new predicates. For example, $\mathtt{FIX}(x,e_{x})\land\mathtt{(}y,e_{y})\Rightarrow\mathtt{FIX}(x+y,\min(e_{x},e_{y}))$. These predicates are especially useful to detect real numbers exactly representable in a given format. In particular, Gappa uses them to find rounded operations that can safely be ignored because they do not contribute to the global rounding error. Let us consider the floating-point subtraction of two floating-point numbers $x\in[3.2,3.3]$ and $y\in[1.4,1.8]$. Note that Sterbenz’ Lemma is not sufficient to prove that the subtraction is actually exact, as $\frac{3.3}{1.4}>2$. Gappa is, however, able to automatically prove that $\circ(x-y)$ is equal to $x-y$.

As $x$ and $y$ are floating-point numbers, Gappa first proves that they can be represented with $24$ bits each (assuming single precision arithmetic). As $x$ is bigger than $3.2$, it can then deduce it is a multiple of $2^{-22}$, or $\mathtt{FIX}(x,-22)$. Similarly, it proves $\mathtt{FIX}(y,-23)$. The property $\mathtt{FIX}(x-y,-23)$ comes then naturally. By computing with intervals, Gappa also proves that $|x-y|$ is bounded by $1.9$. A consequence of these last two properties is $\mathtt{FLT}(x-y,24)$: only $24$ bits are needed to represent $x-y$. So $x-y$ is representable by a single-precision floating-point number.

There are also some specialized predicates for enclosures. The following one expresses the range of an expression $u$ with respect to another expression $v$:

This predicate is seemingly equivalent to the enclosure of the relative error $\frac{u-v}{v}$. It allows, however, to simplify proofs, as the error can now be manipulated even when $v$ is potentially zero. For example, the relative rounding error of a floating-point addition vanishes on subnormal numbers (including zero) and is bounded elsewhere, so the following property holds: $\mathtt{REL}(\circ(x+y),x+y,[-2^{-53},2^{-53}])$ when rounding to nearest in double precision.

Because basic interval evaluations do not keep track of correlations between expressions sharing the same terms, some computed ranges may be too wide to be useful. This is especially true when bounding errors. For example, when bounding the absolute error $\circ(a)-b$ between an approximation $\circ(a)$ and an exact value $b$. The simplest option is to first compute the ranges of $\circ(a)$ and $b$ separately and then subtract them. However, first rewriting the expression as $(\circ(a)-a)+(a-b)$ and then bounding $\circ(a)-a$ (a simple rounding error) and $a-b$ separately before adding their ranges usually gives a much tighter result. These rules are inspired by techniques developers usually apply by hand in order to certify their numerical applications.

Gappa includes a database of such rewriting rules (section 3.5 shows how the user may expand this database). The tool applies them automatically, so that it can bound expressions along various evaluation paths. Since any of these resulting intervals encloses the initial expression, their intersection does, too. Gappa keeps track of the paths that lead to the tightest interval intersection and discards the others, so as to reduce the size of the final proof. It may happen that the resulting intersection is empty; it means that there is a contradiction between the hypotheses of the logical property and Gappa will use it to prove all the goals of the logical property.

Once the logical property has been proved, a formal proof is generated by retracing the paths that were followed when computing the ranges.

When Gappa is not able to satisfy the goal of the logical property, this does not necessarily mean that the property is false. It may just mean that Gappa has not enough information about the expressions it tries to bound.

It is possible to help Gappa in these situations by providing rewriting hints. These are rewriting rules similar to those presented above in the case of rounding, but whose usefulness is specific to the problem at hand.

We consider again the following C code, where we have added the constants:

There is a lot of rounding operations in this code, so the first thing to do is to define Gappa rounding operators for the rounding modes used in the program. In our example, we use the following line to define IEEEdouble as a shortcut for IEEE-compliant rounding to the nearest double, which is the mode used in CRlibm.

Then, if the C code is itself sufficiently simple and clean, the translation step only consists in making explicit the rounding operations implicit in the C source code. To start with, the constants s3, s5 and s7 are given as decimal strings, and the C compilers we use convert them to (binary) double-precision FP numbers with round to nearest. We ensure that Gappa works with these same constants as the compiled C code by inserting explicit rounding operations:

Then we have to do the same for all the roundings hidden behind C arithmetic operations. Adding by hand all the rounding operators, however, would be tedious and error-prone, and would make the Gappa syntax so different from the C syntax that it would degrade confidence and maintainability. Besides, one would have to apply without error the rules (well specified by the C99 standard [32]) governing for instance implicit parentheses in a C expression. For these reasons, Gappa has a syntax that instructs it to perform this task automatically. The following Gappa lines

define the same mathematical relation between their right-hand side and left-hand side as the corresponding lines of the C programs. This, of course, is only true under the following conditions:

• •

  all the C variables are double-precision variables,

• •

  the Gappa variables on the right-hand side represent them,

• •

  the compiler/OS/processor combination used to process the C code respects the C99 and IEEE-754 standards and computes in double-precision arithmetic.

Finally, we have to express in Gappa the Fast2Sum algorithm. Where in C it is a macro or function call, for our purpose we prefer to ignore this complexity and simply express in Gappa the resulting behaviour, which is a sum without error (we have here to trust an external proof of this behaviour [8, 9]):

Note that we are interested in the relative error of the sum s with respect to the exact sine, and for this purpose the fact that s has to be represented as a sum of two doubles in C is irrelevant.

More importantly, this adds another condition for this code translation to be faithful: As the proof of the Fast2Sum has as hypothesis that the exponent of yh is larger than or equal to that of r, we now have to prove that. We will simply add this goal to the theorem to prove.

As a summary, for straight-line program segments with mostly double-precision variables, a set of corresponding Gappa definitions can be obtained straightforwardly by just replacing the C = with Gappa IEEEdouble=, a very safe operation.

To analyse this code, we now need a “mathematically ideal” definition of all the variables, a reference with respect to which the error is computed. This notion of mathematically ideal may be quite subtle: What is the mathematically ideal of yh2? It could be

• •

  the exact square of yh (without rounding), or

• •

  the exact square of yh+yl that yh approximates, or

• •

  the exact square of the ideal reduced argument $y$ , which is usually irrational.

The right choice depends on the properties to prove. Here, the mathematically ideal value for both yh+yl and yh will be the ideal reduced argument, which we note My. Similarly, the purest mathematical value that yh2 approximates is noted My2 and will be defined as My2 = My*My.

For the polynomial approximation ts, we could choose, as mathematical ideal, either the value of the function that the polynomial approximates, or the value of the same polynomial, but computed on My and without rounding error. Here we chose the latter, although it is the lesser ideal.

Here come a few naming conventions. The first was obviously that the Gappa variables that mimick C variables have the same name. We also impose the convention that such variables begin with a lowercase letter. In addition, Gappa variables for mathematically ideal terms will begin with a “M”. The other intermediate Gappa variables should begin with capital letters to distinguish them from variables mimicking the code. Of course, related variables should have related and, wherever possible, explicit names. Again, these are conventions and are part of a proof style, not part of Gappa syntax: the capitalization will give no information to the tool, and neither will the fact that variables have related names.

For instance, it will be convenient to define a variable equal to yh+yl:

Defining mathematically ideal values resumes to defining in Gappa what the C code is supposed to implement. For instance, using our previous conventions, the line for ts was probably evaluating the value of the same polynomial of the ideal My:

We have kept the polynomial coefficients in lower case: As already discussed in Section 2, the polynomial thus defined nevertheless belongs to the set of polynomial with real coefficients, and we have means to compute (outside Gappa) a bound of its relative error with respect to the function it approximates.

The link between ts and the polynomial approximating the sine is also best expressed using mathematically ideal values:

To sum up, PolySinY is the actual polynomial with the same coefficients s3 to s7 as in the C code, but evaluated without rounding error, and evaluated on the ideal value that yh+yl approximates.

Another approach could be to use a Taylor polynomial, in which case the approximation error would be given by the rest in the Taylor formula, the ideal polynomial would be the Taylor one, it would have ideal Taylor coefficients (beginning with M), some of which would have to be rounded to FP numbers to appear in the program (lowercase). Gappa could handle it, too, but it would be less convenient.

Another crucial question is, how do we define the real, ideal, mathematical function which we eventually approximate? Gappa has no builtin sine or logarithm. The current approach can be described in English as: “$\sin(\mathtt{My})$ is a value which, as long as $\mathtt{My}$ is smaller than 6.29e-3, remains within a relative distance of 2.26e-24 of our ideal polynomial”. In Gappa, this will translate to some hypotheses in the property to prove:

Here the interval of Yhl is defined by the range reduction, and its bound has been computed separately: it is $\pi/512$, plus some margin that accounts for the inexactness of argument reduction.

Similarly, the relative distance between the sine and the polynomial on this interval is computed outside Gappa. We used to use Maple’s infinite norm, but it only returns an approximation, so this was in principle a weakness of the proof. We now use a safer, interval-based approach [18], implemented in the Sollya tool²²2http://sollya.gforge.inria.fr/.

Of course this tool provides a theorem which could be expressed as
|Yhl| in [0, 6.29e-03] -> |(PolySinY - SinY)/SinY| <= 2.26e-24. We just inject the property of this theorem as an hypothesis in Gappa.

Concerning style, it will be more convenient to have a variable defined as this approximation error. Similarly, it will make the proof much clearer and concise to add, from the beginning, as many definitions as possible for the various terms and errors involved in the computation:

With the previous introduction, the theorem to prove, expressed as implications using classical first-order logic, is stated as follows:

The full initial Gappa script is given below. It adds a more accurate definition of yh and yl, stating that they are double-precision numbers and that they form a disjoint double-double. It also adds a lower bound on the absolute value of the reduced argument, obtained thanks to Kahan/Douglas algorithm. This lower bound is important because it will enable most values to stay away from zero, which ensures that relative errors are not arbitrarily big due to underflow.

Invoking Gappa on this file produces the following output:

Warning: no path was found for Epstotal.
Warning: no path was found for |r / yh|.

Results for |yl / yh| in [0, 1.11022e-16] and |Yhl| in [6.22302e-61, 0.00629]
        and |Epsargred| in [0, 2.53e-23] and |Epsapprox| in [0, 2.26e-24]:
Warning: some enclosures were not satisfied.

This means that Gappa needs some help, in the form of hints. Where to start? There are several way to interact with the tool to understand where it fails.

• •

  We may add additional goals to obtain enclosures for intermediate variables. For instance, adding the goal |My| in ?, we obtain the following answer

  |My| in [0, 0.00629]
  

  Gappa was able to deduce this enclosure from the enclosure of Yhl (hypothesis) and the definition of Epsargred. Similarly, we may check for instance that the built-in engine is able to build a good enclosure of PolySinY, but not of s.

• •

  We may add additional hypotheses and see what progress they entail. For instance, providing a dummy Epsround as an hypothesis allows Gappa to complete the proof, thanks to its automatic hints.

This way it is possible to track the point where Gappa’s engine gets lost, and provide hints to help it.

In our case, the best thing to do is to express all the approximation layers detailed in Section 2.5. Written as Gappa equations, we get:

Remark again that all these relative errors are defined relatively to the most accurate term.

We may add goals for these new relative errors: Gappa will be unable to bound any of them. We have to provide hints.

Consider only Eps4, the relative difference between S3 and PolySinY – which we saw Gappa is able to bound. Both S3 and PolySinY are polynomial expressions without any rounding, and with identical coefficients. Therefore, the difference between them resumes to the difference between Yhl=yh+yl, used in S3, and My, used in PolySinY. We precisely have a measure of this difference: it is Epsargred. The hint we have to provide to Gappa should therefore express Eps4 as a function of Epsargred which, when evaluated by intervals, will provide a tight enclosure. Here is a generic technique to obtain such an hint. We start with Eps4 -> (S3-PolySinY)/PolySinY, which is just the definition of Eps4, and we rewrite it incrementally until we have obtained an expression involving Epsargred. In the following hint, we have left, for the purpose of this tutorial, the intermediate rewriting steps commented out.

Considering the orders of magnitudes (see Section 2.5), a naive interval evaluation of this last expression will be very accurate. Indeed, Mts as well as ts are very small compared to 1, therefore the first term is close to Epsargred. The second is the relative error of ts with respect to Mts (expected to be no larger than $2^{-52}$), multiplied by Mts which is smaller than $2^{-14}$. We therefore have a sum of two small terms which should provide a small enclosure.

Still, even with this hint, Gappa still fails to provide an enclosure for Eps4. Adding the goals Mts in ? and ts in ?, we observe that Mts is properly enclosed, but not ts. We therefore add a definition of the relative error of ts with respect to Mts:

and we perform the same analysis: we describe the succession of approximation layers between ts and Mts, define intermediate error terms for them, and provide hints for bounding them. This will in turn require to explain to Gappa how to go from My to yh2 through three approximation layers.

We will not detail this process line by line. The final Gappa script is given in appendix, and is available from the distribution of CRlibm³³3http://lipforge.ens-lyon.fr/www/crlibm/.

Writing hints is the most time-consuming part of the proof, because it is the part where the designer’s intelligence is required. However, we hope to have shown that it may be done very incrementally.

The example chosen in this article is actually quite complex: its Gappa proof consists of more than 150 lines, half of which are hints. The bound found on Epstotal is $2^{-67,24}$ and is obtained in a few seconds on a recent machine (the time can be longer when there is a dichotomy). The resulting Coq proof is more than 7000 lines long.

Some functions are simpler. We could write the proof of a logarithm implementation [5] with a few hints only [33]. One reason is that the logarithm never comes close to 0, so the full proof can be handled only with absolute errors, for which writing hints is much lighter.