Érik Martin-Dorel
email: erik.martin-dorel@ens-lyon.org

Guillaume Hanrot
email: guillaume.hanrot@ens-lyon.fr

Micaela Mayero
email: micaela.mayero@lipn.univ-paris13.fr

Laurent Théry
email: laurent.thery@inria.fr

Keywords: formal proofs, certificate checkers, Hensel's lemma, modular arithmetic

In order to derive efficient and robust floating-point implementations of a given function f , it is crucial to compute its hardest-to-round points, i.e. the floatingpoint numbers x such that f (x) is closest to the midpoint of two consecutive floatingpoint numbers. Depending on the floating-point format one is aiming at, this can be highly computationally intensive. In this paper, we show how certificates based on Hensel's lemma can be added to an algorithm using lattice basis reduction so that the result of a computation can be formally checked in the Coq proof assistant.

Formally verified certificate checkers for hardest-to-round computation 1 Introduction and motivations

Hensel's lemma is a very classic mathematical tool. Given a polynomial P over the integers, it starts from a solution to P (x) ≡ 0 (mod p) and constructs a solution to P (x) ≡ 0 (mod p k ) for any k ∈ N * , under mild hypotheses (for instance, that p does not divide the discriminant of P ). This lemma was given its final shape by Hensel, in his study of p-adic numbers, where it plays a central role [START_REF] Hensel | Neue Grundlagen der Arithmetik[END_REF]; see also [START_REF] Steuding | Diophantine Analysis[END_REF]Chapter 13]. The algorithm derived from the proof of Hensel's lemma is usually called Hensel lifting.

The strengths of Hensel's lemma are its versatility (it can be stated in a very general context, being also valid, e.g., for power series rings), the fact that it combines an assertion of existence (there is a lifting) and uniqueness (this lifting is unique) under mild hypotheses, and the fact that it is effective, and provides a simple and efficient algorithm. Hensel's lemma has applications in several areas of mathematics and computer science where one has to solve equations over the integers (or over the ring of polynomials over a field, as in the last step of Guruswami-Sudan list decoding algorithm for Reed-Solomon codes [START_REF] Guruswami | Improved decoding of Reed-Solomon and algebraic-geometry codes[END_REF][START_REF] Augot | A Hensel Lifting to Replace Factorization in List-Decoding of Algebraic-Geometric and Reed-Solomon Codes. Information Theory[END_REF]). In this work, we are interested in some of them that we present just below.

Hensel's lemma in computer algebra

Since the beginning of the development of computer algebra, the algorithmic version of Hensel's lemma (Hensel lifting) has been one of the fundamental tools for computations over the integers. The strategy is to isolate a "nice" prime p where the problem under study over the integers has a "good reduction" to the same problem modulo p k for all k ∈ N * . One then solves the problem modulo p, lifts the solution modulo p m for some m such that p m is slightly larger than an a priori bound on the size of the solution.

Standard applications include linear system solving over the integers (Hensel's technique is very powerful in practice as it avoids manipulating huge rational numbers while keeping the size of integers involved under control), but also 0-dimensional polynomial systems solving over the integers (though the poor quality of the bounds for real solutions as the number of variables grows might limit the applicability in that case) or polynomial factoring over the integers. In all those cases, Hensel's lemma can be used either to compute the solutions or to provide a certificate for a list of solutions (using the uniqueness statement of the lemma).

In the present paper, we are mainly interested in the second aspect, in a setting where Hensel lifting is used at the end of a complex algorithmic chain where formal verification is a major issue; we now turn to the description of the problem for which this algorithmic chain was designed: the Table Maker's Dilemma, a problem which needs to be solved for implementing elementary functions with correct rounding.

The Table Maker's Dilemma

The newly revised IEEE 754-2008 standard for floating-point (FP) arithmetic recommends that some mathematical functions (exp, log, x → 2 x , . . . ) should be correctly rounded: the system must always return the FP number that is the nearest to the exact mathematical result of the operation. Having correctly rounded functions greatly improves the portability of numerical software.

If we consider a correctly-rounded elementary function f in radix 2 and precision n, most implementations use the so-called Ziv's strategy:

1. compute a small enclosing interval of width 2 -n ′ for f (x); 2. if both endpoints of this interval round at precision n to the same FP number v, this number v is the correctly rounded value of f (x); 3. otherwise, increase n ′ , and go to step 1.

The correctness of this strategy follows from the fact that rounding is monotonic-if both endpoints of the interval round to the same value, any point of the interval also does. Note that in order to avoid an infinite loop, one should be able to deal separately with "exact cases", i.e. the FP numbers x such that (for rounding-tonearest) f (x) is the exact midpoint of two consecutive FP numbers, or (for directed rounding modes) f (x) is a FP number, like x = 0 for the exp function.

In order to make this strategy optimal, we want to minimize the number of iterations through step 3, and set the n ′ to a large enough precision to guarantee that the test at step 2 is true, and yet not too large, since the cost of step 1 increases significantly with n ′ . Finding the optimal value (say n ′ 0 ) of this "large enough" precision is the so-called Table Maker's Dilemma (TMD) [30, chap. 12]. Finding an upper bound for n ′ 0 will be called the Approximate TMD. Solving the TMD, for a given function f and a rounding mode (say, rounding-to-nearest) reduces to finding the hardest-to-round (HR) points, that is the FP numbers x such that f (x) is closest to the midpoint of two consecutive FP numbers.

In this paper, we are interested in one of the algorithms that solve this problem: the Stehlé-Lefèvre-Zimmermann (SLZ) algorithm [START_REF] Stehlé | Searching Worst Cases of a One-Variable Function Using Lattice Reduction[END_REF][START_REF] Stehlé | Algorithmique de la réduction des réseaux et application à la recherche de pires cas pour l'arrondi des fonctions mathématiques[END_REF]. This algorithm requires complex and very long calculations (years of cumulated CPU time), and its only output as of today is a list of hard-to-round points (each of which gives a lower bound on n ′ 0 ) and a claim that the worst one is in that list. This is even worse when the SLZ algorithm is applied to the Approximate TMD, as the output is limited to a yes/no answer. This situation thus appears to be unsatisfactory: the initial motivation for the TMD problem is to provide strong guarantees for numerical computations, but the trusted computing base contains highly optimized implementations of complicated algorithms. One of the main contributions of this paper is to propose "logs" for the execution of the SLZ algorithm, allowing one to use aggressively optimized implementations of computer algebra primitives while making it easy to verify the results.

The critical step of the SLZ algorithm consists in generating two polynomials v 1 , v 2 (∈ Z[X, Y ]) as Z-linear combinations of a specific polynomial family. These polynomials are hard to find-they are produced by an aggressively optimized version of the complicated LLL [START_REF] Lenstra | Factoring polynomials with rational coefficients[END_REF] algorithm, which, for standard parameter values is by far the most time-consuming part of the SLZ algorithm. It is then possible to treat the LLL calls as oracles and log their results in a certificate. To be more precise, the SLZ algorithm consists in the following main steps. First, it performs a domainsplitting step and calls a polynomial approximation algorithm. In doing so, the initial (Approximate)-TMD problem is rewritten into myriads of problems that involve integers only. All problems are then instances of the Integer Small Value Problem [START_REF] Stehlé | Searching Worst Cases of a One-Variable Function Using Lattice Reduction[END_REF]Section 3.2]:

Problem 1 (ISValP) For P ∈ Z[X] and (A, B, M ) ∈ N 3 , find all solutions of the system          (x, y) ∈ Z 2 , |x| A, |y| B, P (x) ≡ y (mod M ).
Next, each instance of ISValP is solved by successively using the modular, bivariate version of Coppersmith's technique, and the bivariate version of Hensel lifting. This is where the concept of ISValP certificate comes into play, and the core motivation for our work.

Outline of the paper

In Section 2, we first present the bivariate version of Hensel lifting and then focus on the formalization of the uniqueness statement of bivariate Hensel's lemma in the Coq formal proof assistant. In Section 3, we give the formalization of two interrelated certificate checkers, to address instances of the bivariate small integral roots problem as well as the integer small value problem (ISValP). Next, in Section 4, we present some extra formal material that is used to increase the efficiency and the modularity of these certificate checkers. We then evaluate the performances of our ISValP certificate checker on several examples. In Section 5, we discuss some salient technicalities related to our formalization. Finally we draw some conclusions in Section 6.

2 Formal study of bivariate Hensel's lemma

Preliminary definitions and notations

For any x ∈ R, the largest integer x is denoted by ⌊x⌋ and the smallest integer x by ⌈x⌉, so that we have ⌈x⌉ = -⌊-x⌋. For any (a, b) ∈ Z 2 (assuming a b), the set of integers k satisfying a k b is denoted by a, b . Then for any x ∈ R and q ∈ ]0, +∞[, the "standard modulo" of x by q (also called the remainder of x modulo q) taken in [0, q[ is denoted by

x mod q := x -q x q ( 0). (1)
As a result, we straightforwardly get x mod q ∈ 0, q -1 for any x ∈ Z and q ∈ N * . Finally for any (a, b) ∈ R 2 and q ∈ ]0, +∞[, the usual abbreviation for modular equality is used:

a ≡ b (mod q) ⇐⇒ a mod q = b mod q.
Definition 1 (Small integral root) For any pair (v 1 , v 2 ) ∈ Z[X, Y ] 2 of bivariate polynomials over Z and any pair (A, B) ∈ N 2 , a small integral root of (v 1 , v 2 ) (with respect to bounds A and B) is a pair (x, y) ∈ Z 2 that is solution of the system

         |x| A, |y| B, v 1 (x, y) = 0, v 2 (x, y) = 0.
Then the problem of finding all the solutions of this system is referred to as the bivariate small-integral-roots problem (SIntRootP).

Definition 2 (Modular root) For any pair

(v 1 , v 2 ) ∈ Z[X, Y ] 2 and any q ∈ N * , a modulo-q root of (v 1 , v 2 ) is a pair (x, y) ∈ Z 2 satisfying v 1 (x, y) ≡ 0 (mod q), v 2 (x, y) ≡ 0 (mod q). Note that if (x, y) is a modulo-q root of (v 1 , v 2 ) ∈ Z[X, Y ] 2 , then any (z, t) ∈ (x + qZ) × (y + qZ) is also a modulo-q root of (v 1 , v 2 ).

Definition 3 (Modular inverse of an integer)

For any x ∈ Z coprime with q ∈ N * , the modulo-q inverse of x (taken in 0, q -1 ) is the integer x -1 q ∈ 0, q -1 satisfying ∃y ∈ Z, xx -1 q + yq = 1. This quantity can be computed by using the extended Euclidean algorithm.

Definition 4 (Modular inverse of a matrix)

For any M ∈ M n (Z) such that det M is coprime with q ∈ N * , the modulo-q inverse of M can be defined as

M -1 q := (det M ) -1 q • (adj M ) mod q,
where the modulo operation is applied coordinatewise and where adj M denotes the adjugate matrix of M , that is, the transpose of the cofactor matrix of M , cf. [START_REF] Stewart | On the adjugate matrix[END_REF]. In particular, for order-2 matrices over Z, we get:

a b c d -1 q = (ad -bc) -1 q • d -b -c a mod q, ( 2
)
assuming ad -bc is coprime with q ∈ N * . Note that if q = p 2 k for a given prime p, this coprimality condition is equivalent to ad -bc ≡ 0 (mod p).

A slightly unusual notion in modular arithmetic is that of centered modulo:

Definition 5 (Centered modulo) For any x ∈ R and q ∈ ]0, +∞[, the centered modulo of x by q, denoted by x cmod q, satisfies |x cmod q| q 2 and ∃n ∈ Z, x cmod q = x + nq.

(

) 3
To ensure the uniqueness of this definition, we can consider the centered modulo that has values in -q 2 , q 2 , which amounts to considering the following definition:

x cmod q := x -q x q - 1 2 . ( 4
)
Note that if x and q are integers, so is x cmod q. As a matter of fact, in the following, we will only deal with centered moduli over the integers. Also, the centered modulo can be emulated with the standard one:

x cmod q = x mod q if x mod q q 2 x mod q -q otherwise. (

Furthermore, any centered modulo function trivially enjoys the following equivalence:

∀(z, q, n) ∈ Z 3 , |z cmod q| n ⇐⇒ ∃x ∈ Z, |x| min q 2 , n ∧ z ≡ x (mod q), ( 6
)
which can be useful to "unfold" the centered modulo when it occurs in such an inequality. Note also that in practice we will always have n < q/2, allowing one to further simplify the right-hand side of the equivalence (6).

Finding bivariate small-integral-roots using Hensel lifting

Before focusing on the uniqueness statement of bivariate Hensel's lemma and in order to give more intuition, we start with the algorithmic version of this lemma, called bivariate Hensel lifting. The main step of this iterative algorithm is summarized in Algorithm 1, along with the conditions and post-conditions that actually relate to the existence statement of bivariate Hensel's lemma.

Algorithm 1: Bivariate Hensel lifting: sketch of the computation to be iterated

Input : (v 1 , v 2 ) ∈ Z[X, Y ] 2 , p prime, k ∈ N, q = p 2 k , (a, b) ∈ Z 2 such that det [Jv 1 ,v 2 (a, b)] ≡ 0 (mod p) and v 1 (a, b) ≡ 0 ≡ v 2 (a, b) (mod q). Output : (a ′ , b ′ ) ∈ Z 2 such that a ′ ≡ a (mod q), b ′ ≡ b (mod q), and v 1 (a ′ , b ′ ) ≡ 0 ≡ v 2 (a ′ , b ′ ) (mod q 2 ). a ′ b ′ ← a b -Jv 1 ,v 2 (a, b) -1 q 2 v 1 (a, b) v 2 (a, b) cmod q 2
We can see that Algorithm 1 deals with the modular roots of a pair (v 1 , v 2 ) of bivariate polynomials over Z, and involves the Jacobian matrix 1 of (v 1 , v 2 ) evaluated at (a, b) ∈ Z 2 . Since p is prime, the condition det [J v1,v2 (a, b)] ≡ 0 (mod p) implies that it is legal to consider the modular inverse J v1,v2 (a, b) -1 q 2 (for q 2 = p 2 k+1 ). Note 1 Recall that this matrix is the generalization of the derivative in higher-dimensional calculus and is defined in the bivariate case by:

Jv 1 ,v 2 (a, b) = ∂ v 1 ∂X (a, b) ∂ v 1 ∂Y (a, b) ∂ v 2 ∂X (a, b) ∂ v 2 ∂Y (a, b)
. that this algorithm is very similar to the so-called Newton's method, so its univariate counterpart is called p-adic Newton iteration in the book Modern Computer Algebra [START_REF] Von | Modern Computer Algebra[END_REF]Algorithm 9.22,p. 264]. Note also that the univariate algorithm presented in this reference stores the modulo-(p 2 k ) inverses in an extra variable, which is updated at each iteration using Hensel's iteration for the inverse "on the fly": this optimization could be generalized to the bivariate case too.

To find the bivariate small-integral-roots of (v 1 , v 2 ), we just need to:

(i) Find a "nice" prime p so that Algorithm 1 is applicable; (ii) Iterate Algorithm 1 on each modulo-p root of (v 1 , v 2 ); (iii) Stop as soon as one reaches an a-priori bound on the solutions.

For item (i), we can choose any prime p such that the following condition holds:

∀(x, y) ∈ Z 2 , v 1 (x, y) ≡ 0 ≡ v 2 (x, y) (mod p) =⇒
det J v1,v2 (x, y) ≡ 0 (mod p). [START_REF] Boneh | Finding Smooth Integers in Short Intervals Using CRT Decoding[END_REF] Note that (7) can be decided by checking p values for x and p values for y, that is,

p 2 different pairs (x, y) ∈ 0, p -1 2 .
Algorithm 2 describes how to process item (ii) and item (iii), in order to solve instances of bivariate SIntRootP (according to the terminology of Definition 1).

Algorithm 2: Find the small integral roots of two bivariate polynomials

Input: (v 1 , v 2 ) ∈ Z[X, Y ] 2 , (A, B) ∈ N 2 .
Output: The small integral roots of (v 1 , v 2 ) with respect to bounds A and B. We now give a concrete example to see how Algorithm 2 works. This example comes from an actual execution of the SLZ algorithm and will thus be reused later.

1 p ← a prime that satisfies (7) 2 S ← ∅ 3 foreach (r, s) ∈ 0, p -1 2 such that v 1 (r, s) ≡ 0 ≡ v 2 (r, s) (mod p) do 4 k ← 0 5 K ← p 6 a ← r cmod K 7 b ← s cmod K // Here, K = p 2 k , a ≡ r (mod p), b ≡ s (mod p), |a| K 2 , |b| K 2 8 while K 2 • max(A, B) do 9 k ← k + 1 10 K ← K 2 11 a b ← a b -Jv 1 ,v 2 (a, b) -1 K v 1 (a, b) v 2 (a, b) cmod K 12 end // Here, K = p 2 k , a ≡ r (mod p), b ≡ s (mod p), |a| K 2 , |b| K 2 13 if v 1 (a, b) = 0 = v 2 (
Example 1 Let us consider the following bivariate polynomials:

         v 1 (X, Y ) = Q(X, Y ) 2 + (-6Y -1593703)M Q(X, Y ) + (9Y 2 + 4781109Y + 634972313052)M 2 , v 2 (X, Y ) = -25M Q(X, Y ) + (76Y + 19921289)M 2 ,
where M = 2 92 and

Q(X, Y ) = -X + (36506148256923413Y 3 + 28709603185988793532416Y 2 + 15052100435175692583523319808Y + 3945817816478696756615137147748352).
Suppose we want to determine all the integral roots (x, y) of (v 1 , v 2 ) that satisfy

|x| A := 2 87 and |y| B := 2 4 .
We first compute the Jacobian matrix of (v 1 , v 2 ),

J v1,v2 = ∂ v1 ∂X ∂ v1 ∂Y ∂ v2 ∂X ∂ v2 ∂Y ,
whose determinant cancels for some2 modulo-2 roots of (v 1 , v 2 ). In contrast, we can check that modulo 3, no modular root of (v 1 , v 2 ) is a modular root of det J v1,v2 :

             v 1 (x, y) ≡ 0 (mod 3), v 2 (x, y) ≡ 0 (mod 3)
, and det J v1,v2 (x, y) ≡ 0 (mod 3) if (x mod 3, y mod 3) ∈ {(0, 2), (2, 1)} , v 1 (x, y) ≡ 0 (mod 3) or v 2 (x, y) ≡ 0 (mod 3) otherwise.

This means that p := 3 is the first prime satisfying the key hypothesis [START_REF] Boneh | Finding Smooth Integers in Short Intervals Using CRT Decoding[END_REF]. Following Algorithm 2, we can iterate bivariate Hensel lifting on both modulo-3 roots of (v 1 , v 2 ): since 3 2 5 2 × max(A, B) < 3 2 6 , 6 iterations are performed and lead to the modular roots summarized in Table 1.

Focusing on the values computed at the 6th iteration, we can notice that both (x, y) = (233673222969575457782319, 11) and (33647808062653569260574440, -14)

satisfy v 1 (x, y) = 0 = v 2 (x, y), |x| A and |y| B,
which means we have found two small integral roots of (v 1 , v 2 ). Thanks to the material presented in this paper, it will also be possible to provide a certificate that guarantees that there is no other small integral root of (v 1 , v 2 ), that is, we have found all the solutions of the problem at stake (in particular, see Section 3 for the sequel of this example).

k (a, b) at k-th iteration for (r, s) := (0, 2) (a, b) at k-th iteration for (r, s) := (2, 1) 0 (0, -1) (-1, 1) 1 (0, 2) (-4,
(m -n)| < 2 • d, that is |m -n| < d. Since we also have m -n ≡ 0 (mod d), Lemma 1 implies that m -n = 0, that is, m = n. ⊓ ⊔
We can now present the uniqueness statement of bivariate Hensel's lemma:

Lemma 2 (Hensel) Let v 1 , v 2 be two bivariate polynomials with integer coefficients, and let p be a prime that satisfies [START_REF] Boneh | Finding Smooth Integers in Short Intervals Using CRT Decoding[END_REF]. For all (x, y) ∈ Z 2 and k ∈ N such that

v 1 (x, y) ≡ 0 ≡ v 2 (x, y) (mod p 2 k ), (8)
any (Z 2 )-valued sequence (a n , b n ) 0 n k defined by the relations

(a 0 , b 0 ) ≡ (x, y) (mod p) (9)
and

∀n ∈ 0, k -1 , a n+1 b n+1 ≡ a n b n -J v1,v2 (a n , b n ) -1 p 2 n+1 v 1 (a n , b n ) v 2 (a n , b n ) (mod p 2 n+1 ) (10) satisfies: ∀n ∈ 0, k , (a n , b n ) ≡ (x, y) (mod p 2 n ). ( 11
)
Notice that the modular equalities that appear in conditions ( 9) and ( 10) somewhat generalize the algorithms presented in the previous Section 2.3, which were involving a centered modulo. This generalization is made possible as Hensel lifting does not depend on the representatives chosen at each step.

The following corollary shows that Lemma 2 constitutes a uniqueness result for bivariate Hensel lifting:

Corollary 2 (Uniqueness of Hensel lifting) Assume

(v 1 , v 2 ) ∈ Z[X, Y ] 2 satisfies (7) for a given prime p. For any k ∈ N, if (a, b) ∈ Z 2 and (c, d) ∈ Z 2 are modulo-(p 2 k ) roots of (v 1 , v 2 ), that is if v 1 (a, b) ≡ 0 ≡ v 2 (a, b) (mod p 2 k ) v 1 (c, d) ≡ 0 ≡ v 2 (c, d) (mod p 2 k ),
then the lifting of the corresponding modulo-p roots is unique, that is

a ≡ c (mod p) b ≡ d (mod p) =⇒ a ≡ c (mod p 2 k ) b ≡ d (mod p 2 k ).
Proof (Corollary 2) First, apply Lemma 2 with (x, y) := (a, b) and fix a sequence (a n , b n ) 0 n k satisfying ( 9) and ( 10), e.g., using the standard modulo. Then take n :=

k in (11), so that (a k , b k ) ≡ (a, b) (mod p 2 k ). Second, apply Lemma 2 with (x, y) := (c, d), using the same sequence (a n , b n ) 0 n k which indeed satisfies (a 0 , b 0 ) ≡ (a, b) ≡ (c, d) (mod p). Then take n := k again, so that (a k , b k ) ≡ (c, d) (mod p 2 k ), hence the result. ⊓ ⊔
Before giving a complete mathematical proof of Lemma 2, let us mention that we give in Appendix A the mathematical proof of the correctness of Algorithm 2 that we deduce from this lemma. Although this proposition in appendix is not directly related to the main topic of the paper, it can be useful to illustrate the usefulness of Lemma 2 beyond the formal verification of certificates checkers based on Hensel's lemma. (In particular, the verification of these checkers does not require Algorithm 2 to be proved correct, as we will focus on a certifying variant of Algorithm 2 instead in Section 3.1.) Proof (Lemma 2) Let us assume that all the hypotheses of the lemma hold for

(v 1 , v 2 ) ∈ Z[X, Y ] 2 , p prime, (x, y) ∈ Z 2 ,
and (a n , b n ) 0 n k . To shorten the formulas involved in the sequel, we denote by (x n , y n ) n∈N the sequence defined for all n ∈ N by (x n , y n ) := (x mod p 2 n , y mod p 2 n ). Proving [START_REF] Cohen | Refinements for Free! In[END_REF] thus amounts to showing that ∀n ∈ N, n k ⇒ (a n , b n ) ≡ (x n , y n ) (mod p 2 n ). We prove it by induction on n:

-First, we have p = p 2 0 , so that (9) implies (a 0 , b 0 ) ≡ (x 0 , y 0 ) (mod p 2 0 ), which proves the base case. -For the inductive case, let us show that for any integer n < k that satisfies

(a n , b n ) ≡ (x n , y n ) (mod p 2 n ), we have (a n+1 , b n+1 ) ≡ (x n+1 , y n+1 ) (mod p 2 n+1 ). We can write x n = x mod p 2 n = [x mod p 2 n+1 ] mod p 2 n = x n+1 mod p 2 n , so that x n+1 ≡ x n (mod p 2 n ), hence x n+1 ≡ a n (mod p 2 n
) by using the induction hypothesis. This implies that

∃λ ∈ Z, x n+1 = a n + λp 2 n .
Likewise, we deduce that

∃µ ∈ Z, y n+1 = b n + µp 2 n .
For any bivariate polynomial P , we consider the operator

∆ i,j (P ) := 1 i!j! ∂ i+j ∂X i ∂Y j P . ( 12
) Despite the division, ∆ i,j maps Z[X, Y ] to Z[X, Y ].
This can be easily checked on monomials where we have ∆

i,j (X k Y ℓ ) = k i ℓ j X k-i Y ℓ-j .
Using this operator, we apply Taylor's theorem to each bivariate polynomial v l (l = 1, 2):

v l (x n+1 , y n+1 ) = i,j∈N λp 2 n i µp 2 n j ∆ i,j (v l )(a n , b n ). ( 13
)
The left-hand side of ( 13) is zero modulo p 2 n+1 , since

v l (x n+1 , y n+1 ) = v l (x mod p 2 n+1 , y mod p 2 n+1 ) ≡ v l (x, y) = 0 (mod p 2 n+1 ).
Concerning the right-hand side, we can notice that

∀i, j ∈ N, i + j 2 =⇒ λp 2 n i µp 2 n j = λ i µ j p 2 n i+j-2 • p 2 n+1
and

∆ i,j (v l )(a n , b n ) ∈ Z,
therefore all terms in the summation involved in ( 13) are zero modulo p 2 n+1 whenever i + j 2. As a result, (13) becomes

0 ≡ v l (a n , b n ) + λp 2 n ∂ ∂X v l (a n , b n ) + µp 2 n ∂ ∂Y v l (a n , b n ) (mod p 2 n+1 ). ( 14
)
We can combine ( 14) for both values l = 1, 2 to obtain

0 0 = v 1 (a n , b n ) v 2 (a n , b n ) + J v1,v2 (a n , b n ) λp 2 n µp 2 n (mod p 2 n+1 ), (15)
where the modulo is taken coordinatewise. Then, we can easily combine [START_REF] Brisebarre | Rigorous Polynomial Approximation Using Taylor Models in Coq[END_REF] with the induction hypothesis to deduce that

v 1 (a n , b n ) ≡ 0 ≡ v 2 (a n , b n ) (mod p), hence by (7), det J v1,v2 (a n , b n ) ≡ 0 (mod p)
and thereby det J v1,v2 (a n , b n ) ≡ 0 (mod p 2 n+1 ), so that (15) can be rewritten to

-J v1,v2 (a n , b n ) -1 p 2 n+1 v 1 (a n , b n ) v 2 (a n , b n ) ≡ λp 2 n µp 2 n (mod p 2 n+1 ). ( 16
)
Then we replace λp 2 n with x n+1 -a n (resp. µp 2 n with y n+1 -b n ) and we obtain

a n b n -J v1,v2 (a n , b n ) -1 p 2 n+1 v 1 (a n , b n ) v 2 (a n , b n ) ≡ x n+1 y n+1 (mod p 2 n+1 ).
Now we can use [START_REF] Chrząszcz | Modules in Coq Are and Will Be Correct[END_REF] to finally deduce that

a n+1 b n+1 ≡ x n+1 y n+1 (mod p 2 n+1 ). ⊓ ⊔

Bivariate Hensel's lemma: formal proof in Coq

The code of our CoqHensel library, along with some documentation, is available online at http://tamadi.gforge.inria.fr/CoqHensel/.

Let us now focus on the formalization of bivariate Hensel's lemma in Coq [START_REF] Coq | The Coq Proof Assistant: Reference Manual: version 8.4pl4[END_REF][START_REF] Bertot | Interactive Theorem Proving and Program Development. Coq'Art: The Calculus of Inductive Constructions[END_REF]. As regards the basic types, functions and lemmas required for this formalization, we can see from the previous Section 2.3 that we need:

(i)
N with the usual operations (including exponentiation) as well as the primality and divisibility predicates; (ii) Z with the usual operations (including absolute value) and a comprehensive formalization of modular arithmetic; (iii) Z[X, Y ] with polynomial evaluation, derivatives, and Taylor's theorem for bivariate polynomials; (iv) Some facilities to handle iterated operators, and in particular the summations that appear when proving (and using) Taylor's theorem; (v) ExtendedGCD-based modular inversion over M 2 (Z).

A large part of these notions are covered by the libraries distributed with the SSReflect [START_REF] Gonthier | A Small Scale Reflection extension for the Coq system[END_REF][START_REF] Gonthier | An introduction to small scale reflection in Coq[END_REF] extension of Coq. First, item (i) above is fully handled using libraries ssrnat, prime and div. For item (ii), we rely on the SSReflect libraries ssrint, intdiv and ssrnum which define integers as

Inductive int := Posz (n : nat) | Negz (n : nat).
0 is represented by (Posz 0), 1 by (Posz 1), -1 by (Negz 0), and so on. As this definition involves Peano integers (nat), it is definitely not adapted to largescale computation, but quite convenient to develop proofs and going back-and-forth between signed integers and natural numbers. For computing effectively with integers we will rely on other data structures, which will be the topic of Section 4.1.

For item (iii), we simply define bivariate polynomials as polynomials with univariate polynomials as coefficients. The corresponding Coq code is:

Notation "{ 'bipoly' R }" := {poly {poly R}} : type_scope.

Note that this definition is asymmetric with respect to the indeterminates X and Y (which are respectively defined as the outermost indeterminate and the innermost one). One of the benefits of this definition is that we directly inherit all ring operations and support theorems from the SSReflect library poly. We just need to define the n-th partial derivatives and bivariate Horner evaluation, then prove the following version of Taylor's theorem: Theorem 1 (Taylor) For any unital ring R, for any polynomial P ∈ R[X, Y ] and for any x 0 , y 0 , h, k ∈ R that commute pairwise 3 we have

P (x 0 + h, y 0 + k) = i j ∆ i,j (P )(x 0 , y 0 ) • h i • k j . ( 17
)
The operator ∆ i,j above has been previously defined in Equation ( 12), and we recall that it maps the ring R[X, Y ] to R[X, Y ], despite the presence of the division 1 i!j! . The proof of Theorem 1 essentially relies on some manipulations of summations, 3 that is, x 0 y 0 = y 0 x 0 , x 0 h = hx 0 , x 0 k = kx 0 , y 0 h = hy 0 , y 0 k = ky 0 and hk = kh.

and uses twice the univariate version of Taylor's theorem. Note that our version of univariate Taylor's theorem has been included to the standard SSReflect library poly.

For item (iv), we benefited from the formalization of generic "big operators" provided in the SSReflect library bigop [START_REF] Bertot | Canonical Big Operators[END_REF].

For item (v), we rely on the formalization of rings Z/qZ provided in the SSReflect library zmodp to address Definition 3. Then, we benefit from the SSReflect library matrix to deduce Definition 4. We formalize this definition in a generic way (for order-n matrices) before specializing it to order-2 matrices within the proof of Lemma 2. Although not essential to derive the proof of bivariate Hensel's lemma, the use of matrices leads to more concise expressions and proof steps. In particular, a key ingredient about order-2 matrices is given by the following version of Cramer rule with moduli:

∀p prime, ∀k ∈ N, ∀A ∈ M 2 (Z), ∀t ∈ Z 2 = M 2,1 (Z), det(A) ≡ 0 (mod p) =⇒ A -1 p 2 k+1 × (A × t) ≡ t (mod p 2 k+1
). ( 18) The formal proof of Lemma 2 closely follows the pencil-and-paper proof presented in Section 2.3. In particular, we rely on the result [START_REF] Gonthier | An introduction to small scale reflection in Coq[END_REF] to deduce ( 16) from [START_REF] Coq | The Coq Proof Assistant: Reference Manual: version 8.4pl4[END_REF].

Formal verification of certificates using Hensel's lemma

We now follow the certificate-based approach to address the Integer Small Value Problem (ISValP) described in Section 1.2 : -In Section 3.1, we devise some bivariate small integral roots certificates that can be generated by using a variant of Algorithm 2. These certificates can then be verified by a checker that does not need to iterate Hensel lifting all over again, and whose correctness proof relies on Corollary 2. -Then in Section 3.2, we devise some ISValP certificates that are built upon the previous certificates, and that can be generated using the modular, bivariate version of Coppersmith's technique, based on the LLL algorithm [START_REF] Lenstra | Factoring polynomials with rational coefficients[END_REF]. Again, the checker does not need any call to the LLL algorithm to verify the certificates.

3.1 Certificates for the bivariate small-integral-roots problem

The problem at stake

The problem we tackle in this section is to formally verify that we know all the small integral roots (x, y) ∈ Z × Z of a pair of bivariate polynomials (v 1 , v 2 ) on Z, with respect to bounds A and B: Problem 2 (Certifying bivariate small integral roots) Given

             v 1 ∈ Z[X, Y ], v 2 ∈ Z[X, Y ], A ∈ N, B ∈ N, S ⊂ Z × Z, (19)
we want to formally verify that the following property holds:

∀(x, y) ∈ Z 2 , (x, y) ∈ S ⇐⇒ v 1 (x, y) = 0 = v 2 (x, y) ∧ |x| A ∧ |y| B. ( 20
)
For this purpose, we follow the following steps:

1. Define the type for bivariate small integral roots certificates as a Record; 2. Define the checker for these certificates as a Coq function; 3. Prove the correctness of this checker (if it returns true, then (20) holds), using the uniqueness result given by Corollary 2; 4. Optimize the implementation of the checker by using efficient data structures for the computation.

We detail below the main definitions or statements that are involved in the three first steps, while the fourth step will be described in Section 4.

The type of certificates

First, we need to define the type of certificates for the bivariate small-integral-roots problem. Typically, it should gather at the same time the "input/output" of the problem-here given in Equation ( 19)-, as well as some additional data that allow one to deduce that the particular instance of the problem is valid-here it amounts to showing that Equation ( 20) holds. We thus consider the following type:

Record bivCert := BivCert { bc_v1 : {bipoly int} ; bc_v2 : {bipoly int} ; bc_A : nat ; bc_B : nat ; bc_p : nat ; bc_k : nat ; bc_L : list (int * int * bool) }.

In addition to the polynomials v 1 , v 2 , and the bounds A and B-which constitute the "input" of the problem-, the certificate stores the prime p and the final value of k that both occur in Algorithm 2.

We recall that p is a prime number that fulfills [START_REF] Boneh | Finding Smooth Integers in Short Intervals Using CRT Decoding[END_REF], and k is the number of iterations of Hensel lifting performed to reach the specified bounds (i.e., we have p 2 k > 2 • max(A, B)). Finally, the list L involved in the certificate is supposed to store all the modulo-(p 2 k ) roots4 of (v 1 , v 2 ), along with a Boolean value that indicates whether each modular root is an actual integral root or not.

Example 2 (Bivariate small integral roots certificate) A possible certificate built on

Example 1 is the following one: where v 1 and v 2 are not shown here for brevity (see Example 1 for their definition).

C := v 1 , v 2

The certificate checker

Let us now present the algorithm to be implemented in the certificate checker. We say that a certificate (v 1 , v 2 , A, B, p, k, L) is valid if the following conditions are fulfilled:

p is a prime number, (21a)
if the list L is not empty, the integer k 0 satisfies

p 2 k > 2 • max(A, B), ( 21b
)
for L p := {(r mod p, s mod p) | ∃t ∈ {true, false} , (r, s, t) ∈ L} , we have:

∀(x, y) ∈ 0, p -1 2 , (x, y) ∈ L p ⇐⇒ v 1 (x, y) ≡ 0 ≡ v 2 (x, y) (mod p), (21c) the elements of L p are pairwise distinct, (21d)
∀(x, y) ∈ L p , det J v1,v2 (x, y) ≡ 0 (mod p), (21e)
and for all (r, s, t) ∈ L, we have:

v 1 (r, s) ≡ 0 ≡ v 2 (r, s) (mod p 2 k ), (21f)
|2 • r| p 2 k and |2 • s| p 2 k , (21g) t = true ⇐⇒ v 1 (r, s) = 0 = v 2 (r, s) ∧ |r| A ∧ |s| B . ( 21h
)
These conditions are implemented by a Coq function biv_check that takes a certificate of type bivCert and returns a boolean that indicates whether the certificate is valid or not.

The formal correctness proof

The correctness proof of the bivariate checker biv_check consists in proving that any certificate that is accepted by the checker is valid, i.e., contains all the small integral roots of the considered pair of bivariate polynomials. We thus prove the following Theorem 2 (Correctness of biv_check) For any bivariate small-integral-roots certificate C := (v 1 , v 2 , A, B, p, k, L), if (biv_check C) returns true, then for all (x, y) ∈ Z 2 we have the equivalence

v 1 (x, y) = 0 = v 2 (x, y) ∧ |x| A ∧ |y| B ⇐⇒ (x, y, true) ∈ L. (22)
Remark 1 (Invertibility hypothesis) Equation (21e) can be viewed as a "distributed version" of the main invertibility hypothesis [START_REF] Boneh | Finding Smooth Integers in Short Intervals Using CRT Decoding[END_REF] of bivariate Hensel's lemma. Indeed, according to (21c), the values that are given to variable (x, y) in (21e) correspond to all the modulo-p roots of (v 1 , v 2 ), so that (21c) and (21e) imply

∀(x, y) ∈ 0, p -1 2 , v 1 (x, y) ≡ 0 ≡ v 2 (x, y) (mod p) =⇒ det J v1,v2 (x, y) ≡ 0 (mod p).
Let us now present the main steps of this correctness proof, which has been mechanized in Coq.

Proof (Theorem 2) Suppose that (biv_check C) returns true for a given certificate C := (v 1 , v 2 , A, B, p, k, L), and let (x, y) ∈ Z 2 . Note that modular equalities over Z 2 below are taken coordinatewise. -Let us prove the "=⇒" part of equivalence [START_REF] Karatsuba | Multiplication of Many-Digital Numbers by Automatic Computers[END_REF], relying on Corollary 2. To begin with, the assumption v 1 (x, y) = 0 = v 2 (x, y) implies

v 1 (x mod p, y mod p) ≡ 0 ≡ v 2 (x mod p, y mod p) (mod p). ( 23
)
So applying (21c) with (x 0 , y 0 ) := (x mod p, y mod p) leads to (x 0 , y 0 ) ∈ L p , that is by definition of L p , ∃(r, s, t) ∈ L, (x 0 , y 0 ) = (r mod p, s mod p), hence (x, y) ≡ (r, s) (mod p).

(24) Now let us show that (r, s, t) = (x, y, true), which will allow one to deduce that (x, y, true) ∈ L.

-First, we want to show that (x, y) = (r, s). We can use Corollary 2: the main hypothesis ( 7) is fulfilled thanks to (21e), as explained in Remark 1. Moreover, we immediately have

v 1 (x, y) ≡ 0 ≡ v 2 (x, y) (mod p 2 k )
as well as v 1 (r, s) ≡ 0 ≡ v 2 (r, s) (mod p 2 k ) by using (21f). So we can apply Corollary 2 and use [START_REF] Lammich | Automatic Data Refinement[END_REF] to deduce that (x, y) ≡ (r, s) (mod p 2 k ).

In addition, with K := p 2 k , we have

         |x| A < K/2 by (21b), |r| K/2 by (21g), |y| B < K/2 by (21b), |s| K/2 by (21g).
Hence by applying twice Corollary 1, we get x = r and y = s.

-Second, we want to show that t = true, knowing that the equivalence (21h) holds. By hypothesis, we have

v 1 (x, y) = 0 = v 2 (x, y) ∧ |x| A ∧ |y| B.
It then suffices to combine this result with the fact that (x, y) = (r, s) that we have proved just now to get v 1 (r, s) = 0 = v 2 (r, s) ∧ |r| A ∧ |s| B, which is equivalent to the desired result t = true. -For the "⇐=" part of [START_REF] Karatsuba | Multiplication of Many-Digital Numbers by Automatic Computers[END_REF], we need to verify that the values stored in L are actual integral roots. This directly follows from (21h). ⊓ ⊔

Certificates for the Integer Small

Value Problem (ISValP)

The problem at stake

In the previous Section 3.1, we have seen that we can take advantage of Hensel's lemma to devise some bivariate small integral roots certificates. We will now see how this kind of certificates can be reused to address the ISValP problem. For this, we introduce the list S and rephrase Problem 1 as:

Problem 3 (Certifying ISValP) Given              P ∈ Z[X], A ∈ N, B ∈ N, M ∈ N, S ⊂ Z × Z, (25)
we want to formally verify that the following property holds:

∀(x, y) ∈ Z 2 , |x| A ∧ |y| B ∧ P (y) ≡ x (mod M ) =⇒ (x, y) ∈ S. ( 26
)
Notice that the order of variables x and y in ( 26) is reversed with respect to Problem 1 and [START_REF] Stehlé | Searching Worst Cases of a One-Variable Function Using Lattice Reduction[END_REF]. This different wording is due to some Coq technicalities, on which we will elaborate in Section 5. Note also that unlike Problem 2 where the required property [START_REF] Haftmann | Data Refinement in Isabelle/HOL[END_REF] was an equivalence, here we only require an implication [START_REF] Martin-Dorel | Contributions to the Formal Verification of Arithmetic Algorithms[END_REF] saying that S is a superset of the ISValP solutions. Indeed we are especially interested in formally proving that no ISValP solution has been forgotten. And for computing the exact ISValP solutions, it suffices to check the ISValP conditions on each solution candidate (i.e., on each element of S).

Focusing on the SLZ algorithm

In order to understand which kind of certificates could be suitable for ISValP, we have to give a closer look to what the SLZ algorithm is doing. The main steps are summarized in Figure 1 and explained below.

Integer Small Value Problem

P ∈ Z[X], find all y ∈ -B, B such that |P (y) cmod M | A Q(X, Y ) := P (Y ) -X ∈ Z[X, Y ]
Bivariate Small Modular Roots Find all (x, y) ∈ -A, A × -B, B s.t. Q(x, y) ≡ 0 (mod M )

Coppersmith's technique with parameter α > 0:

Consider Qi,j (X, Y ) = M α-i Q(X, Y ) i Y j (i α).
Heuristically, find two Z-linear combinations v1, v2 of (Qi,j ) s.t.

∀x, y ∈ Z, |x| A and |y| B =⇒ |v l (x, y)| < M α . Notice that the small modulo-M roots of Q also satisfy v l (x, y) ≡ 0 (mod M α ). The algorithm starts from an ISValP instance and reduces it to a bivariate small integral root problem (to be solved by using Hensel lifting, as described in Section 3.1). Here is the proof of correctness of this reduction:

Bivariate Small Integral Roots

Proof (Reduction from an ISValP instance) Assume that y ∈ -B, B satisfies |P (y) cmod M |

A. According to Equation ( 6) in Section 2.1, this amounts to saying that we have (x, y) ∈ -A, A × -B, B such that

P (y) ≡ x (mod M ). Then, pose Q(X, Y ) := P (Y ) -X ∈ Z[X, Y ], implying that Q(x, y) ≡ 0 (mod M ). (27)
Then, choose an integer α > 0 (the Coppersmith parameter) and consider the family of polynomials

Q i,j (X, Y ) = Q(X, Y ) i M α-i Y j (i α).
Under heuristic assumptions, we may find two Z-linear

combinations 5 v 1 , v 2 of (Q i,j ) such that ∀(z, t) ∈ -A, A × -B, B , |v 1 (z, t)| < M α and |v 2 (z, t)| < M α . ( 28
)
By definition of (Q i,j ), v 1 and v 2 , Equation ( 27) implies

v 1 (x, y) ≡ 0 ≡ v 2 (x, y) (mod M α ). ( 29
)
Combining ( 28) and ( 29) leads to

v 1 (x, y) = 0 = v 2 (x, y). ⊓ ⊔
Note that the small integral roots of the reduced problem may happen to be a strict superset of the solutions of ISValP.

Designing the certificate

The certificate must obviously contain the initial values P , A, B, M . We also need to include the Coppersmith parameter α. As regards the linear combinations v 1 , v 2 , we focus on their coefficients with respect to the family of polynomials {Q i,j (X, Y )} i,j , and we store these coefficients as mere polynomials u 1 , u 2 in the monomial basis. In other words, for both l ∈ {1, 2}, the certificate will contains the polynomial

u l = i α j u l;i,j X i Y j .
which will represent the linear combination

i α j u l;i,j Q i,j (X, Y ) = v l .
We where

P (X) = 36506148256923413X 3 + 28709603185988793532416X 2 + 15052100435175692583523319808X
+ 3945817816478696756615137147748352, and

u 1 (X, Y ) = X 2 + (-6Y -1593703)X + (9Y 2 + 4781109Y + 634972313052), u 2 (X, Y ) = -25X + (76Y + 19921289).

The certificate checker

The algorithm of our ISValP certificate checker is described in Figure 2.

check_ISValP (P, A, B, M, α, u 1 , u 2 , p, k, L) ∈ {true, false} 1. Compute Q(X,Y ) := P (Y ) -X ∈ Z[X, Y ]; 2. Compute v l (X, Y ) = M α • u l Q(X, Y ) M , Y for l ∈ {1, 2};
3. Return true if and only if:

           i,j |v 1;i,j | A i B j < M α , i,j |v 2;i,j | A i B j < M α , (v 1 , v 2 , A, B, p, k, L
) is a valid SIntRootP certificate (as specified in Section 3.1.3). In this algorithm, Line 1 computes the bivariate polynomial Q, then Line 2 performs a change of polynomial basis: the polynomials v 1 (X, Y ) and v 2 (X, Y ) are initially expressed in the polynomial basis {Q i,j (X, Y )} i,j and we need to compute their coefficients in the usual monomial basis. In the formalization, we perform this change of polynomial basis in two steps:

1. First, we "lift" the coefficients of u = i α j u ij X i Y j from u ij to u ij M α-i ; 2. Then in the polynomial i α j u ij M α-i X i Y j , we replace the indeterminate X with the polynomial Q to get v = i α j u ij M α-i Q(X, Y ) i Y j ;
this may be accomplished by using a Horner-based polynomial composition.

Weighted Norm-1. Now let us explain Line 3 in Figure 2. A crucial ingredient for building our ISValP checker consists in formalizing a simple and efficient test to ensure that (28) is fulfilled. We introduce the notion of "weighted norm-1," defined for any P ∈ Z[X, Y ] and (A, B) ∈ N 2 by:

|P |(A, B) := i,j |P i,j |A i B j ∈ N.
This notion enjoys the following result:

Lemma 3 For any P ∈ Z[X, Y ] and any (A, B) ∈ N 2 , we have ∀x, y ∈ Z, |x| A ∧ |y| B =⇒ |P (x, y)| |P |(A, B).
Thanks to this result (and considering P := v l ), we can notice that it is sufficient to have |v l |(A, B) < M α to make sure that Equation (28) holds. (Note that the design of the SLZ algorithm internally relies on a similar notion, so the fact that it only provides a sufficient condition does not raise any issue.)

The formal correctness proof

The correctness claim for our ISValP checker is given by the following statement:

Theorem 3 (Correctness of check_ISValP) For any ISValP certificate C := (P, A, B, M, α, u 1 , u 2 , p, k, L) that is accepted, we have the following property:

∀(x, y) ∈ -A, A × -B, B , P (y) ≡ x (mod M ) =⇒ (x, y, true) ∈ L.
As pointed out in Section 3.2.2, Coppersmith's technique proceeds by implication, so it finally leads to an implication in Theorem 3. In other words, all the solutions of the ISValP instance are gathered in the list L. At the very end (to compute the hardest-to-round points of the considered elementary function), it will suffice to check the solution candidates gathered in the list L.

In order to prove Theorem 3, we closely follow the reduction proof given in Section 3.2.2. We formally prove Lemma 3, then we combine it with Lemma 1 to obtain the following result: Lemma 4 For any P ∈ Z[X, Y ] and any (N, A, B) ∈ Z 3 , we have:

∀x, y ∈ Z, |x| A ∧ |y| B ∧ |P | (A, B) < N ∧ P (x, y) ≡ 0 (mod N ) =⇒ P (x, y) = 0.
We use this result twice, one for each polynomial v 1 , v 2 . Then we apply Theorem 2 to get Theorem 3.

Effective certificate checkers

The certificate checkers that we have presented in Sections 3.1 and 3.2 have been developed using some data structures that are suitable for proving, but not for computing. To get effective computation, we need to re-implement our certificate checkers with effective data structures. This is a standard approach, see e.g. [START_REF] Dénès | A Refinement-Based Approach to Computational Algebra in Coq[END_REF].

Focusing on modularity and efficiency

We implement the effective version of our certificate checkers in a modular way, relying on the Module system of Coq. These checkers are then proved correct with respect to the reference checkers presented in Section 3. A key aspect of our modular implementation is that we can easily change the arithmetic with which these checkers compute. Three instantiations are available:

binary integers from Coq's standard library ZArith that represent numbers as list of booleans with linear access; machine-efficient integers from Coq's standard library BigZ that represents numbers as complete binary trees of machine integers with logarithmic access; -Integers Plus Positive Exponent (IPPE) numbers that are unevaluated dyadic numbers m×2 e with e 0. We detail below the motivation and the formalization of this latter implementation.

This last instantiation has been motivated by the fact that the coefficients of the approximation polynomials that typically occur in the final verification chain for SLZ are floating-point numbers in radix-2. In order to switch to integer arithmetic, these polynomials are then scaled by a large power of the radix to generate an instance of ISValP. For instance, in some example presented in Section 4.2.1, the power of 2 to be considered is 2 13660 , leading to a polynomial on Z whose degree-0 coefficient is M × 2 10629 where M is an odd integer in 2 3032 , 2 3033 -1 . Factoring out the powers of 2 in all computation seems a reasonable optimization. In fact, it leads to an appreciable twofold speedup as we will see in Section 4.2.1. Some obvious optimizations were mandatory in order to get reasonable checking time:

-First, the quantity det J v1,v2 (x, y) involved in Condition (21e) can be computed in several ways: we can either compute the determinant

∆(X, Y ) = ∂ v1 ∂X × ∂ v2 ∂Y - ∂ v2 ∂X × ∂ v1
∂Y and deduce the quantity ∆(x, y), or perform the Horner evaluation beforehand on each of the derivatives and compute the determinant of the resulting matrix. The latter calculation has the benefit of not requiring multiplication of polynomials, since the determinant acts on a matrix in M 2 (Z) instead of M 2 (Z[X, Y ]). For implementing this calculation, we redefine the type of order-2 matrices as a polymorphic record with four projections, which is better suited to effective computation than the SSReflect type 'M_2 for order-2 matrices. Then we refine the corresponding algorithms with respect to 'M_2.

-Second, the case where L is empty is very frequent so worth optimizing. In particular, when this list is empty, some computations like that of the Jacobian determinant above are unnecessary.

-The last optimization concerns Condition (21c) (and (21e) as well). It consists in checking whether the list Since the coefficients of v 1 and v 2 are huge, the polynomials v 1 (X, Y ) (mod p) and v 2 (X, Y ) (mod p) are computed once before performing the different Horner evaluations.

L p = {(

Evaluation of the performances of the ISValP checker

The current section is organized as follows. First, we evaluate the performances of the current version of our library by focusing on four ISValP instances that were first presented in [26, Section 5.4.5] and which cover a broad spectrum of situations related to our target application (solving the TMD using SLZ). Second, we present a complete verification of an ISValP computation that addresses the full range of an exponent, in order to demonstrate the scalability of our Coq implementation.

The certificates were generated using SAGE (which relies, mostly, on the MPIR, FLINT and fpLLL libraries). All the Coq benchmarks have been computed on an Intel R Xeon R E5-2667 CPU clocked at 2.90 GHz, using OCaml 4.00.1 and the following two versions of Coq: the experimental native-coq6 branch of Coq [START_REF] Boespflug | Full Reduction at Full Throttle[END_REF] and the current stable version of Coq (version 8.4pl4). Note that the native-coq version is obviously faster but it has also a larger trusted computing base since it dynamically links new code generated by the OCaml compiler.

Four typical instances of ISValP

For each instance of ISValP, Table 2 gives the function at stake, the precision variables n and n ′ that were presented in Section 1.2, the normalized (or scaled) approximation polynomial P ∈ Z[X], the modulo M , and the bounds A and B. These instances have been chosen so that -#1, #2, and #3 deal with the binary64 format, while #4 deals with binary128; -#1 contains a hardest-to-round point, unlike other instances; -#3 and #4 deal with the Approximate TMD (cf. Section 1.2), the value of n ′ being much larger than 2 × n. For a given parameter α, SLZ produces the polynomials v 1 and v 2 , and runs the bivariate Hensel lifting to generate the ISValP certificates. Table 3 summarizes the data stored in the certificates, and the timings for their verification inside Coq. The last four columns give the timings using the compiled native_compute reduction tactic of native-coq (using IPPE numbers, then BigZ numbers), and the timings using the interpreted vm_compute reduction tactic of Coq 8.4pl4 In what follows, we focus on the last two rows of Table 3 (as using native-coq yields no appreciable speedup for the small examples #1 and #2, because of the compilation overhead). These figures indicate that using IPPE integers in place of the standard BigZ integers yields a 1.7x to 4x speedup, and using the compiled native_compute in place of the interpreted vm_compute of Coq 8.4pl4 yields a 4.5x to 22x speedup. Also, the fastest version of the Coq checker is only 6 to 12 times slower than Maple. This is reasonable if we consider that the formal checker is executed in a trusted framework with restricted computing power, unlike Maple that notably takes advantage of the GMP library.

Large-scale benchmarks for our formalized ISValP checker

We now focus on sets of certificates that are needed to address the full range of an exponent in binary64 (formerly known as the double precision format). We recall that the bound B involved in the ISValP certificates directly indicates the number of floating-point numbers that are taken into account. As a result, we would need 2 53-1 /B = 2 20 certificates similar to the ISValP instance #3 of Table 3 to address all the floating-point numbers 7 of one exponent of binary64. We can then extrapolate the overall timing that would be required in this context: 0.154 s × 2 20 ≈ 45 hours. And the total size of the corresponding certificates would be of at least 2 806 bytes×2 20 ≈ 2.9 GB of ASCII to be fed to Coq. As a matter of fact, we can do much better by taking a larger value of α for the SLZ algorithm. For the parameters given in 5. We can see that using the native_compute reduction tactic in place of vm_compute yields a speedup of 5x on average. Almost all certificates contain no 7 Technically speaking, we only consider the normal binary floating-point numbers here, i.e., their leading bit is assumed to be 1. root and are relatively easy to verify. Also, there is a ratio of 6.6x (= 140 s/21.2 s) between the generation using SAGE and the verification in Coq. This means that using certificates, the computation has been formally verified with an overhead of 15 % only. As shown by our implementation of the checker in Maple, progress could be made to further reduce this overhead. Note that for the very big certificate #4 of Table 3, the formalized ISValP checker is only one order of magnitude slower than Maple.

f n n ′ range B α

Discussion

In order to prove the correctness of our checkers, a variety of mathematical notions has to be formalized: Taylor's theorem for bivariate polynomials, Jacobian matrix for pairs of bivariate polynomials, Cramer's rule for order-n matrices in modular arithmetic, and weighted norm-1 for bivariate polynomials among others.

Taylor's theorem already existed in Coq for real analysis. It has been reproved independently so as to be applicable to polynomials over an arbitrary ring (this is integrated in the standard distribution of SSReflect). This leads to the notion of normalized derivatives that uses binomial coefficients to internalize the division by 1 n! that is not available on a ring. Doing so, we also avoid a dependency with the real analysis and its axiomatic approach. A similar dependency could also occur for the IPPE integers, since they rely on the Coq.Interval library. We carefully checked that no dependency was introduced also here. So, our development is axiom-free.

A careful reader would have noticed that in Section 3.2, we have interchanged the two variables X and Y with respect to [35, Section 3.2], which was considering Q i,j (X, Y ) = X i Q j (X, Y )M α-j (j α). This technical detail comes from the way we have represented bivariate polynomials:

Z[X, Y ] is encoded as (Z[Y ])[X].
A consequence of the asymmetry between the two variables X and Y is that we get for free the substitution of a polynomial Q(X, Y ) for the variable X in P (X, Y ) with the simple univariate composition P • Q. The other substitution with the variable Y is less direct. The ISValP checker is then encoded using the first substitution.

Formalizing effective implementations of algorithms in a proof assistant is often a challenging task. A key ingredient that we have been using here is to develop modular code. The Coq proof assistant provides three mechanisms for modularity: type classes [START_REF] Sozeau | First-Class Type Classes[END_REF], canonical structures [START_REF] Saïbi | Typing Algorithm in Type Theory with Inheritance[END_REF], and modules [START_REF] Chrząszcz | Modules in Coq Are and Will Be Correct[END_REF][START_REF] Chrząszcz | Implementing Modules in the Coq System[END_REF]. Modules are less generic than the other two (which are first-class citizens) but they have a better computational behavior. Indeed, module applications are performed statically, so the code that is executed is often more compact. We have been using modules to develop the generic implementation presented in Section 4.1. In particular, this has clearly simplified the formalization of IPPE integers, since the Coq.Interval library is also relying on modules.

An aspect of our work that could be largely automated is the refinement from the naive checker that works on a simple data-structure to the more efficient one. For the moment, this is done in an ad-hoc way using morphisms between datastructures and explicitly applying rewriting rules for the correctness proof. The use of a more systematic approach such as the ones proposed in [START_REF] Cohen | Refinements for Free! In[END_REF][START_REF] Haftmann | Data Refinement in Isabelle/HOL[END_REF][START_REF] Lammich | Automatic Data Refinement[END_REF] would be a clear improvement.

A key part of the mathematical background that underlies the formal verification of our ISValP checker is Hensel's lemma, in the bivariate case. We recall that ISValP can be solved by using a variant of Coppersmith's technique, which exists in two main versions, univariate and bivariate [START_REF] Coppersmith | Finding a Small Root of a Univariate Modular Equation[END_REF][START_REF] Coppersmith | Finding a Small Root of a Bivariate Integer Equation; Factoring with High Bits Known[END_REF]. We thus started our formalization effort by mechanizing the univariate case of Hensel's lemma and generalizing it to the bivariate case for the needs of ISValP. But beyond ISValP, further generalizations of Hensel's lemma could be envisioned, for instance to address the small integral roots problem for n bounds and n polynomials. To this end, we could rely on our formalization of Definition 4, which already handles order-n matrices. But a key building block that would probably require more work consists of the availability of a comprehensive library of multivariate polynomials.

Conclusion and perspectives

In this work, we have built several effective verified certificate checkers upon Hensel's lemma. In particular, our ISValP certificate checker is suitable to validate that all the solutions of given instances of the Integer Small Value Problem (ISValP) have been found. This problem consists of finding all the small integer entries on which a univariate polynomial P ∈ Z[X] has small values modulo a large integer M . It can be solved by using the SLZ algorithm, which is based on a variant of Coppersmith's technique [START_REF] Coppersmith | Small Solutions to Polynomial Equations, and Low Exponent RSA Vulnerabilities[END_REF]. We recall that this algorithmic chain consists of four main steps: (i) domain splitting, then on each sub-domain: (ii) polynomial approximation, (iii) the modular, bivariate version of Coppersmith's technique, and (iv) bivariate Hensel lifting. The formalization presented here constitutes a key formal component for certifying the SLZ algorithm. In this paper, we have devised some certificates that address steps (iii) and (iv) in a formal setting. To demonstrate the scalability of this component, we have evaluated its performances within the Coq proof assistant on realistic examples generated by the SLZ algorithm.

We have been using the Coq proof assistant to design certificates that make it possible to validate a long computation. We have derived some reference checkers that have been proved correct. In our case, they were efficient enough to be used for validation. It is always debatable how much trust one can have in a Coq computation. What can be said is that all the components that we have been using (starting from arbitrary precision arithmetic) have been designed with a very strong concern about correctness.

Hensel's lemma has already been formalized in an abstract setting in the Isabelle/HOL proof assistant [START_REF] Kobayashi | Formalization of Hensel's lemma[END_REF]. To the best of our knowledge, our work constitutes the first formalization of bivariate Hensel's uniqueness lemma along with verified, Hensel-based certificate checkers. Further optimizations of our certified library could be performed. For example, our formalization of quadratic Hensel lifting naturally yields an arithmetic modulo p 2 k , while we could restrict ourselves to the smallest power p m > 2 • max(A, B) in the final certificates. Even if this would have no impact on most certificates (those whose list L is empty), it would be interesting to measure the performance gain for the other few certificates. Finally we would like to improve the multiplication algorithm (and thereby, the bivariate polynomial composition one) by implementing a Karatsuba-based approach [START_REF] Karatsuba | Multiplication of Many-Digital Numbers by Automatic Computers[END_REF]. This may lead to a further speedup for verifying some "extreme" certificates such as our binary128 example in Section 4.2.1, while keeping the confidence we have in our formally verified ISValP certificate checker.

Our next priority will be to combine the CoqHensel machinery presented here with that of the CoqApprox library for certified polynomial approximation [START_REF] Brisebarre | Rigorous Polynomial Approximation Using Taylor Models in Coq[END_REF][START_REF] Martin-Dorel | Certified, Efficient and Sharp Univariate Taylor Models in COQ[END_REF]. We will then get a validation for the whole SLZ algorithm for solving the Table Maker's Dilemma. It will then be possible to provide certificates that formally guarantee that we know all the bad cases for correctly rounding a given mathematical function.

Finally the part of this work which, beyond Hensel's lemma, addresses Coppersmith's technique, might have, up to slight modifications, some use in applications of the latter, for instance the list decoding algorithm for CRT-type codes, see [START_REF] Boneh | Finding Smooth Integers in Short Intervals Using CRT Decoding[END_REF][START_REF] Daniel | Simplified High-Speed High-Distance List Decoding for Alternant Codes[END_REF], or for RS-codes already pointed out in Section 1.

  a, b) and |a| A and |b| B then S ← S ∪ {(a, b)} 14 end 15 return S

  , A := 2 87 , B := 2 4 , (p := 3) , (k := 6) , L := (233673222969575457782319, 11, true), (33647808062653569260574440, -14, true).

FindFig. 1

 1 Fig. 1 Solving ISValP using the SLZ algorithm.

C

  ′′ := P, A := 2 87 , B := 2 4 , M := 2 92 , (α := 2) , u 1 , u 2 , (p := 3) , (k := 6) , L := (233673222969575457782319, 11, true), (33647808062653569260574440, -14, true) ,

Fig. 2

 2 Fig. 2 ISValP certificate checker.

  r, s) mod p | ∃t ∈ {true, false} , (r, s, t) ∈ L} satisfies ∀(x, y) ∈ 0, p -1 , (x, y) ∈ L p ⇐⇒ v 1 (x, y) ≡ 0 ≡ v 2 (x, y) (mod p).

Table 1

 1 Modular roots obtained by iterating bivariate Hensel lifting (with a centered modulo) for initial values (0, 2) and (2, 1), which are modulo-3 roots of (v 1 , v 2 ). The variables k, r, s, a, b are those of Algorithm 2.2.3 Bivariate Hensel's lemma: pencil and paper proofWe start this section by introducing two elementary arithmetic results that are used several times in what follows.

	4)

Lemma 1 For all d ∈ N and z ∈ Z, if |z| < d and z ≡ 0 (mod d), then z = 0.

Corollary 1 For all (d, m, n) ∈ Z 3 , if |2 • m| d, |2 • n| < d,

and m ≡ n (mod d), then we have m = n. Proof (Corollary 1) Suppose d, m, n ∈ Z satisfy the hypotheses of Corollary 1. By the triangle inequality, we have |2 •

  then need to include the values p, k, and L that correspond to the bivariate small integral root certificate. Altogether this gives the following Coq record:

	Record cert_ISValP := Cert_ISValP
	{ ic_P : {poly int}
	; ic_A : nat
	; ic_B : nat
	; ic_M : nat
	; ic_alpha : nat
	; ic_u1 : {bipoly int}
	; ic_u2 : {bipoly int}
	; ic_p : nat
	; ic_k : nat
	; ic_L : list (int * int * bool)
	}.
	Let us give an example of ISValP certificate.

Example 3 (ISValP certificate) This example is based on the same numerical data as in Examples 1 and 2 :

Table 2

 2 Short description of four instances of ISValP.

	Instance of ISValP #1	f exp	n 53	n ′ 100	deg(P ) 2	max i (|P i |) 1.68×2 237	M 2 185	A 2 139	B 2 12
	#2	exp	53	100	2	1.22×2 237	2 185	2 139	2 12
	#3	exp	53	300	12	1.36×2 996	2 942	2 696	2 32
	#4	exp	113 3000	90	1.36×2 13661	2 13547	2 10661	2 72

Table 3

 3 (again using IPPE numbers, then BigZ numbers). Short description of several ISValP certificates with their verification timing in Coq. The first column refers to Table2.

	Inst.	α	M α	p	size(L)	Maple 16	CPU time to return true native-coq Coq 8.4pl4 IPPE BigZ IPPE BigZ
	#1	2	2 370	5	1	0.002 s	0.069 s	0.077 s	0.070 s	0.066 s
	#2	2	2 370	23	0	0.005 s	0.086 s	0.089 s	0.095 s	0.100 s
	#3	4	2 3768	5	0	0.024 s	0.154 s	0.261 s	0.705 s	1.641 s
	#4	6	2 81282	5	0	29.91 s	6 ′ 6 s	11 ′ 14 s	1h 3 ′ 39 s	4h 16 ′ 29 s

Table 4 ,

 4 only 2 53-1 /B = 212 certificates are necessary for addressing the same range. The generation of each certificate took approximately 140 s, which leads to around 20 h of real time on a 8-core Intel R Xeon R X5550 CPU clocked at 2.67 GHz. Roughly 95 % of this time was spent in the lattice reduction calls, on a 105 × 3056 matrix with 700-bit integer entries.All these certificates were verified on the same machine as in Section 4.2.1, by compiling 8 Coq files, each containing the evaluation (i.e., Eval native_compute or Eval vm_compute) of the ISValP checker on 512 certificates. The average verification time (with standard deviation) for the considered population of ISValP certificates is given in Table

Table 4

 4 Short description of the common values of 4096 ISValP certificates.

						total size of input files
	exp	53 300	[1, 2]	2 40	13	100 MB
	Class of certificates	Average CPU time ± std deviation to accept one certificate native-coq & native_compute Coq 8.4pl4 & vm_compute
	4090 certs. s.t. size(L) = 0		21.2 s ± 4.5 s	109.0 s ± 22.1 s
	6 certs. s.t. size(L) = 1		42.8 s ± 10.9 s	249.3 s ± 62.0 s
	all the 4096 certificates		21.2 s ± 4.6 s	109.2 s ± 22.9 s

Table 5

 5 Some statistics on the formal verification of 4096 ISValP certificates within Coq, using the "bigZ × bigN" implementation of IPPE integers. The biggest chosen prime is pmax := 47.

Namely, (0, 0) and (1, 1) are modulo-2 roots of both (v 1 , v 2 ) and det Jv 1 ,v 2 .

These modular roots being considered in an appropriate range, say inp 2 k 2 , p 2 k

These polynomials v 1 and v 2 can be found by means of the LLL algorithm for Euclidean lattice basis reduction, which is the most time consuming part of the overall algorithm SLZ. More details on the way the polynomials v 1 and v 2 are computed may be found in[START_REF] Stehlé | Algorithmique de la réduction des réseaux et application à la recherche de pires cas pour l'arrondi des fonctions mathématiques[END_REF][START_REF] Stehlé | Searching Worst Cases of a One-Variable Function Using Lattice Reduction[END_REF][START_REF] Stehlé | On the Randomness of Bits Generated by Sufficiently Smooth Functions[END_REF].

URL: https://github.com/maximedenes/native-coq

Acknowledgements This research was partly funded by the TaMaDi project of the Agence Nationale de la Recherche (ref. ANR-2010-BLAN-0203-01). It was mainly done while the first author was PhD student in École Normale Supérieure de Lyon, with the LIP research laboratory. The authors are very grateful to the anonymous reviewers, whose suggestions have been very helpful for revising this paper.

A One additional mathematical proof

The proposition below illustrates the usefulness of Lemma 2 beyond the formal proof of the certificate checkers presented in the paper.

Proposition 1 (Correctness of Algorithm 2) For any

if there exists a prime p satisfying [START_REF] Boneh | Finding Smooth Integers in Short Intervals Using CRT Decoding[END_REF], then the output of Algorithm 2 satisfies

Proof (Proposition 1) First, let us notate

In order to prove that S = R, we proceed by double-inclusion:

-First, we have S ⊂ R since S is initially empty, and each element added to S by Algorithm 2 on Line 13 does belong to R;

So the foreach loop of Algorithm 2 is executed for (r, s) := (x mod p, y mod p). In particular, it iterates bivariate Hensel lifting on this modulo-p root of (v 1 , v 2 ), and the corresponding while loop stops as soon as we have K > 2 • max(A, B). Let us denote by ℓ the final value of k (at the end of the while loop). Then notice that the successive values of variables (a, b) in this loop correspond to a sequence (a k , b k ) 0 k ℓ that satisfies relations ( 9) and [START_REF] Chrząszcz | Modules in Coq Are and Will Be Correct[END_REF]. Given that (7) holds by hypothesis, and that ( 8) is trivially satisfied by the integral root (x, y), we can apply Lemma 2. To sum up, we have

Then we deduce that x = a ℓ and y = b ℓ by applying twice Corollary 1. Thence, the small integral root (x, y) of (v 1 , v 2 ) has indeed been found and added to S by Algorithm 2. ⊓ ⊔


virgile. Virgile Lacharnay
email: lacharnay@lurpa.ens-cachan.fr

Christophe Tournier
email: christophe.tournier@lurpa.ens-cachan.fr

Claire Lartigue
email: claire.lartigue@lurpa.ens-cachan.fr

A potential field approach for collision avoidance in 5-axis milling

Keywords: 5-axis machining, Collision detection, Potential field, tool paths optimization

Introduction

5-axis surface machining is an essential process in the field of aerospace, molds and dies industries. 5-axis milling is required for the realization of difficult parts such as blades and impellers and is also very convenientto improve quality for the machining of deep molds in plastic injection and casting by reducing tool length. Despite the evolution of CAM software, 5-axis tool path programming requires advanced skills and collision detection remains a challenge during the tool path computation. One can distinguish two kinds of tool collision when addressing machining issues. The first one is local and only involves the active part of the tool, the second one is global and in that case the whole body of the tool, the tool holder and even the spindle can be considered.

Literature shows that collisions are often avoided using a geometric point of view. Indeed, the study traditionally starts with a representation of the tool geometry [START_REF] Bi | Generating collision-free tool orientation for 5-axis NC Machining[END_REF][2], then by a geometrical representation of the environment [3][4], thirdly by a collision test between the obstacle and the tool [START_REF] Monies | Determining a collision-free domain for the tool in five-axis machining[END_REF][6] and finally by a correction and optimization of the tool axis direction to avoid the obstacle [START_REF] Ho | Five axis tool orientation smoothing using quaternion interpolation algorithm[END_REF]. In the field of robotics, other methods based on potential fields virtually attached to the obstacle are developed. These fields emit a repulsive force on the robot when it enters the neighborhood of the obstacle [START_REF] Ge | New Potential Functions for Mobile Robot Path Planning[END_REF] [START_REF] Khatib | Real-Time Obstacle Avoidance for Manipulators and Mobile Robots[END_REF]. The aim of this paper is to show the benefit of using potential fields in order to prevent collisions during the computation of the tool axis orientation along a given tool path followed by the center of a ball-end milling tool.

Main idea

In the proposed approach, 5-axis collision avoidance is managed thanks to repulsive forces deriving from a potential field. The tool is considered as a rigid body moving in 3D space on which repulsive and attractive forces are acting. Repulsive forces are due to potential fields attached to check surfaces and an attractive force exerted by a spring is introduced to restore the tool axis orientation in the programmed configuration. Moreover, a viscous damper is used to allow the system returns to steady state without oscillating. The tool center follows the programmed path on the part surface whereas the tool axis orientation is modified to avoid the obstacle by resolving the fundamental principle of dynamics. Differential equations of the tool motion are solved using an Ordinary Differential Equation solver (ODE).

The repulsive force comes from robotics [START_REF] Khatib | Real-Time Obstacle Avoidance for Manipulators and Mobile Robots[END_REF] and is formulated as follows: with  η: Repulsive force coefficient  r: Distance between obstacle and tool gravity center  r 0 : Programmed neighborhood for the obstacle Two different applications are proposed. First, a tool moving along a straight line enters successively two spherical potential fields created by two single points. The tool axis orientation is automatically modified to avoid the collision and go back to its programmed orientation as fast as possible with a smooth response. This first case allows tuning simulation parameters such as mass, spring constant and damping as well as tool velocity during simulationto avoid oscillations.



r  r 0  r  r 0  F rep = h.( 1 r - 1 r 0 ) 2 . 1 r 2 .grad    (r)  F rep =  0 ì í ï î ï
The second study deals with a 5-axis open pocket application such as impeller machining. Part surface and check surfaces are modeled as Bezier patches. The check surfacesare discretized and every point of the mesh is considered as a repulsive point. The main purpose is to avoid oscillations by choosing an optimum mesh size and by programming the return in the initial orientation with a critical or over damped behavior. Fig. 1 and Fig. 2 illustrate the results of both studied cases.

Conclusions

Research works introduced in this paper show an original method for 5-axis collision avoidance by using potential fields commonly used in robotics. Using potential fields allows to be sure that no collision will ever occur with the obstacles as the neighborhood emits a repulsive force growing to infinity when the tool gets closer to the obstacle. As the tool behaves like a damped harmonic oscillator, investigations are carried out to prevent the tool axis orientation from oscillating. Numerical investigations show that the proposed approach is efficient and perspectives brought by these works are numerous. In particular, the reduction of the tool length to minimize tool deflection and the smoothness of the tool axis orientation variation to ensure High Speed Machining are under investigations.

Fig. 1 :

 1 Fig. 1: Avoidance of 2 spheres neighborhood Fig. 2: Avoidance of anImpeller-type surface