Florent Bréhard

Assia Mahboubi

Damien Pous

A certificate-based approach to formally verified approximations

Keywords: 2012 ACM Subject Classification Theory of computation → Logic approximation theory, Chebyshev polynomials, Banach fixed-point theorem, interval arithmetic, Coq

We present a library to verify rigorous approximations of univariate functions on real numbers, with the Coq proof assistant. Based on interval arithmetic, this library also implements a technique of validation a posteriori based on the Banach fixed-point theorem. We illustrate this technique on the case of operations of division and square root. This library features a collection of abstract structures that organise the specfication of rigorous approximations, and modularise the related proofs. Finally, we provide an implementation of verified Chebyshev approximations, and we discuss a few examples of computations.

Introduction

While numerical analysis offers sophisticated computational methods to solve various function space problems, the numerical errors caused by floating-point computations, discretisations or finite iterations, are a major concern in domains like safety-critical engineering or computer assisted proofs in mathematics. To address these issues, rigorous numerics [START_REF] Tucker | Validated numerics: a short introduction to rigorous computations[END_REF] provides algorithms to compute validated enclosures of the exact solution. However, their correctness is ensured by pen-and-paper mathematical proofs. In particular, there is no guarantee concerning their implementations.

In this regard, formal proof offers the highest level of confidence. Several noteworthy works use formally proved rigorous numerics to completely formalise highly nontrivial mathematical results, like the Flyspeck project [START_REF] Thomas | A formal proof of the Kepler conjecture[END_REF] for the Kepler conjecture or the formal verification [START_REF] Immler | A verified ODE solver and the Lorenz attractor[END_REF] of the computer-aided proof of the Lorenz attractor [START_REF] Tucker | A rigorous ODE solver and Smale's 14th problem[END_REF]. However, those methods often require intensive computations, which rapidly becomes restrictive inside proof assistants. In the context of formal verification, certificate-based methods is an appealing strategy [START_REF] Harrison | A skeptic's approach to combining HOL and maple[END_REF][START_REF] Caprotti | Formal and efficient primality proofs by use of computer algebra oracles[END_REF][START_REF] Barendregt | Autarkic computations in formal proofs[END_REF]. It consists in discharging part of the computation work load to external oracles, while correctness remains guaranteed via a posteriori validation steps performed inside the proof assistant. This approach has mostly been used for the purpose of verifying symbolic computations, e.g. primality proofs [START_REF] Grégoire | A purely functional library for modular arithmetic and its application to certifying large prime numbers[END_REF], but we illustrate here how it can also by used in the context of rigorous numerical analysis.

Interval arithmetic. Invented in the 60s by Moore [START_REF] Moore | Interval Analysis[END_REF] and significantly developed in the 80s by Kulisch et al., interval arithmetic is an essential building block of rigorous numerics.

The key idea consists in using real intervals with representable endpoints (e.g., floatingpoint numbers) as rigorous enclosures of real numbers, and providing operations preserving correctness. For example, from π ∈ [3.1415, 3.1416] and e ∈ [2.7182, 2.7183], one obtains π + e ∈ [3.1415, 3.1416] ⊕ [2.7182, 2.7183] = [5.8597, 5.8599]. Efficient implementations are available, as MPFI [START_REF] Revol | Motivations for an arbitrary precision interval arithmetic and the mpfi library[END_REF], IntLab [START_REF] Siegfried | Intlab-interval laboratory[END_REF], C-XSC [START_REF] Klatte | C-XSC: A C++ class library for extended scientific computing[END_REF], ARB [START_REF] Johansson | Arb: efficient arbitrary-precision midpoint-radius interval arithmetic[END_REF]. The CoqInterval library [START_REF] Melquiond | Proving bounds on real-valued functions with computations[END_REF] moreover provides a fully verified implementation inside the Coq proof assistant.

Rigorous Chebyshev approximations. Interval arithmetic is however not a panacea, and replacing all operations on real numbers by interval ones should always be considered with caution: the dependency phenomenon may lead to disastrous over-approximations. In such cases, higher order methods such as rigorous polynomial approximations (RPAs) are preferable. A pioneer work is that of Berz and Makino on Taylor models [START_REF] Berz | Verified integration of odes and flows using differential algebraic methods on high-order Taylor models[END_REF]. Those provide not only a polynomial, but also a remainder s.t. the latter contains the difference between the former and the represented function. Since then, efforts were made to clarify the definition of RPAs and extend them to other bases, in particular the Chebyshev basis [START_REF] Brisebarre | Chebyshev interpolation polynomial-based tools for rigorous computing[END_REF][START_REF] Joldeş | Rigorous Polynomial Approximations and Applications[END_REF], due to their far better approximation properties than Taylor expansions [START_REF] Nicholas | Approximation Theory and Approximation Practice[END_REF].

On the formal proof side, the CoqInterval library includes an implementation of Taylor models called CoqApprox [START_REF] Martin | Proving tight bounds on univariate expressions with elementary functions in coq[END_REF], allowing in particular for an automated rigorous evaluation procedure of definite integrals inside Coq [START_REF] Mahboubi | Formally verified approximations of definite integrals[END_REF]. Unfortunately, an equally accomplished equivalent with Chebyshev approximations does not exist now. Our contribution is a first step towards a formally proved counterpart of the popular Chebfun package [START_REF] Driscoll | [END_REF] for Matlab.

Fixed-point based a posteriori validation. Some operations in function spaces admit straightforward self-validating algorithms by replacing all operations in R by interval ones.

Unfortunately, more complicated operations (e.g., division, square root, differential equations) face several obstructions: the intervals may fail to give sufficiently tight enclosures, bounds for the remainders may be unknown, or only asymptotic, or depend on noneffective quantities.

In such cases, a posteriori validation techniques are an attractive alternative, widely used in rigorous numerics. They consist in reconstructing afterwards an error bound for a candidate approximation. Dating back from the works of Kantorovich about Newton's method, they gained prominence with the rise of modern computers and were applied to numerous functional analysis problems [START_REF] Edgar | Self-validating numerics for function space problems: Computation with guarantees for differential and integral equations[END_REF][START_REF] Yamamoto | A numerical verification method for solutions of boundary value problems with local uniqueness by Banach's fixed-point theorem[END_REF][START_REF] Van | Chaotic braided solutions via rigorous numerics: Chaos in the Swift-Hohenberg equation[END_REF][START_REF] Lessard | Rigorous numerics for nonlinear differential equations using Chebyshev series[END_REF]. Even more recently, those methods were used to compute RPAs for solutions of linear ODEs [START_REF] Benoit | Rigorous uniform approximation of D-finite functions using Chebyshev expansions[END_REF][START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF]. Broadly speaking, the function of interest is characterised as a fixed-point of a contracting operator, from which an error bound is recovered thanks to the Banach fixed-point theorem [START_REF] Vasile Berinde | Iterative approximation of fixed points[END_REF]Thm. 2.1]. Such techniques are of special interest for formal verification, for they allow one to rely on efficient but untrusted external tools while keeping the trusted codebase small: it suffices to formalise the theory about contracting operators and provide means of computing with those operators.

Contributions and outline.

We present a Coq library that makes it possible to compute rigorous Chebyshev approximations of functions on reals. We support basic operations like multiplication or integration in the standard way. For more complex operations like division and square root, we resort to a posteriori validation techniques, thus making a first step towards a potential cooperation between external numerical tools and Coq.

We use the interval arithmetic provided by CoqInterval, but we design our abstractions for RPAs from scratch: this allows us to experiment with different design choices, with more flexibility. We first describe the main lines of the hierarchy (Section 2): we rely on canonical structures to abstract over the concrete implementation details of interval arithmetic, and we use them to denote both real valued functions and their rigorous approximations. We also abstract away from the concrete basis for approximations, to work in the future with different bases, even non polynomial ones (e.g., Bessel functions). We provide instances for the monomial and Chebyshev bases, the latter being described in Section 3.

The main theorem we need to perform a posteriori validation is the Banach fixed-point theorem, whose formalisation is described in Section 4. We show in Section 5 how to apply this theorem to compute rigorous approximations for division and square root using Newtonlike operators. We finally discuss the benefits of our approach on two examples (Section 6): RPAs for the absolute value function, and verified computation of integrals related to the second part of Hilbert's 16th problem.

Approximating real numbers and functions

Numerical errors come from the estimation of both real numbers, e.g. using floating-point numbers, and real functions, e.g. using polynomials. Rigorous estimations must take all these uncertainties into account. For this purpose, interval arithmetic provides an explicit enclosure and rigorous polynomial approximations attach an interval to a polynomial approximant, which bounds the method error on a given domain. Note that the coefficients of polynomial approximations are usually themselves obtained from evaluations of the function or of its derivatives, and therefore also subject to numerical errors. A formal library about rigorous approximation thus implements several variants of each operation, on real numbers, floats, intervals, mathematical functions, approximants, etc., whose relationships are made precise in the various layers of specifications. Our library features a small hierarchy of structures which formalises and organises the dependencies between these variants.

Reals and Intervals

At the bottom of the hierarchy, structure Ops0 collects the operations available on reals, floats, intervals, but also on polynomials and rigorous approximations. It provides the signature of a ring structure, with symbols +, -, * , 1 and 0 shared by all instances thanks to Coq's system of canonical structures. Yet the ring equational theory is a priori only available for real numbers. These operations are also those trivially self-validating. A super-structure Ops1 collects other operations required on data-structures used for scalars: reals, floats, interval endpoints, intervals, etc. They are not meant to be implemented on polynomial approximations. Structure Rel0 specifies the relationship between the operations of Ops0 on reals and those on intervals. The field rel is a relation between the two instances C and D, which share overloaded notations. The relation will eventually be instantiated with the containment relation between intervals and reals. When doing so, the requirements on the relation precisely correspond to the fact that interval operations properly approximate real operations. A record Rel1 is defined in the very same way for Ops1.

As much as possible, we will work with polymorphic functions like the following one:

Definition f (C: Ops1) (x: C): C := 1 / (1 + sqrt x).
First of all, this allows us to define at once a function on real numbers (here, x → 1 1+ √

x ) and a function on intervals, whatever the implementation of intervals. Second, and even more importantly, the corresponding approximation correctness theorem will always hold-by a parametricity meta-result, such a function f will always satisfy the following lemma:

Lemma rf: ∀ C D (T: Rel1 C D), ∀ x y, T x y → T (f x) (f y).
This is only a meta-result: we need to provide a proof for each function f ; but the proof is always trivial, and we automatise it.

There are however operations which cannot be implemented at this level of abstraction, even if we were to add some operations to the record Ops1. This is typically the case for division and square root of rigorous approximations, which require operations on intervals that do not make sense on real numbers (e.g., computing the range of a function and checking that it is bounded). In order to define those operations while remaining rather agnostic about the choice of interval implementation, we setup an intermediate layer of abstraction using the structure NBH (for neighbourhood):

: ∀ X x, contains X x → ∀ Y y, contains Y y → ∀ z, x≤z≤y →contains (bnd X Y) z; is_ltE: ∀ X Y, wreflect (∀ x y, contains X x → contains Y y → x<y) (is_lt X Y); minE, maxE, botE: ... }.
We will also make use of the two following derived operations:

Definition mag (N: NBH) (X: II): option II := max (abs X). Definition sym (N: NBH) (X: II): II := let X := abs X in bnd (-X) X.

The first one approximates the magnitude as an interval, if possible; the second one returns an interval centered in 0 that contains the argument. Note that we assume that intervals are convex. We provide an instance of this structure using the CoqInterval library, using intervals of floating point numbers from the Flocq library [START_REF] Boldo | Verifying Floating-point Algorithms with the Coq System[END_REF]. It is actually a family of instances indexed by the desired precision.

Abstract functions

The structure FunOps describes inductively the catalogue of expressions that the library can approximate. It is parameterised by a type C of ground values (typically, reals or intervals); it packages a set of basic operations on some abstract type for functions (pointwise addition, multiplication. . . ), together with operations specific to functions: identity, constant function, evaluation, integration. It also asks for division and square root operations; those have an additional argument which is used to pass parameters to the oracles used in the implementation of those operations (for now, the degree of the interpolants). When C is R, the type of real numbers, this structure is instantiated with the standard operations on R→R (ignoring the extra parameters for division and square root); our main goal is to provide instances with intervals for C, with which it is possible to perform computations.

Like for ground values, the structure FunOps makes it possible to write polymorphic functions like:

Definition g (C: Ops1) (F: FunOps C): F := let f: F := div' 33 1 (1 + sqrt' 33 id ) in let a: C := integrate f 0 1 in pi + id * cst a
Such a declaration defines at the same time a function on reals (x

→ π + x 1 0 dt 1+ √ t
) and approximations of it, which will be obvious to prove correct whenever the chosen instance F satisfies appropriate properties. Those instances are obtained using rigorous approximations.

Rigorous Approximations

Approximating a function usually consists in projecting this function onto a finite dimension vector space, by expansion on a basis with appropriate properties. For instance, so-called Taylor models [START_REF] Berz | Verified integration of odes and flows using differential algebraic methods on high-order Taylor models[END_REF], are an instance of rigorous polynomial approximation. They attach an interval bounding the remainder to a certain polynomial, in this case represented in monomial basis, so as to describe a set of functions containing the one to be approximated. More generally in this section, a rigorous approximation refers to a linear combination of elements in a suitable basis, packaged with an interval remainder. In the code, we will also use the shorter term model, by analogy with Taylor models.

A basis is described by a family of functions, non necessarily polynomials, indexed by natural numbers, that is a term T: nat→R→R. The structure BasisOps_on below describes the signature required on a basis T. It is parameterised by the type C of coefficients; sequences of such coefficients (seq C) represent linear combinations of elements of T. Linear operations (+, -, 0) need not be provided since they can be implemented independently from the basis. The range operation is important: its role is to bound the range on the given domain; it should be as accurate as possible since it is used at many places to compute error bounds in rigorous approximations (e.g., for multiplication and a posteriori validation). We define BasisOps to be a polymorphic function so that we capture with a single object the idealised operations on reals and their concrete implementation with intervals. Given such operations, we equip type seq C with the basic operations in Ops0. Then we can define rigorous approximations:

Record Model C := { pol: seq C; rem: C }.
Like with seq C, we equip Model C with the basic operations in Ops0, and then with those from FunOps. For instance, addition, evaluation and integration are defined as follows: This is also the case for division and square root, which we will discuss in Section 5. All in all, we obtain instances FunOps through a construction of the following type: It finally remains to show that those operations defined on rigorous approximations properly match the idealised operations on functions over reals. We fix in the sequel an instance N: NBH of neighbourhood and basis operations B: BasisOps, and we write Model for Model II). The central definition to establish this correspondence is the following one, where the function eval is the obvious evaluation function for linear combinations of elements of T.

Definition mcontains (F: Model

) (f: R → R) := ∃ p: seq R, scontains (pol F) p /\ ∀ x, lo≤x≤hi → contains (rem F) (f x -eval T p x)
Intuitively, a model contains a real-valued function f if it contains a generalised polynomial which is close enough to f on the domain of the basis. (The binary predicate scontains denotes the pointwise extension of the relation contains to sequences: in the definition, the real coefficients of p should be pointwise contained in the interval coefficients of pol F.)

Equipped with this definition, we prove lemmas like:

Lemma rmmul: ∀ F f G g, mcontains F f → mcontains G g → mcontains (F*G) (f*g). Lemma rmdiv: ∀ n F f G g, mcontains F f → mcontains G g → mcontains (div' n F G) (div' n f g). Lemma rmintegrate: ∀ F f A a B b, (∀ x, lo≤x≤hi → continuous_at f x) → lo≤a≤hi → lo≤b≤hi → mcontains F f → contains A a → contains B b → contains (integrate F A B) (integrate f a b).
Of course, we need assumptions on the basis operations in order to do so. Those assumptions are summarised in the following structure. Recall that a B: BasisOps provides us with operations B ROps1 on reals and operations B II on intervals. The structure assumes: 1) the expected properties on the operations on reals, i.e, efficient evaluation corresponds to evaluation with T, multiplication indeed corresponds to pointwise multiplication under evaluation, etc.; and 2) a relationship between the operations on reals and on intervals. This separation of concerns is very convenient: the latter containment lemmas are always proved in a trivial way (i.e., automatically), and the former properties do not involve intervals at all, but only real numbers and functions, for which usual mathematical intuitions apply.

Record ValidBasisOps (N: NBH) (B: BasisOps) := { (* properties of operations on reals (B ROps1) *) lohi: lo < hi; bevalE: ∀ p x, beval p x = eval T p x; eval_cont: ∀ p x, continuity_pt (eval T p) x; eval_mul: ∀ p q x, eval T (bmul p q) x = eval T p x * eval T q x; eval_prim: ∀ p a b, eval T (bprim p) b -eval T (bprim p) a = RInt (eval T p) a b; ... (* relationship between operations on intervals (B II) and on reals (B ROps1) *) rbeval: ∀ P p X x, scontains P p → contains X x → contains (beval P X) (beval p x); rbmul: ∀ P p Q q, scontains P p → scontains Q q → scontains (bmul P Q) (bmul p q); rbprim: ∀ P p, scontains P p → scontains (bprim P) (bprim p); ... }.

Arithmetic on Chebyshev polynomials

In order to use the previously described rigorous approximations, it remains to provide implementation of operations (BasisOps) for certain families T of functions. We provide two instances of them: one for the standard monomial basis, where T n x = x^n, and one described in this section for Chebyshev basis, where T n is the n-th Chebyshev polynomial. Chebyshev polynomials are defined by the following recurrence, which immediately translates to a recursive definition in Coq.

T 0 = 1 T 1 = X T n+2 = 2XT n+1 -T n
We can then prove simple properties of those polynomials, for instance:

T n T m = (T n+m + T m-n )/2 (n ≤ m) ( 1
)
T 0 = T 1 T 1 = T 2 4 T n+3 = T n+3 2(n + 3) - T n+1 2(n + 1)
(2)

T n (cos t) = cos(nt) (3)
Those are proved in a few lines using existing lemmas about derivation and cosine.

Clenshaw's evaluation algorithm

The first operation we must implement for BasisOps is the evaluation function (beval). This operation should be polymorphic and as efficient as possible: it will be executed repeatedly when constructing and using rigorous approximations. We use Horner evaluation scheme for the monomial basis, and Clenshaw's algorithm [START_REF] Fox | Chebyshev polynomials in numerical analysis[END_REF] for Chebyshev, which are both linear in the number of elementary operations. The latter is usually presented as a dynamic programming routine. We fix abstract operations C: Ops1 for the remaining Coq snippets in this section, and we translate this routine into a recursive function with two accumulators: In the right-hand side, catrev is the function that reverses its first argument and catenate it with the second one. The proof is done by induction in just three lines, using the Coq tactic for ring equations [START_REF] Grégoire | Proving equalities in a commutative ring done right in coq[END_REF]. Correctness (i.e., field bevalE from structure ValidBasisOps) follows.

Note that while the definition of beval can be used with any Ops1 structure, its correctness is proved only on reals: the lemma ClenshawR does not hold in every Ops1 structure. The behaviour of beval on those structures is specified only through the fact that it respects containments (field rbeval from structure ValidBasisOps, which is proved automatically.)

Multiplication

Another important operation is multiplication. Again, this operation should be polymorphic, and efficient. A difficulty here is that due to Equation ( 1), the n-th coefficient of a multiplication potentially depends on all coefficients of its arguments, not only on the coefficient of smaller rank. We use the following definition, with two auxiliary recursive functions: (sadd (sscal a q') (sscal b p'))) (mul_mns p' q') end. Definition smul (p q: seq C): seq C := sscal (1/2) (sadd (mul_mns p q) (mul_pls p q)) where sscal is the scalar multiplication for polynomials-we cannot yet use the standard notation for this operation since we are in the process of defining an Ops0 structure on seq C. The function mul_pls actually corresponds to multiplication in the monomial basis, it covers the first summand in the right-hand side of (1). The function mul_mns differs only in the fact that the recursive call is not pushed away using two 'cons' operations; it covers the second summand in the right-hand side of (1). Like previously, that smul preserves containments (field rbmul of structure ValidBasisOps) is obvious: this operation only performs a finite sequence of operations preserving containments. Proving that it behaves correctly on reals numbers is more interesting; the key invariant is the following one:

Lemma eval_mul_: ∀ (p q: seq R) n x, eval_ n p x * eval_ n q x = (eval (mul_mns p q) x + eval_ (n+n) (mul_pls p q) x)/2.

Here, eval_ n p evaluates P padded with n zeros in front of it. Again, the difficulty is to find the lemma: it is proved in six lines using [START_REF] Barendregt | Autarkic computations in formal proofs[END_REF], and correctness of smul on reals immediately follows. Taking primitives in Chebyshev basis follows the same pattern (see Appendix A).

Range

As mentioned above, we need accurate estimations of the range of a given polynomial in order to be able to compute precise rigorous approximations. This range can always be estimated by evaluating the polynomial on the interval representing the domain (i.e., given p: seq C, compute beval p (bnd lo hi)). This technique is however not sufficient in practice: this tends to produce largely over-estimated bounds. With Chebyshev basis we can proceed differently: indeed, thanks to Equation (3), T n ranges over [-1; 1] on [-1; 1]. Therefore, the range of a polynomial on [-1; 1] can be estimated by using the sum of the absolute values of the coefficients in Chebyshev basis (and actually, we do not need to take the absolute value of the first coefficient since T 0 = 1).

Rescaling

Putting everything together, we obtain the polymorphic operations chebyshev.basis: BasisOps, which can readily be used to construct rigorous approximations, with the instance MFunOps from Section 2.3. This basis however requires to work on the domain [-1; 1] (for estimating the range as explained in the previous section, but also to perform interpolation, see Section 5.1).

In order to use it on other domains, we provide a rescaling function that takes a B: BasisOps and rescales it to a given interval [a; b] using the obvious affine function. We show that this operation preserves validity of basis operations, so that we can use it whenever needed.

Formalisation of Banach fixed-point theorem

Banach fixed-point theorem is the cornerstone of the method discussed here.

Theorem 1 (Banach fixed-point). Let (X, • ) be a Banach space, an operator F : X → X, h • ∈ X, and µ, b, r ∈ R + , satisfying the following conditions:

(i) h • -F • h • ≤ b; (ii) F is µ-Lipschitz over the closed ball B(h • , r) := {h ∈ X | h -h • ≤ r}: ∀h 1 , h 2 ∈ X, h 1 ∈ B(h • , r) ∧ h 2 ∈ B(h • , r) ⇒ F • h 1 -F • h 2 ≤ µ h 1 -h 2 ; (iii) µ < 1 : F is contracting over B(h • , r); (iv) b + µr ≤ r.
Then F admits a unique fixed-point h * in B(h • , r).

This classic result has been formalised in various flavours of logic and proof assistants. In particular, Boldo et al. have provided a formal proof of a version of this fixed-point theorem, based on the Coquelicot library, for the purpose of the formalisation of the Lax-Milgram theorem [START_REF] Boldo | A Coq formal proof of the Lax-Milgram theorem[END_REF]. Using the same backbone library, we formalise an alternative version of the theorem: our version is significantly more concise, and closer to the computational content of the result. We describe below this formalisation.

The Coquelicot library formalises topological concepts using filters [START_REF] Bourbaki | General Topology[END_REF][START_REF] Hölzl | Type classes and filters for mathematical analysis in isabelle/hol[END_REF], which we briefly recall here. A filter on a type T is a collection of collections of inhabitants of T which is non-empty, upward closed and stable under finite intersections:

Record Filter (T : Type) (F : (T → Prop) → Prop) := { filter_true : F (fun _ => True) ; filter_and :

∀ P Q : T → Prop, F P → F Q → F (fun x => P x /\ Q x) ; filter_imp : ∀ P Q : T → Prop, (∀ x, P x → Q x) → F P → F Q }.
While filters are used to formalise neighbourhoods, balls allow for expressing the relative closeness of points in the space. Balls are formalised using a ternary relation between two points in the carrier type, and a real number, with the following axioms: A filter is called a Cauchy filter when it contains balls of arbitrary (small) radius:

Definition cauchy (T : UniformSpace) (F : (T → Prop) → Prop) := ∀ eps > 0, ∃ x, F (ball x eps).

Finally, a uniform space is a type equipped with a ball relation and a complete space moreover has a limit operation on filters, which ensures the convergence of Cauchy sequences via the following axioms (where ProperFilter F is equivalent to Filter F /\ ∀ P, F P → ∃ x, P x): lim : ((T → Prop) → Prop) → T ; ax1 : ∀ F, ProperFilter F → cauchy F → ∀ eps > 0, F (ball (lim F) eps) ; ax2 :

∀ F1 F2, F1 ⊆ F2 → F2 ⊆ F1 → close (lim F1) (lim F2)
The above formal definition of balls does not enforce closedness nor openness. We thus introduced the relation associated with the closure of balls, so as to model closed neighbourhoods:

Definition cball x r y := ∀ e > 0, ball x (r+e) y.

Equipped with this definition, hypothesis (ii) of Theorem 1 is formalised as follows:

Definition lipschitz_on (F : U → U) (mu : R) (P : U → Prop) := ∀ x y : U, ∀ r ≥ 0, P x → P y → cball x r y → cball (F x) (mu*r) (F y).

We now sketch our formalised proof, using mathematical notations. We consider a complete space X and we write y ∈ B(x, r) for the formal (ball x r y), and y ∈ B(x, r) for (cball x r y). The key notion is that of strongly stable ball (see Figure 1): u,r) and if there is a non-negative real number s, called the offset, s.t.:

Definition 2 (Strongly stable ball). A ball B(u, r) is µ-strongly stable for F if F is µ-Lipschitz on B(
F • u ∈ B(u, r)
and s + µr ≤ r.

Remark 3 (Stability). For any x in B(u, r), a strongly stable ball for F , F • x ∈ B(u, r).

Remark 4 (Contracting case). When 0 ≤ µ < 1, for any µ-strongly stable ball B(v, ρ), with offset σ, B(F • v, µρ) is also a strongly stable ball, with offset µσ. Moreover, Assume that F has a µ-strongly stable ball B(u, r) of offset s, with µ < 1. In particular, F is contracting on B(u, r). Consider the sequence of balls defined as follows:

B(F • v, µρ) is included in B(v, ρ). u F • u F 2 • u r s μ r μ s B ( u , r )
B n = B(u n , r n ) with u n = F n • u and r n = rµ n
where F n • u denotes the iterated images of u under F . By Remark 4, (B n ) n∈N is a nested sequence of µ-strongly stable ball for F , with offset sµ n . Let F be the family of collections of points in U defined as:

F = {P ⊆ U | ∃ n, B n ⊆ P }.
It is a proper filter: F contains U , it is obviously upward closed, and for P, Q ∈ F, P ∩ Q is also in F because (B n ) n∈N is decreasing for inclusion. Thus F has a limit w, such that for any ε > 0, balls B n are eventually included in B(w, ε). We provide a formal proof of Theorem 5, a reformulation of Theorem 1 using the vocabulary of the Coquelicot library: Theorem 5. The limit w of the filter F is in B 0 , and w is a fixed point of F . Moreover, w is close to every other fixed point of F in B 0 . Proof. In this statement "w is a fixed point of F " means "w is close to F • w". First, w ∈ B n for all n. Indeed, for any ε > 0, there is an m ≥ n s.t. B m ⊆ B(w, ε), and since B m ⊆ B n , u m ∈ B n ∩ B(w, ε). In particular, w ∈ B 0 . It is also clear by stability that F • w ∈ B n for all n. Moreover, w is close to any point v s.t. v ∈ B n for all n (for any ε > 0, choose n s.t. 2µr n < ε). Taking v := F • w proves that w is a fixed point of F .

Finally, if w ∈ B 0 is another fixed point of F , then it follows from an easy induction that w ∈ B n for all n. Hence, the foregoing shows that w is close to w .

Strongly stable balls model the requirements set on the untrusted data to be formally verified. They can also be seen as balls centered at the initial point, and large enough to include all its successive iterates, i.e. as instances of the locus at stake in classical presentations of the proof. The version proved by Boldo et al. has a slightly more technical wording, which seems to be made necessary by its further usage in the verification of the Lax-Milgram theorem. Our proof script is significantly shorter, partly because we automate proofs of positivity conditions (for radii of balls) using canonical structures for manifestly positive expressions. But the key ingredient for concision is to make most of the filter device in the proof, and to refrain from resorting to low-level properties of geometric sequences. To the best of our knowledge, the other libraries of formalised analysis featuring a proof of this result, notably Isabelle/HOL and HOL-Light, are based on variant of proof strategy closer to the approach of Boldo et al. than to ours.

Newton-like validation operators

The purpose of this section is twofold. We first present the general principle of fixed-point based a posteriori validation methods, and more particularly, the use of Newton-like validation operators. Then we apply it to the division and square root of models. Throughout this section, let (X, • ) denote a Banach space, and h * the exact solution of an equation in X. In this article, X stands for the space C(I) of real-valued continuous functions defined over a compact segment I = [a, b], with the uniform norm h := sup x∈I |h(x)|. The division and square root of functions are simple examples of solutions of equations in C(I), but there are also differential equations, integral equations, delay equation, etc. The general scheme for Banach fixed-point based a posteriori validation methods follows two steps:

1. Approximation step. A numerical approximation h • ∈ X of h * is obtained by an oracle, which may resort to any approximation method. In particular, this step requires no mathematical assumption and can be executed purely numerically outside the proof assistant, good approximation properties being only desirable for efficiency. In our setting with X = C(I), interpolation at Chebyshev nodes (Section 5.1) is an efficient and accurate oracle for a wide range of function space problems. 2. Validation step. The initial problem is rephrased in such a way that h * is a fixed point of a (locally) contracting operator F : X → X. An a posteriori error bound on h • -h * is deduced from the Banach fixed-point theorem (Theorem 1).

We thus need to find a contracting operator F of which h * is a fixed point. To this end, we use Newton-like validation methods, which transform an equation M • h = 0 into an equivalent fixed-point equation F • h = h with F contracting. More specifically, suppose that M : X → Y is differentiable; we use a Newton-like operator F : X → X defined as:

F • h = h -A • M • h, h ∈ X,
with A : Y → X an injective bounded linear operator, intended to be close to (DM h • ) -1 .

The operator A may be given by an oracle and does not need to be this exact inverse, which anyway might be non representable on computers exactly. The mean value theorem yields a Lipschitz ratio µ for F over any convex subset S of X:

∀h 1 , h 2 ∈ S, F •h 1 -F •h 2 ≤ µ h 1 -h 2 , with µ = sup h∈S DF h = sup h∈S 1 X -A•DM h ,
which is expected to be small over some neighbourhood of h • . Concretely, in order to apply Theorem 1, one needs to compute the following quantities:

a bound b ≥ A • M • h • = h • -F • h • ; a bound µ 0 ≥ 1 X -A • DM h • = DF h • ; a bound µ (r) ≥ A • (DM h -DM h • ) = DF h -DF h • valid for any h ∈ B(h • , r),
and parameterised by a radius r ∈ R + .

If we are able to find a radius r ∈ R + satisfying:

µ(r) := µ 0 + µ (r) < 1, and b + rµ(r) ≤ r, (4)
then Theorem 1 guarantees the existence and uniqueness of a root h * of M in B(h • , r).

Remark 6. Finding an r as small as possible while satisfying (4) may be an nontrivial task for automated validation procedures. For many problems, µ (r) is polynomial, hence conditions (4) are polynomial inequalities over r: this is called the radii polynomial approach [START_REF] Hungria | Rigorous numerics for analytic solutions of differential equations: the radii polynomial approach[END_REF] in rigorous numerics. In our case, division (resp. square root) induces an affine (resp. quadratic) equation, which admits closed form solutions.

Approximation step: interpolation

Since they are certified a posteriori, (non-rigorous) approximations for division and square root of given models can be obtained using arbitrary numerical techniques. We use interpolation at Chebyshev nodes of the second kind for its efficiency and excellent approximation properties [START_REF] Nicholas | Approximation Theory and Approximation Practice[END_REF].

Ideally, we would implement this operation outside of the proof assistant, in order not to pay the price of an interpreted language. This would however require a of work in order to design a proper interface between Coq values and external values (e.g., converting Coq representation of floating points numbers into machine level floating points, and back). Instead, and for now, we implement the oracles inside Coq, as unspecified functions. To this end, we add a field to the structure BasisOps_on, to compute interpolants of a given degree:

interpolate: nat → (C → C) → seq C;
We implement this operation for Chebyshev basis using the discrete orthogonality relations on Chebyshev polynomials.

To reduce the price of staying inside Coq for those computations, we exploit the polymorphism built in our framework to perform those computations on floating-point numbers rather than intervals. To this end, we add the following fields to the structure NBH: The field FF, of type Ops1, will make it possible to call the functions interpolate and beval from the basis with C as FF, i.e., to let them operate on floating point numbers. By doing so we do not have to reimplement Clenshaw's evaluation scheme on floating point numbers.

Validation step for division

For f, g ∈ C(I) with g nonvanishing over I, the quotient f /g is the unique root of M : h → gh -f . Let h • be a candidate approximation given by the approximation step. Constructing the Newton-like operator F requires an approximation A of (DM h • ) -1 : k → k/g. For that purpose, suppose w ≈ 1/g ∈ C(I) is also given by an oracle, and define:

F • h = h -w(gh -f ). (5)
The next proposition computes an upper bound for h • -f /g ; it is implemented in div.newton.

Proposition 7. Let f, g, h • , w ∈ C(I), and µ, b ∈ R + such that:

(7 i) w(gh • -f ) ≤ b, (7 ii) 1 -wg ≤ µ, (7 iii) µ < 1.
Then g does not vanish over I and h

• -f /g ≤ b/(1 -µ).
Proof. Conditions (7 ii) and (7 iii) imply that F (Equation ( 5)) is contracting over C(I) with ratio µ. The radius r := b 1-µ makes the ball B(h • , r) strongly stable with offset b (7 i), since b + µr = r. Therefore, h * is the (global) unique root of M , and h • -h * ≤ r.

Finally, w and g do not vanish because 1 -wg ≤ µ < 1. Hence, h * = f /g over I.

The concrete division of models is implemented as follows: let p, q := mcf F, mcf G in mdiv_aux F G (mfc (interpolate n (fun x => beval p x / beval q x))) (mfc (interpolate n (fun x => 1 / beval q x))).

Note that we use the trivial model mbot={|pol:=[];rem:=bot|} as a default value, when the concrete computations fail to validate the guess of the oracle (either because this guess is just wrong, or because of over-approximations in the computations). The correctness lemmas use the properties of operations on models to prove the assumptions of div.newton.

Lemma rmdiv_aux F f G g H h W w: mcontains F f → mcontains G g → mcontains H h → mcontains W w → mcontains (mdiv_aux F G H W) (f/g). Lemma rmdiv n F f G g: mcontains F f → mcontains G g → mcontains (mdiv' n F G) (f/g).

Validation step for square root

Let f ∈ C(I) be strictly positive over I. The square root √ f is one of the two roots of the quadratic equation M • h := h 2 -f = 0 (the other being -√ f ). Let h • be a candidate approximation. Since DM h : k → 2hk, one also needs an approximation w ≈ 1/(2h • ) ≈ 1/(2 √ f ) ∈ C(I) in order to define A : k → wk, approximating (DM h • ) -1 . Then:

F : h → h -w(h 2 -f ). (6)
The next proposition (implemented by sqrt.newton), computes an upper bound for h • -√ f . Proposition 8. Let f, h • , w ∈ C(I), µ 0 , µ 1 , b ∈ R + and t 0 ∈ I such that:

(8 i) w h •2 -f ≤ b, (8 ii) 1 -2wh • ≤ µ 0 , (8 iii) w ≤ µ 1 , (8 iv) µ 0 < 1, ( 8
v) (1 -µ 0 ) 2 -8bµ 1 ≥ 0, ( 8 vi
) w(t 0 ) > 0.

Then f > 0 over I and h • -√ f ≤ r * , where:

r * := 1 -µ 0 -(1 -µ 0 ) 2 -8bµ 1 4µ 1 .
Proof. First, since 1 -2wh • ≤ µ 0 < 1 (by (8 ii) and (8 iv)) and w(t 0 ) > 0 (8 vi), w and h • are strictly positive over I, by continuity. Using (8 iii), µ 1 > 0.

If b = 0, then r * = 0 and h • = √ f over I, because w(h •2 -f ) = 0 (8 i) and w, h • > 0. Hence the conclusion holds.

From now on, we assume b > 0. F is Lipschitz of ratio µ(r) := µ 0 + 2µ 1 r over B(h • , r)

for any r ∈ R + , because:

F •h 1 -F •h 2 = (h 1 -h 2 )-w(h 2 1 -h 2 2 ) = [(1 -2wh • ) + w(h • -h 1 ) + w(h • -h 2 )] (h 1 -h 2 )
. Therefore, satisfying b + µ(r)r ≤ r is equivalent to the quadratic inequality:

2µ 1 r 2 + (µ 0 -1)r + b ≤ 0. (7)
Condition (8 v) implies that (7) admits solutions, and r * is the smallest one. Moreover, since b, µ 1 > 0, we get r * > 0, so that b + µ(r * )r * = r * also implies µ(r * ) < 1. Now, all the assumptions of Theorem 1 are fulfilled. Hence, F has a unique fixed point h * in B(h • , r * ). To obtain h * = √ f over I, it remains to show that h * > 0. This follows from w > 0 and:

1 -2wh * ≤ 1 -2wh • + 2w(h * -h • ) ≤ µ 0 + 2µ 1 r * = µ(r * ) < 1.
Remark 9. Contrary to the case of division where continuity was not needed at all, it is here used for w. Therefore, sqrt.newton requires w to be continuous over I.

The Coq code for the corresponding operations on models msqrt_aux and msqrt, together with the statements of their correctness lemmas, are given in Appendix B.

Examples

Playing with approximations of the absolute value function

Consider the function f ε : x → √ ε + x 2 over [-1, 1], with ε > 0. When ε → 0, f ε converges uniformly to the absolute value function x → |x| (which is not analytic at 0), with:

|f (x) -|x|| = ε + x 2 - √ x 2 = ε √ ε + x 2 + √ x 2 ≤ √ ε. ( 8
)
Rigorous uniform approximations. Approximating f ε with polynomials becomes harder for small ε, due to the complex singularities ±i √ ε getting closer to the interval [-1, 1]. Nevertheless, Chebyshev interpolation still works and our implementation computes rigorous approximations as accurate as desired (see Figure 2b), of exponential convergence with ratio determined by ε. Note that for too small degree, the computed approximation of the square root is too far from the solution, and the a posteriori validation returns an infinite remainder.

In order to provide a comparison with CoqApprox's Taylor models, we used the tactic interval with (i_depth 1, i_bisect_taylor x N, i_prec p) to build a Taylor model of degree N with precision p. Timings given in Table 2c reveal a significant advantage of our implementation (there we use ε = 2 to avoid convergence issues of Taylor models). Concerning accuracy, our experiments tend to show that when ε ≤ 1, CoqApprox fails to compute converging Taylor models. Indeed, even with large L, a goal like:

Goal Fail : ∀ x : R, -1 ≤ x ≤ 1 → sqrt (1/100+x*x) ≤ L
is not solved when the degree N becomes too large, probably indicating that the Taylor models diverge due to complex singularities inside the unit disk. (Note that the interval tactic can solve this goal, but only by resorting to subdivision techniques.)

Error bounding. We want to bound |f ε (x) -|x|| for x ∈ [-1, 1] without making use of any symbolic manipulation like [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF]. At first glance, one can choose to use the rigorous approximations over [-1, 1] obtained previously, and evaluate f ε (x) -x (resp. f ε (x) + x) over [0, 1] (resp. [-1, 0]) using Clenshaw's algorithm. However, even if the approximations are quite good, this evaluation strategy gives huge overestimations because [0, 1] and [-1, 0] are not small intervals. Instead, we compute separately two approximations for f ε : one over [0, 1] and one over [-1, 0], and we evaluate f ε (x) -x (resp. f ε (x) + x) over [1, 0] (resp. [-1, 0]) using the Chebyshev range function. This approach yields bounds that are rather close to the optimal √ ε (see Figure 2d). However, this does not allow for arbitrary accuracy: a subdivision procedure would be necessary here.

Evaluating an Abelian integral

Abelian integrals naturally appear when computing the number of limit cycles bifurcating from a Hamiltonian polynomial vector field in the plane. Indeed, the number of sign alternations of those contour integrals (parameterised by the energy level of the potential function) gives a lower bound on the number of limit cycles of the perturbed system, which is a hard question related to Hilbert's 16th problem.

In [START_REF] Johnson | A quartic system with twenty-six limit cycles[END_REF], the author claims to prove the existence of 26 limit cycles for a well-constructed quartic system, whereas the previous record for degree 4 was 22 [START_REF] Christopher | Estimating limit cycle bifurcations from centers[END_REF]. However, the implementation with which the Abelian integrals were "rigorously" computed was erroneous, which led to apparently more sign alternations than in reality. By tuning the coefficients and computing the integrals with another rigorous numerics library, the authors of the ongoing work [START_REF] Bréhard | A New Lower Bound on the Hilbert Number for Quartic Systems[END_REF] obtain 24 limit cycles, which, if not 26, is still greater than the current record 22.

To conclude this article, we rigorously evaluate some of these integrals inside Coq to show how our implementation behaves on non-crafted examples. Below are the formulas defining a family of integral I ij (r) which need to be computed precisely for several values of r. Table 1 summarises the results of our computational experiments. In each line, we chose parameters that were enough to obtain the desired precision. These encouraging results give us hope that it will be possible to fully verify the critical computations involved in recent work of the first author [START_REF] Bréhard | A New Lower Bound on the Hilbert Number for Quartic Systems[END_REF].

I ij (r) = x+ x-
x i (y + (x) j-1 -y -(x) j-1 ) dx + y+ y-(x -(y) i-1 + x + (y) i-1 )y j y 2 -y 0 δ x (y) dy.

x 0 = 9 10 , x ± = x 0 ± r/ √ 2 , δ y (x) = r 2 -(x 2 -x 0 ) 2 , x ± (y) = x 0 ± δ x (y) , y 0 = 11 10 , y ± = y 0 ± r/ √ 2 , δ x (y) = r 2 -(y 2 -y 0 ) 2 , y ± (x) = y 0 ± δ y (x) . 1 Reached precision for Iij(r) for different values of r, computed with degree-N Chebyshev models and floating-point precision p (in each cell we display the width of the computed enclosure).

Conclusion and future work

The Coq development is available online [START_REF] Bréhard | Web appendix to the present paper[END_REF]. It consists of around 1300 lines of specifications and 1500 lines of proofs. We leave several directions for future work: integrate it with CoqInterval to benefit from its automatic subdivision techniques; interface the library with external tools for the approximation steps; implement other bases; address higher-dimensional problems. Applying this approach to verify solutions of linear ODEs in a systematic way [START_REF] Benoit | Rigorous uniform approximation of D-finite functions using Chebyshev expansions[END_REF][START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF] is also a longer-term perspective.

  Record BasisOps_on (C: Type) := { lo, hi: C; (* bounds for the domain *) beval: seq C → C → C; (* (efficient) evaluation *) bmul: seq C → seq C → seq C; (* multiplication *) bone, bid: seq C; (* constant to 1, identity *) bprim: seq C → seq C; }. (* primitive *) brange: seq C → C*C; }. (* range *) Definition BasisOps := ∀ C: Ops1, BasisOps_on C.

  Definition madd (C: Ops1) (M N: Model C): Model C := {| pol := pol M + pol N; rem := rem M + rem N |}. Definition meval (C: Ops1) (M: Model C) (X: C): C := beval (pol M) X + rem M. Definition mintegrate (C: Ops1) (M: Model C) (a b: C): C := let N := bprim (pol M) in beval N b -beval N a + (b-a)*rem M. For those relatively simple operations, it suffices to have the basic operations (Ops1) on C. For other operations like the range of a model, we actually need the additional operations on intervals provided by the structure NBH: Definition mrange (N: NBH) (M: Model II) := let (a,b) := brange (pol M) in bnd a b + rem M.

Canonical

  Structure MFunOps (N: NBH) (B: BasisOps): FunOps II. (* with carrier [Model II] *)

Fixpoint

  Clenshaw b c (p: seq C) x := match p with | [] => c -x*b | a::q => Clenshaw c (a + 2*x*c -b) q x end. Definition beval (p: seq C) x := Clenshaw 0 0 (rev p) x. This code might look mysterious; it is justified by the following invariant on real numbers: Lemma ClenshawR b c p x: Clenshaw b c p x = eval T (catrev p [c -2*x*b; b]) x.

Fixpoint

  mul_pls (p q: seq C): seq C := match p,q with | [],_ | _,[] => [] | a::p', b::q' => sadd (a*b::(sadd (sscal a q') (sscal b p'))) (0::0::mul_pls p' q') end. Fixpoint mul_mns (p q: seq C): seq C := match p,q with | [],_ | _,[] => [] | a::p', b::q' => sadd (a*b::

  Definition range_: seq C → C := foldr (fun A X => abs A + X) 0. Definition range (P: seq C): C*C := match p with | [] => (0,0) | A::Q => let R := range_ Q in (A-R,A+R) end.

  ball : M → R → M → Prop ; ax1 : ∀ x (e > 0), ball x e x ; ax2 : ∀ x y e, ball x e y → ball y e x ; ax3 : ∀ x y z e1 e2, ball x e1 y → ball y e2 z → ball x (e1 + e2) z Two points are called close when they cannot be separated by balls: Definition close (x y : M) : Prop := ∀ eps > 0, ball x eps y.

Figure 1

 1 Figure 1 Balls B0 and B1.

FF

  type for floating points and their operations *) I2F: II → FF;(* conversion from intervals to floating points to (midpoint) *) F2I: FF → II;(* conversion from floating points to intervals (singleton) *)Equipped with these operations, we can define conversion operations between models (on intervals) and polynomials with floating point coefficients:Definition mcf (M: Model): seq FF := map I2F (pol M). Definition mfc (p: seq FF): Model := {| pol := map F2I p; rem := 0 |}.

  Definition mdiv_aux (F G H W: Model): Model := let K1 := 1-W*G in let K2 := W*(G*H -F) in match mag (mrange K1), mag (mrange K2) with | Some mu, Some b when is_lt mu 1 => {| pol := pol H; rem := rem H + sym (b/(1-mu)) |} | _ => mbot end. Definition mdiv n (F G: Model): Model :=

  Functions fε and x → |x| over [-1, 1]. Magnitude of the remainders of degree-N Chebyshev models approximating fε over [-1, 1].

( c )

 c Timings of degree-N models for f2. Overapproximation ratio of the remainder of the degree-N Chebyshev model for fε -1 over [0, 1], compared to the optimal bound √ ε.

Figure 2

 2 Figure 2 Approximating functions fε and x → |x| with Chebyshev models.

  r

Funding This work has been funded by the European Research Council (ERC) under the European Union's Horizon 2020 programme (CoVeCe, grant agreement No 678157), and was supported by the LABEX MILYON (ANR-10-LABX-0070) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR). This work was supported in part by the project FastRelax ANR-14-CE25-0018-01 1 Record Rel0 (C D: Ops0) := { rel:> C → D → Prop; radd: ∀ x y, rel x y → ∀ x' y', rel x' y' → rel (x+x') (y+y'); rsub: ∀ x y, rel x y → ∀ x' y', rel x' y' → rel (x-x') (y-y'); rmul: ∀ x y, rel x y → ∀ x' y', rel x' y' → rel (x*x') (y*y'); rzer: rel 0 0; rone: rel 1 1 }.

gitlab.inria.fr/amahboub/approx-models

A Coq code for primitive in Chebyshev basis

Fixpoint prim_ (C: Ops1) (n: nat) (p: seq C): seq C := match P,n with | [],_ => [] | a::q,0 => sadd [0; a] (prim_ 1 q) | a::q,1 => 0 :: (sadd [0; a/4] (prim_ 2 q)) | a::q,_ => sadd [-a/(2(n-1)); 0; a/(2(n+1))] (0 :: prim_ (n+1) Q) end. Definition prim (C: Ops1) (P: seq C) := prim_ 0 P.

The key lemma is the following one:

Lemma eval_prim_ n (p: seq R) x: Derive (eval_ T (n-1) (prim_ n P)) x = eval_ T n P x.

B

Coq code for the square root of a model Lemma rmsqrt_aux (F H W: Model) (X: II) (f h w : R → R) (x: R): mcontains F f → mcontains H h → mcontains W w → contains X x → lo≤x≤hi → (∀ x, lo≤x≤hi → continuity_pt w x) → mcontains (msqrt_aux F H W X) (sqrt f).

Lemma rmsqrt n F f: mcontains F f → mcontains (msqrt' n F) (sqrt f).