Frédéric Chyzak

Guillaume Melquiond

Éric Schost

Pierre-Yves Strub

Georges Gonthier

Santiago Zanella-Beguelin

Sylvie Boldo

Laurent Théry

Pierre- Yves Strub

Cyril

Hélène Philippe Damien

Pierre

Laurent

Marc, Rémi, Rémy, Suzy Louis

Daniel Hirschkoff

Guillaume Hanrot

Nicolas Brisebarre

Mathieu Sozeau

Antoine P Meriem

Natalia Fosca

Juliette Mathieu

. (aka Doudou Antoine A Fanny

Jean-Florent Il

Déborah

Adrien Katia

Jessica Laure

Alix Marion

Héloïse J'ai

Merci À Toute

La Famille

Mamie Suze

Papi Jean

A Complexity of Some Matrix Multiplication Algorithms

fait le chemin de la thèse en même temps que moi pour l'avoir rendue plus agréable par leur compagnie, je pense en particulier à Sonia, Salomé, Loïc et Alice.

Je tiens à remercier tous ceux qui rentrent moins facilement dans une case mais qui tiennent une place dans mon coeur depuis des années,

Context

Computers and computer science have progressively made themselves indispensable to mathematicians. This is of course true of communication with the ubiquitous use of search engines, of L A T E X [START_REF] Lamport | LATEX: a document preparation system: user's guide and reference manual[END_REF] to write papers, of Math Overflow [START_REF] Spolsky | Math Overflow[END_REF] to ask and answer questions from peers, of blogs [START_REF] Gowers | Massively collaborative mathematics[END_REF] [START_REF] Polymath | The "bounded gaps between primes" polymath projecta retrospective[END_REF] to discuss and produce new results.

More directly, the vast computational power available today has allowed mathematicians to extend their ability to do mental computations: they can experiment further than before, try to find more complicated counter-examples to conjectures, visualize patterns on more terms of a sequence. In An Introduction to Mathematics [START_REF] North | An introduction to mathematics[END_REF], Whitehead claims about mathematical symbolism that "civilization advances by extending the number of important operations which we can perform without thinking about them". Specialized software (for numerical computations, computer algebra, etc...) have arguably participated in this extension.

Not only have computations at the service of mathematicians become more ambitious, some recent developments have even attempted to cut out the middleman in the mathematician's main product: proofs. The Boolean Satisfiability Problem (SAT) consists in determining, given a formula φ involving boolean variables x 1 ,• • • , x n and the and (∧) and or (∨) connectives, whether there exist boolean values b 1 ,• • • , b n such that when substituting b i for x i , φ evaluates to true. For example, x 1 ∧(¬x 2 ∨ ¬x 1 ) is satisfiable with the assignment x 1 := true and x 2 := false. Thus, a successful assignment of boolean values to the variables of φ constitutes, up to a routine verification which can be accomplished by a simple program, a proof of the satisfiability of φ. Powerful solvers for the SAT problem (called SAT-solvers) have been developed. Surprisingly, complex problems can be reduced to the satisfiability of a (potentially large) boolean formula. The Pythagorean Triples Problem asks whether for all partitions of the set {1, 2,• • • , n} into 2 parts, at least one partition contains a Pythagorean triple. A recent paper [START_REF] Marijn | Solving and verifying the boolean pythagorean triples problem via cube-and-conquer[END_REF] shows that it holds using a SAT solver, producing a 200terabyte file. Finding proofs based on this method raises at least two questions. A first one is that of whether we can trust the proof checker as a program. This is a small instance of an issue which will be addressed in the next paragraph.

Secondly, the boolean encoding changes the statement of the problem so that if a proof is found, it does not follow the structure of the original problem. Thus, the proof may carry little to no insight as to why the theorem is true. About an unintelligible 1500 page computer-generated computer proof [START_REF] Phillips | Using automated theorem provers in nonassociative algebra[END_REF], Hales argues that "eventually, we will have to content ourselves with fables that approximate the content of a computer proof in terms that humans can comprehend" [START_REF] Thomas | Mathematics in the age of the Turing machine[END_REF].

If we stay on the side of proofs made and checked by humans, many mathematical proofs in the past century have come to rely decisively on computations; we list several examples from various domains of mathematics. The Four Color Theorem, stating that every planar graph can be colored with at most four colors so that no two adjacent nodes bear the same color, was first proved in 1977 [START_REF] Appel | Every planar map is four colorable. part i: Discharging[END_REF] by Appel and Haken by examining 1, 936 possible configurations using a computer. The Double Bubble conjecture, characterizing the shape which encloses two equal volumes with minimal surface as being made of two pieces of spheres separated by a flat disk, was first proved in 1995 by Hass, Hutchings and Schlafly. They reduced the proof to showing thousands of inequalities on simple integrals and wrote a custom C++ program doing interval arithmetic which ran in 10 seconds at the time. Tucker solved Smale's 14th problem in 1999, proving that the Lorenz attractor is strange [120], notably by writing and using a program computing solutions to ordinary differential equations. Hales and Ferguson solved the Kepler conjecture in 1998 [START_REF] Thomas | The Kepler conjecture[END_REF], establishing the densest way to pack spheres in three dimensions, originally needing 2000 hours of computations to solve huge nonlinear optimization problems. After four years, the 12 referees gave up and decided they could not possibly review the whole proof and computer code. In 2013, Helfgott proved the Ternary Goldbach Conjecture [START_REF] Helfgott | The ternary Goldbach conjecture is true[END_REF]. It states that every odd number n greater than 5 is the sum of three primes. The proof relies on two different natures of computations: the "manual" verification that the conjecture holds for n ≤ 10 27 , and the approximation of some numerical integrals. In order to compute the latter, he posted a question on Math Overflow [START_REF] Helfgott | Rigorous numerical integration[END_REF] to ask for a tool doing rigorous integration. Such proofs relying on heavy computations are controversial: programs tend to have bugs, and the number of bugs is reputed to be proportional to the number of lines of code. One single bug could render the whole proof invalid. For this reason, efforts have been made to dramatically increase confidence in these proofs. One possible route has been to make a new, "computer-free" proof of the same result, as a later paper [START_REF] Hutchings | Proof of the double bubble conjecture[END_REF] on the Double Bubble conjecture boasts. Another way has been to make computations as rigorous as possible, like Tucker who used interval arithmetic systematically in his proof. When the last two solutions are not suitable, there remains an ultimate step: building a complete mechanized formal proof of the result from elementary mathematical axioms to remove all possible doubt. This is the way chosen by Gonthier [START_REF] Gonthier | The four colour theorem: Engineering of a formal proof[END_REF] with the help of Werner when they proved the Four Color Theorem in the Coq proof assistant in 2005 and by Hales [START_REF] Hales | A formal proof of the Kepler conjecture[END_REF] who formalized with his team his proof of the Kepler conjecture in the Isabelle and HOL Light proof assistants.

These proof assistants are pieces of software which allow to build proofs from first principles, i.e. logical axioms and inference rules. Despite Bourbaki's opinion that "formalized mathematics cannot in practice be written down in full", the increase in computational power has in fact made this practical and numerous theorem provers have flourished in the past three to four decades.

It is possible to adopt a skeptical approach to certifying computations in a proof assistant. Suppose for example that we want to prove that some large number n is not prime. We could have a program produce two numbers 1 < p < n and q such that n = p • q.

(1.1)

Once we have the result, we can forget about the program and simply build a proof of Equation (1.1) inside a proof assistant. Here, the representation of natural numbers, the multiplication operation on them, and equality are underlying objects of the logic of this proof assistant.

The internal logic of the proof assistant can however also be made efficient enough to directly deal with non-trivial computations. In this case, the program doing the computation is written in a language internal to the logic and its result is natively accepted as valid. This is called an autarkic approach and it is the path chosen by Gonthier and Werner.

Whether using an autarkic or a skeptical approach to certifying computations, there still remains to formalize the remaining parts of a proof, which do not necessarily rely on computations. Furthermore, formalizations of mathematical proofs have not been limited to controversial, computation-heavy results; they have also been used for more traditional (but highly non-trivial) results, such as the Odd Order Theorem [START_REF] Gonthier | A Machine-Checked Proof of the Odd Order Theorem[END_REF]. Though this was not a contested result, one otherwise strong motivation for doing formal proofs in general is to reduce the risk of human mistakes while writing proofs. This aspect convinced Fields medalist Voevodsky to turn to formal proofs after a false result which he published in 1991 [START_REF] Mikhail | ∞-groupoids and homotopy types[END_REF] remained unquestioned until 1998, when Carlos Simpson showed that the result could not be true [114] [122]. Voevodsky finally found the mistake in his paper in 2013 [START_REF] Voevodsky | Univalent foundations[END_REF].

Whether large proofs are computation-heavy or not, mathematical results from all the theories involved in the proof must be systematically formalized in order to mechanize them in proof assistants. Much like when programmers write complex software, this involves careful choices of data structures and statements. Tools for helping the users of proof assistants in their tasks must be developed. Empirical experience shows that what may seem like a "good enough" design choice at a small scale can become unbearable technical debt for larger projects.

In fact, whereas one might think that formalizing results consists in copying paper proofs in a slightly more rigorous way, this activity leads to fresh insights on sometimes well-established mathematical theories. The formalization of the Odd Order theorem led Gonthier to come up with an original way to model abstract linear algebra using only matrices as a data structure to represent spaces, families of vectors [START_REF] Gonthier | Point-free, set-free concrete linear algebra[END_REF]. In his paper about the formalization of the Prime Number Theorem, Avigad [START_REF] Avigad | A formally verified proof of the prime number theorem[END_REF] gives several examples of "dumbing down" proofs, that is of coming up with more elementary proofs of a result in the absence of a library providing the necessary theory to directly adapt a textbook proof.

Building tools to help users develop formal proofs is not only about building libraries of theorems, it is also about automating away part of the tedious work of proving every little detail in a proof. Having a routine which can automatically prove that two rational fractions are equal, or establish simple inequalities between linear expressions can greatly improve the experience of the user. This is another way of "extending the number of operations which can be performed without thinking about them".

In this thesis, we illustrate three different ways to use computers to build mathematical proofs. We describe how experiments with existing matrix multiplication algorithms led to the discovery of a new family of parallel matrix multiplication algorithms. We then turn to formal proofs. We fully formalize a theorem by Apéry, in big part by certifying the outputs of a Computer Algebra System. Finally, we build an extension of the CoqInterval library which can compute fully verified approximations of integrals in a certified way.

Contributions

In this section, we list the contributions made in this thesis, by order of appearance and grouped following the structure of this document. Whenever our work has been published or refers to code stored online, we include the appropriate references.

Experiment

The contents of this chapter are partly the result of unpublished work done during the author's Master's degree in an internship supervised by Éric Schost. Part of this work was presented at the Journées Nationales du Calcul Formel1 in 2015 in Cluny and continued during a three-week visit to University of Waterloo in November 2016. The theorems of Section 4.3 and the Maple code were done as a collaboration with Éric Schost. The Ocaml code is mostly our own work, though it benefited along the way from insights and contributions from Éric Schost. The whole Specfun team is to be thanked for helping us come up with a hopefully clearer presentation of these results.

A New Class of Parallel Matrix Multiplication Algorithms We discover a new class of parallel matrix multiplication algorithms (Theorem 3). These algorithms can be obtained by applying a practical, finite version of Schönhage's τ -theorem to existing "approximation" algorithms by Pan and Schönhage and observing that one can remove the additional ε variable which had been introduced in these algorithms. We also provide a constructive proof of the τ -theorem which is inspired from this result.

A Software Tool for Manipulating Matrix Multiplication Algorithms

We provide a software tool for the symbolic manipulation of tensors encoding matrix multiplication algorithms. Our software can parse such tensors, compose them, apply the τ -theorem following our contribution in the previous paragraph, and export them to various formats including Maple, L A T E X and C++. This software requires git, and a recent version of Ocaml which can be installed using opam. It can be obtained and installed on Linux or Mac2 by typing the following three commands: git clone -b public https://gitlab.com/matrix-product-experiments/mpe cd mpe/caml sh build.sh

The specific options for using our tool are described in Section 4.4.2.

Certifying Computations

A Full Proof of Apéry's Theorem in Coq As a result of a collaboration with Frédéric Chyzak, Assia Mahboubi and Enrico Tassi, we provide a complete formal proof of Apéry's theorem stating the irrationality of the constant ζ(3) (Theorem 5). Part of this work has been published at ITP 2014 [START_REF] Chyzak | A Computer-Algebra-Based Formal Proof of the Irrationality of ζ(3)[END_REF] and a longer version, which also includes a simplification and formalization of a proof by Hanson regarding the asymptotic behavior of the sequence lcm(1,• • • , n), is currently in the works. Our personal contribution is the development of a library of arithmetic and number theory results concerning the sequences of Apéry's proof, the simplification and formalization of Hanson's proof, and assisting Frédéric Chyzak and Assia Mahboubi in coming up with a paradigm to state sound creative telescoping statements, to develop tactics both to handle big integers and to automatically resolve some linear inequalities over Z. Our Coq development can be found at https://specfun.inria.fr/ ~tsibutpi/ code/apery.zip.

Questioning the status of outputs of creative telescoping algorithms as straightforward proofs The main part of our proof of Theorem 5 is a validation of the outputs of creative telescoping algorithms which were intended to be used as a certificate. We question their status as certificates providing complete proofs of identities (Section 8.1).

Certified Computations

This work was done in close collaboration with Assia Mahboubi and Guillaume Melquiond.

An extension of CoqInterval with proper and improper integrals

We provide an extension of the CoqInterval library in order to approximate both improper and proper integrals using a mix of floating-point computations, interval arithmetic and rigorous polynomial approximations. All computations happen inside the logic of the Coq proof assistant. Using this tool, we highlight several irregularities in the published literature concerning integral computations. It has been published at ITP 2016 concerning proper integrals [START_REF] Mahboubi | Formally verified approximations of definite integrals[END_REF], and a longer article has been accepted and will appear in a special edition of the Journal of Automated Reasoning to be published in 2017 [START_REF] Mahboubi | Formally verified approximations of definite integrals[END_REF].

An extension of Coquelicot with a generalized Riemann integral In order to build the tool described in the previous paragraph, we extended the Coquelicot [START_REF] Lelay | Repenser la bibliothèque réelle de Coq: vers une formalisation de l'analyse classique mieux adaptée[END_REF] [START_REF] Boldo | Coquelicot: A user-friendly library of real analysis for Coq[END_REF] library so as to improve the handling of improper integrals (Section 12.2.5).

Notations and Conventions

In this section, we present some general notations. More context-specific notations will be introduced in each of the independent three parts.

In the context of both Ocaml and Coq code, as is usual, we will write function application in what is called a currified manner, that is f x instead of f (x). This generalizes to more parameters. For example, if f is a function with three parameters x, y and z, then we will write the application of f to x, y and z as f x y z. The partial application f x y is a function which maps z to f x y z.

If k ∈ N and if D ⊆ R, we will denote by C k (D) the set of functions which are k times continuously differentiable on D, and by C k the set C k (R).

We will denote by R[X] the set of polynomials with real coefficients, and if n ∈ N, R n [X] will be the set of polynomials P ∈ R[X] of degree at most n.

If X is a finite set, we denote by |X| the cardinal of X.

For n, k 1 ,• • • , k l ∈ N, with k 1 +• • • + k l = n, the quantity n k1,•••,k l := n! k1!•••k l !
is called a multinomial coefficient. If µ ∈ N l , we will write n µ for n µ1,•••,µ l . We also have

n k 1 ,• • • , k l = l i=1 k 1 +• • • + k i k i . (1.2)
In particular, n k1,•••,k l ∈ N. The multinomial theorem states that

n k1,•••,k l is the coefficient of x k1 1 • • • x k l l in the formal expansion of (x 1 +• • • + x l ) n :
Theorem 1 (Multinomial Theorem). Let A be a commutative ring and let x 1 ,• • • , x l ∈ A. Then for n ∈ N,

(x 1 +• • • + x l ) n = k1+•••+k l =n n k 1 ,• • • , k l x k1 1 • • • x k l l .
where we take the convention 0 0 = 1.

Programming Languages and Libraries

In this part, we present the various software tools used in the rest of the thesis. First, Ocaml and Coq are two pieces of software from Inria whose histories are closely bound together, and which were used extensively in the works presented in this thesis. The functional programming language Ocaml was instrumental in Part I in order to build our software tool to manipulate tensors representing matrix product algorithms and in Part III to build prototypes which were later implemented in Coq. The Coq theorem prover was used to prove Apéry's theorem in Part II and to build and prove correct our extension of CoqInterval to approximate integrals in Part III.

Ocaml

Ocaml [START_REF] Leroy | The OCaml system (release 3.12): Documentation and user's manual[END_REF] is a general purpose programming language of the family of ML languages. ML is the meta-language of Milner's proof assistant LCF (Logic for Computable Functions) [START_REF] Milner | The definition of standard ML: revised[END_REF] [START_REF] Gordon | Edinburgh LCF: a mechanized logic of computation[END_REF]. One core idea of LCF was to represent provable theorems as elements of an abstract type theorem which could only be instantiated with predefined axioms and elementary inference rules. This way, any element of the type theorem would be provable by construction, since it had to have been built from first principles.

On top of predefined base types such as nat for natural numbers, bool for booleans and float for floating-point numbers, Ocaml provides the possibility to define inductive types, which are recursive data types well-suited to encode elementary mathematical and computer science structures. For example, a tree could be defined as type tree = | Leaf | Node of tree * tree and an example of an object of type tree would be Node(Node(Leaf,Leaf),Leaf).

Ocaml has strong typing. This means two things. First, every object has a type. Secondly, in a valid Ocaml program, the only way that the object f x is a valid object of type B is if the function f can be given the type A -> B and if x has the type A. Although this may seem a trivial property, many languages do not have it. For example in Python, the math.exp function has type float -> float, but math.exp('foo') is accepted by the Python compiler in a function definition. However, if this line of code is ever executed it will result in a failure at runtime.

Typing in Ocaml is also static, meaning that programs are typechecked as they are compiled and not at runtime. If the compiler accepts an Ocaml program, then it must be well typed.

Experience shows that strong and static typing, which were initially motivated by the theorem proving usecase of ML, allow for the compiler to automatically detect many of the mistakes which are easy to make in more loosely typed or non-typed programming languages.

Another distinguishing aspect of Ocaml is that it is a functional programming language: functions are first-class objects, and they can be given as input or returned as outputs of programs. On top of this, Ocaml also has more common features such as imperative primitives (while and for loops) and an object-oriented paradigm.

Coq

Coq is a proof assistant. This means that it provides a work environment which one can use to write mathematical statements and mathematical proofs of these statements. It can also be used to write programs, program specifications (the description of the properties a program's outputs should have) and proofs that a program satisfies a given specification.

Coq is based on a formal language called the Calculus of Inductive Constructions. This language consists of objects called terms, which are assigned a type. Types are themselves terms, and as such they themselves have a type. This is unlike Ocaml where only base objects have types. Furthermore, types are dependent: they are allowed to depend on a term. For example, one can build the type of matrices m × n for given m, n ∈ N.

Elementary mathematical objects such as booleans, natural numbers or rationals can be represented in Coq using inductive types (similarly to Ocaml). Through the Curry-Howard isomorphism, there is a correspondence between logical statements in intuitionistic logic and types. If t is a type and A is the corresponding logical statement, then a term p of type t can be seen as a proof of A. Thus, in Coq, proving the statement A amounts to building a term p which typechecks to t: in other words, proving and programming are the same activity in this context.

The main activity of the Coq software consists in type-checking, which as explained above is how proofs are checked to be correct. Crucially, one of the typing rules in particular called the conversion rule allows for computations to be done inside Coq's logic. The part of Coq's code which handles typechecking is very sensitive: a bug could allow one to prove false statements. For this reason, this code is maintained in a relatively small trusted kernel which is kept separated from the rest of the code and to which any modification is closely examined.

As seen in Figure 1.4.2, a Coq session consists mainly of two windows: a text editor where one writes definitions, statements and proofs and a goal window where the remaining goals (i.e. statements which have not yet been proved) are listed. The part of the text highlighted in blue represents the lines of code already processed by Coq, while the remaining lines have not. When doing proofs, one typically writes and executes one line at a time. In order to make proofs less tedious, one can write and use procedures which partially automate the task ot theorem proving. This is done in a high-level language called Ltac which resembles Ocaml. Crucially, these tactics can only participate in building a proof, but are completely absent from checking the proof: thus they are allowed to do arbitrary modifications of the current state of the proof without hindering the trust in the kernel of Coq. Several examples of tactics, including of our own, will be described in later chapters.

More contextual information about Coq and Coq libraries will appear about the encoding of mathematical structures and the handling of computations (Part II) and about Real Analysis and tactics (Part III).

Finally, it should be noted that there exist many theorem provers, that are actively developed and used. Not all of them share the same logic as Coq; however many of them share a common ancestor called LCF. The choice of a theorem prover can have a significant impact on the way proofs are developed. The Isabelle theorem prover [START_REF] Paulson | Isabelle: A generic theorem prover[END_REF] may use classical rather than intuitionistic logic. The Lean Theorem Prover [START_REF] De Moura | The Lean Theorem Prover (System Description)[END_REF] is very similar to Coq, but comes with fewer mathematical libraries as it is more recent.

Part I

Experiment

Chapter 2

Introduction

When one tries to gain intuition on a mathematical problem or a conjecture, one naturally turns to a mix of mental and paper computations. In many cases, using or writing computer programs can significantly extend the range of the achievable. Computers can help visualize data, as in the simple case of plotting a function. They can also be used to check that a conjecture isn't trivially false. Pierre de Fermat famously conjectured that all numbers of the sequence F n := 2 2 n + 1, for n ∈ N, are prime numbers [START_REF] Euler | Observations on a certain theorem of Fermat and on others concerning prime numbers[END_REF]. He manually checked the first five instances 3, 5, 17, 257, and 65537.

Had he had access to a computer, Fermat hardly would have missed the chance to try out the very next element of the sequence, 2 2 5 + 1 = 4294967297. He could have written a basic primality checking function, which runs through all integers k less than the square root of its input n to see if n ≡ 0 [k]. The Ocaml code in Listing 2.1, once compiled, runs in less than 2 milliseconds on a standard laptop:

let is_prime n = let rec aux k = if k*k > n then true else if n mod k == 0 then (Printf.
printf "not prime: %d.\n" k; false) else aux (k+1) in aux 2;; This was in fact pointed out by Leonhard Euler [START_REF] Euler | Observations on a certain theorem of Fermat and on others concerning prime numbers[END_REF] without the help of a computer: "I have observed after thinking about this for many days that this number can be divided by 641".

let rec pow a = function | 0 -> 1 | 1 -> a | n -> let b = pow a (n / 2)
Attempting to disprove a conjecture with a program does not always work, if only because some conjectures may actually be true, or at least counterexamples may be harder to find. Alan Turing used a computer with only 25, 600 bits of memory to check that the first 1, 104 zeros of the Riemann zeta function were consistent with the Riemann hypothesis [START_REF] Hejhal | Alan Turing and the Riemann zeta function[END_REF]. The latest such effort of which we are aware, by Xavier Gourdon [START_REF] Gourdon | The 10 13 first zeros of the Riemann Zeta function, and zeros computation at very large height[END_REF], computes 10 13 zeros of the zeta function, all consistent with Riemann's conjecture.

Computers may also be used to come up with conjectures. One example of this is the technique of guessing employed to find linear recurrence relations with polynomial coefficients satisfied by a sequence [10] [112].

This small part has two chapters. First we focus on a combinatorics conjecture which states that a square of size 2 n-2 × 2 n-2 can be filled without overlap using tiles made up of pairs of paths of size n. We illustrate how, using the right software tool, it can be simple to toy with a nontrivial conjecture and we increase the best known value up to which this conjecture is true from n = 5 to n = 7.

Then, we consider the well-known problem of the complexity of matrix multiplication: how many arithmetic operations C(n) in a field K does it take to multiply two square matrices of size n with coefficients in K? For a long time, it was believed that the usual naive algorithm in O(n 3 ) was optimal. In 1969, Strassen [START_REF] Volker Strassen | Gaussian elimination is not optimal[END_REF] first proved that this is not the case by introducing an algorithm with complexity O(n log 2 7 ). He started a race to find lower bounds on the smallest constant ω such that C(n) = O(n ω ). The idea of introducing a variable ε in the algorithms by Bini [START_REF] Bini | O(n 2.7799 ) complexity for n*n approximate matrix multiplication[END_REF] led to the result ω < 2.79. Improvements quickly became of a purely theoretical nature: Schönhage's τ -theorem [START_REF] Schönhage | Partial and total matrix multiplication[END_REF], which allowed for the impressive improvement to O(n 2.55 ), is typically presented more as a result on the real number ω than as a concrete algorithm (with perhaps one exception [2]).

By contrast, in this work, we are interested in the explicit description of concrete algorithms. We obtain implementations of some of these more "theoretical" algorithms and make observations, for matrices of reasonable size. The main contributions in this chapter are of two kinds.

First, we provide a constructive proof of Schönhage's τ -theorem (Section 4.3) and introduce a new way to obtain "ε-free" algorithms to compute independent parallel matrix multiplications (Section 4.3.2), and to pre-compute their complexity without directly building them.

Secondly, we present an implementation in Ocaml of a new tool for the symbolic manipulation of tensors representing matrix product algorithms (Section 4.4). This tool allows us to actually build the tensors mentioned in the previous paragraph, and even to export them to C++ or as L A T E X pseudo-code. Our hope is that this tool can in turn be used for experiments.

Chapter 3

Experimenting with a Tiling Conjecture

In this chapter, we experiment with a conjecture submitted to us by Alin Bostan, using the Combine [START_REF] El | Combine: an ocaml library for combinatorics[END_REF] library developed in Ocaml by Rémy El Sibaïe and Jean-Christophe Filliâtre. First we state the problem (Section 3.1). Then we describe the theory underlying the Combine library's algorithm for solving tiling problems (Section 3.2). We show the results we found (Section 3.3). Finally, we conclude (Section 3.4).

Problem Description

In the following, a step is a vector of Z 2 . The only steps we consider are East (1, 0) and North (0, 1). A path p is a finite sequence of steps; p is said to go from A ∈ Z 2 to B ∈ Z 2 if A + s∈p s = B. If p 1 and p 2 are two paths which start and end at the same points A and B of Z 2 , we call tile the finite area enclosed by the two paths.

The problem we consider is stated by Woan, Shapiro and Rogers [START_REF] Woan | The catalan numbers, the lebesgue integral, and 4n-2[END_REF] in the following way. Definition 1. A path-pair of length n ≥ 2 is a pair of paths of Z 2 which both start at the origin (0, 0), consist of n steps and meet again for the first time after exactly n steps. All steps in these paths go either East or North. We denote by P n the set of such path-pairs.

An example of such a tile is given in Figure 3.1: The paths are respectively NNNEEEEE and EEENNEEN where the letter N stands for North and E stands for East.

The authors go on to state and prove the following two relevant lemmas:

Lemma 1 (cardinality of P n ). The cardinality of P n is C n-1 , where (C n ) n∈N is the sequence of Catalan numbers defined by C n = 1 n+1 2n
n . Lemma 2 (total area of P n ). The sum of the areas of the path-pairs of P n is 4 n-2 .

The conjecture we are interested in is the following: Lemma 2 justifies that the tiling is not completely implausible, but aside from the fact that the areas of the tiles and the square match we are not aware of any intuition as to why Conjecture 1 may be true. However, if we draw the tiles for n = 4 or n = 5 (Figure 3.2), we see that Conjecture 1 holds for these values. The authorship of Conjecture 1 is attributed to a student called Schwärzler [START_REF] Fürlinger | q-catalan numbers[END_REF]. We could not find references to solutions for n > 5.

Exact Matrix Cover, The Dancing Links Algorithm and Combine

In this section, we describe some of the theory underlying the Combine library and the algorithm it uses to look for tilings such as the ones displayed in Figure 3.2. First, we describe the general Exact Cover problem:

Definition 2 (Exact Cover Problem). Let X be a set and let S be a collection of subsets of X. An exact cover of X is a set S * ⊆ S such that every element of X appears in exactly one element of S * . In other words, S * is a partition of X.

When X is finite, finding an exact cover can be reformulated in the following way. Build a matrix M with |X| columns and |S| rows, indexing the columns (respectively the rows) with the elements x i of X (respectively the elements S j of S), so that m i,j is 1 if x i ∈ S j and 0 otherwise. Then a set S * ⊆ S is a solution to the exact cover problem if and only if the sub-matrix obtained by keeping only the rows indexed by the elements of S * has exactly one 1 per column.

This reformulation is called the exact matrix cover (EMC) problem. It is NPcomplete [START_REF] Michael | Computers and intractability[END_REF]. In 2000, Knuth published the elegant "Dancing Links" algorithm for solving this problem [START_REF] Knuth | Dancing links[END_REF]. The algorithm's efficiency relies on the clever use of a doubly-linked-list based data-structure to represent matrices in memory.

Many well-known problems can be reformulated as instances of the EMC problem: for instance the n-queens problem (how does one place n queens on a n × n chessboard so that they don't attack one another?) or solving a Sudoku puzzle. The tiling problem of Conjecture 1 is in fact also an instance of the EMC problem.

The reformulation of the problem in Conjecture 1 as an EMC problem consists in building a matrix which describes all possible positions of the tiles on the grid. The grid G is formed of 4 n-2 squares g. Let T be the set of tiles (in our case, P n ) and let T be the set of shapes which can be obtained from T : the set T is obtained by applying all permitted rotations on T , so that each element of T is identical to (one element of T and) at most three other elements of T , up to a rotation. The "at most" here is because duplicates due to the symmetries of a tile are removed for efficiency. To each element t of T , we associate a set P t of subsets of G. Each element of P t is a legal position of t on G.

We form a matrix M with |G| + |T | columns, the 4 n-2 first being indexed by the elements of G and the last ones by the elements of T . M has one line for each position p on G of each shape t in T . Hence M has t ∈T |P t | lines. Thus, each of the first set of columns has a 1 on every line indexed by a position which covers the corresponding square. Each of the second set of columns, indexed by a tile t, has a 1 on every line indexed by a shape obtained by a rotation of t. Each line indexed by a shape t and a position p has a 1 on every square that is covered by putting t at position p, and a 1 on the tile t from which it was obtained as a rotation.

Suppose we find a solution to this EMC problem. Since this solution consists in a set of lines with exactly one 1 on each column, it can be seen to be a tiling of G with the elements of T up to rotation.

The Combine library provides a powerful solver for tiling problems using the Exact Matrix Cover encoding and the Dancing Links algorithm [START_REF] Knuth | Dancing links[END_REF].

Solutions for n = 6 and n = 7

We used Combine to confirm the conjecture for n ≤ 7, by first generating the tiles in Ocaml for each n.

Conclusion

We have improved the best known value up to which Conjecture 1 holds from n = 5 to n = 7, using the pre-existing Combine library to find tilings of a square. Unfortunately, our computer does not have enough memory for the problem to be stored in memory on our laptop for n = 8. A parallel approach could be justified and will be attempted later, but it would need some modifications to the Combine library.

Our point was to illustrate that often, experimenting is not only about having a computer at hand and typing a simple program such as in Listing 2.1. Having software tools tailored for specific computations can go a long way towards effectively tackling a problem.

We also tried a different (not documented here) approach for finding tilings using a linear programming encoding on a suggestion of Marc Mezzarobba, but unfortunately this turned out to yield worse times by a factor 10.

Chapter 4

Experiments with Matrix Multiplication Algorithms

It may seem like a paradox that a topic so typical of Computer Science as the complexity of multiplying matrices should serve as an illustration of how experimenting with a computer can help gain insight on results. In fact, there is a widely held belief that past some point, "fast" matrix product algorithms are a purely theoretical amusement which has no bearing on matrices of sizes typically encountered in the real world. We do not categorically dispute this notion, but the work we present here is an attempt to reconcile these theoretical results with more practical approaches to the design of algorithms.

We start by defining the problem of the complexity of Matrix Multiplication and related concepts (Section 4.1). We recall the major milestones in improving the ω constant which underlies it (Section 4.2). We present a constructive proof of Schönhage's τ -theorem; we also present a result which allows to build more practical matrix product algorithms from seemingly impractical ones (Section 4.3). This result was preceded and followed by custom software implementations for manipulating matrix product algorithms symbolically; we describe these implementations and the observations we were able to make thanks to them (Section 4.4).

Definitions

Matrix Multiplication

We fix a field K. We call M m,n (K) the space of matrices of size m × n over the field K. The matrix multiplication mapping

Let A = (a i,j ) 1≤i≤m 1≤j≤n ∈ M m,n (K) and B = (b i,j ) 1≤i≤n 1≤j≤p ∈ M n,p (
m, n, p : M m,n (K) × M n,p (K) -→ M m,p (K) (A, B) → A • B is a bilinear map.
We define the complexity of a matrix multiplication algorithm A as the number of arithmetic operations (additions and multiplications) in K computed by A. Whenever we focus on the asymptotic complexity of matrix multiplication, we will tend to only count multiplications as they asymptotically dominate additions in all the constructs we use.

It is a conjecture that the complexity of matrix multiplication is independent of K. Until it is proved, all results are implicitly parameterized by K. When there is an ambiguity, we will be specific about the ambient field K.

Rank of a bilinear map

If U is a vector space over K, we denote by U * the space of linear forms over U . Definition 3 (tensor product of two vector spaces). Let U , V be two vector spaces over K. There exists a unique (up to isomorphism) vector space over K, U ⊗ V called the tensor product of the spaces U and V , and a bilinear map

Φ : U × V → U ⊗ V
such that for any vector space W over K and bilinear map f : U × V → W , there exists a unique linear map f :

U ⊗ V → W such that f = f • Φ. The map Φ is denoted by Φ(x, y) = x ⊗ y for x ∈ U , y ∈ V . Furthermore, ⊗ is associative, but not commutative.
Let U , V and W be finite dimensional vector spaces over K and let f : U × V -→ W be a bilinear map. Definition 4. The tensor rank of f is the smallest integer r such that there exist x i ∈ U * , y i ∈ V * and z i ∈ W satisfying:

∀ u ∈ U, v ∈ V, f (u, v) = r i=1 x i (u)y i (v)z i . (4.1)
We denote the tensor rank r by R(f ).

For any x i , y i and z i satisfying Equation (4.1), the element

t f = r i=1 x i ⊗ y i ⊗ z i of the space U * ⊗ V * ⊗ W is called a tensorial representation of f .
Let us fix bases for U , V and W . Up to the canonical isomorphism between U and U * (respectively V and V * , W and W * ), t f can be seen as an element of U ⊗ V ⊗ W , a product which gives a similar role to U , V and W . Thus, more generally, we define the tensor rank R(t) of t ∈ U ⊗ V ⊗ W as the minimum r such that t can be written as a sum of r elementary tensors x i ⊗ y i ⊗ z i , i.e. that t can be written

t = r i=1 x i ⊗ y i ⊗ z i . with x i ∈ U , y i ∈ V and z i ∈ W for 1 ≤ i ≤ r.
It is easy to see that for any t f as in Definition 4, R(t f ) = R(f ) by definition, so that we can now indifferently speak of the tensor rank of f or the tensor rank of t f .

Note that the tensor rank of f or t depends on the field K and should be written R K (f ) or R K (t); when there is no ambiguity we will keep the lighter notation R(f ) and R(t).

Example: Identifying row and column matrices with K n , the mapping

1, 2, 1 :K 2 × K 2 → K ((u, v), (x, y)) → u • x + v • y has tensor rank ≤ 2. Indeed, if we take π 1 : K 2 → K to be the first projec- tion (u, v) → u and π 2 : K 2 → K to be the second projection (u, v) → v, then if A, B ∈ K 2 1, 2, 1 (A, B) = π 1 (A)π 1 (B) 1 + π 2 (A)π 2 (B) 1 .
In fact, as soon as

|K| ≥ 2, R K ( 1, 2, 1 ) = 2.

The particular case of the matrix multiplication

Computing the tensor rank of m, n, p gives a hint of the complexity of matrix multiplication. In fact, it can be seen that the tensor rank of n, n, n is less than twice the arithmetic complexity of multiplying two n × n matrices in K. This is proved for example in a book by Burgisser et al. [START_REF] Bürgisser | Algebraic complexity theory[END_REF].

Each algorithm computing the matrix multiplication m, n, p that we consider corresponds to a particular way of writing m, n, p as a sum of elementary tensors. For instance, the naive matrix product algorithm can be written

i,j,k a i,k ⊗ b k,j ⊗ c i,j = m, n, p (4.2)
so that in particular R( m, n, p ) ≤ mnp.

(4.
3)

The substitution of tensor rank for arithmetic complexity leads us to define the following constant:

ω := inf{τ ∈ R | R( n, n, n ) = O(n τ )}.
In this part, we are thus interested in minimizing the rank R( n, n, n ). Equation (4.3) yields a first upper bound for ω, as the sum has mnp terms, i.e

n 3 if m = n = p : ω ≤ 3. (4.4)
Obtaining a lower bound requires introducing more technical material outside of the scope of the present text, but the reader can refer to the book [START_REF] Bürgisser | Algebraic complexity theory[END_REF] by Burgisser et al. to see how to obtain:

ω ≥ 2. (4.5)
Let us consider some aspects of the tensor notation as regards matrix multiplication algorithms.

First, consider the tensor product ⊗ of two tensors: the object m 1 , n 1 , p 1 ⊗ m 2 , n 2 , p 2 can be seen as the tensor obtained when associating to each of the elementary tensors of m 1 , n 1 , p 1 a copy of the tensor m 2 , n 2 , p 2 . This is exactly what we do when we perform a block matrix product: for a moment, we see each block as a coefficient, only to remember later that it is actually a matrix.

It should hence be no surprise that for all m 1 , n 1 , p 1 , m 2 , n 2 , p 2 we have:

m 1 , n 1 , p 1 ⊗ m 2 , n 2 , p 2 ∼ = m 1 m 2 , n 1 n 2 , p 1 p 2 , (4.6)
where the meaning of ∼ = is the following: let n denote {1,• • • , n}. There exists a bijective function

f : m 1 × n 1 × p 1 × m 2 × n 2 × p 2 -→ m 1 m 2 × n 1 n 2 × p 1 p 2 ((i 1 , j 1 , k 1 ), (i 2 , j 2 , k 2 )) → (m 2 (i 1 -1) + i 2 , n 2 (j 1 -1) + j 2 , p 2 (k 1 -1) + k 2 ) such that if we replace a i,k ⊗ b k,j ⊗ c i,j ⊗ a i1,k1 ⊗ b k1,j1 ⊗ c i1,j1 with a i2,k2 ⊗ b k2,j2 ⊗ c i2,j2 , where (i 2 , j 2 , k 2 = f ((i, j, k), (i 1 , j 1 , k 1
)), then the tensor product on the left becomes the tensor on the right. Another important operation for our purposes is taking a direct sum of tensors. Let

t 1 ∈ U 1 ⊗ V 1 ⊗ W 1 and t 2 ∈ U 2 ⊗ V 2 ⊗ W 2 . Then the direct sum t 1 ⊕ t 2 , which lives in the space (U 1 ⊗ V 1 ⊗ W 1 ) ⊕ (U 2 ⊗ V 2 ⊗ W 2 )
can be interpreted as the independent parallel computation of t 1 and t 2 . We will denote the sum t ⊕ t• • • ⊕ t s times by s t. In practice, as long as t 1 and t 2 do not share any variables, we will simply obtain t 1 ⊕ t 2 as t 1 + t 2 .

Algorithmic interpretation: In this context, any elementary tensor can be seen as a description of a program computing a matrix product, where the first two components of the tensor provide the coefficients to be multiplied, and the third one, a weighted list of the result matrix; i.e. the one where the results will eventually be stored.

Hence the elementary tensor (a

1,2 + a 3,5 ) ⊗ b 2,4 ⊗ (c 1,4 + 2 • c 2,4 ) reads as the algorithm tmp ← (a 1,2 + a 3,5 ) • b 2,4 c 1,4 ← tmp c 2,4 ← 2 • tmp
We consider a few important properties of the tensor rank, which also justify why it is called this way.

Usual properties of the rank

The following properties of the rank are stated without their proofs, which can be found for example in Burgisser et al. [START_REF] Bürgisser | Algebraic complexity theory[END_REF]. We will give an intuition of what these properties mean in the context of tensors representing matrix multiplication algorithms. Definition 5. Let π be a permutation of {1, 2, 3} and

t ∈ U 1 ⊗ U 2 ⊗ U 3 . We define π(t) ∈ U π -1 (1) ⊗ U π -1 (2) ⊗ U π -1
(3) to be the corresponding permutation of components on elementary tensors.

We have

π ( m 1 , m 2 , m 3 ) = m π -1 (1) , m π -1 (2) , m π -1 (3) .
Lemma 3. Let f : U × V -→ W and f : U × V -→ W be bilinear maps, and t ∈ U ⊗ V ⊗ W , then we have:

• R(f ⊕ f ) ≤ R(f ) + R(f ); • R(f ⊗ f ) ≤ R(f ) R(f ). • R(π(t)) = R(t)
The first item expresses that it is not more complicated to compute the sum of two independent tensors than to compute one and then the other. The second item simply states that the number of terms in the expansion of the product of two sums is the product of the number of terms in each sum. In terms of matrix multiplication, it implies that if we can compute m 1 , n 1 , p 1 in r 1 multiplications and m 2 , n 2 , p 2 in r 2 multiplications, then we can compute m 1 m 2 , n 1 n 2 , p 1 p 2 in r 1 r 2 multiplications: divide matrices into blocks so that you can do r 1 block products, and do every product between two blocks using r 2 multiplications.

Finally, we prove the third point in the particular case of interest for us, that of matrix multiplication. Consider the trace trilinear form (A, B, C) → T r(ABC) on matrices of sizes (m, n), (n, p) and (p, m) for some positive integers m, n and p. Notice that for 1 ≤ i ≤ m and 1 ≤ j ≤ p, the coefficient of c j,i in T r(ABC) (seen as a polynomial expression) is the matrix coefficient of index (i, j) of the matrix product AB (corresponding to m, n, p ). In the same way, the coefficient of b k,j , 1 ≤ k ≤ n, 1 ≤ j ≤ p in T r(CAB) is the matrix coefficient of index (j, k) of the product CA (corresponding to p, m, n ). A similar property holds for T r(BCA). Any permutation of m, n, p can be obtained, possibly by transposing A, B or C. By observing that T r(ABC) = T r(CAB) = T r(BCA), we can then convert any explicit tensor for m, n, p into a tensor for a permutation of m, n, p by a simple syntactic rewriting which does not change the number of terms, thus all six permutations have the same rank.

Lemma 4. R( n, n, n ) ≤ r ⇒ n ω ≤ r ⇒ ω ≤ log n (r)
This a consequence of the second item of Lemma 3. From Lemmas 3 and 4 we can deduce

Lemma 5. R( m, n, p ) ≤ r ⇒ (mnp) ω ≤ r 3 .
Proof. Let π be the cyclic permutation mapping 1 to 2, 2 to 3, and 3 to 1.

mnp, mnp, mnp ∼ = m, n, p ⊗ π ( m, n, p ) ⊗ π 2 ( m, n, p ) so that using the second and third items of Lemma 3, R( mnp, mnp, mnp ) ≤ r 3 which yields the result by Lemma 4.

Degeneration and border rank

The notion of border rank was crucial in improving the upper bound on ω. The general idea is to relax tensors to be, not just a particular way of writing a matrix product algorithm, but some perturbation of one containing a special variable ε. Such algorithms are usually called approximation algorithms.

In the following, U , V , W designate vector spaces over K. Definition 6. We define U K[ε] to be the extension of scalars of U to the ring

K[ε].
From now on, we will use the notation T for U ⊗ V ⊗ W (whose elements will be called plain tensors) and

T ε for U K[ε] ⊗ V K[ε] ⊗ W K[ε]
(whose elements will be called tensors with ε's or polynomial tensors).

The following definition may seem a bit dry, but the idea is actually simple: suppose the plain tensor we are trying to write as a minimal sum of elementary tensors is t 1 ∈ T . The previous problem we considered until now was to find a way to directly write t 1 as a sum of elementary tensors. Now, we relax this by finding some tensor t 2 ∈ T ε such that t 1 is the "tensor of lowest degree of t 2 ":

t 2 = ε q-1 t 1 + ε q t 3
where t 3 ∈ T ε . Definition 7. Let t 1 ∈ T and t 2 ∈ T ε . We say that t 2 is a degeneration of order q of t 1 (denoted t 1 q t 2 ) if there exists t 3 ∈ T ε such that

t 2 = ε q-1 t 1 + ε q t 3 .
The border rank R(t) of t ∈ T is the smallest r such that there exist

u i ∈ U K[ε] , v i ∈ V K[ε] , w i ∈ W K[ε]
and q ≥ 1 such that:

t q r i=1 u i ⊗ v i ⊗ w i .
When there is no ambiguity, we will leave q implicit and write t 1 t 2 .

We want to show that we can extract a tensor of low rank computing a tensor t using such a degeneration of t. First we state that the degree q of approximation does not grow too fast (linearly) when we take N -th powers: Lemma 6. Let t ∈ T , t 1 ∈ T ε and N ∈ N.

t q t 1 ⇒ t ⊗N N (q-1)+1 t ⊗N 1 .
This enables us to establish a result analogous to Lemma 4, but for tensors with ε's:

Lemma 7. Suppose R ( m, n, p ) ≤ r.
Then

ω ≤ 3 log(r) log(mnp) .
This means that if we can find a degeneration of a matrix product with r (polynomial) multiplications, it yields the same bound on ω as if it were a plain tensor. The proof of Lemma 7 relies on the following lemma:

Lemma 8. Let t ∈ T and r ∈ N. Then if R(t) ≤ r, R(t) ≤ r q(q + 1) 2
Proof. Suppose R(t) ≤ r, and write t q r i=1 u i ⊗ v i ⊗ w i as in Definition 7

(with possibly some terms equal to zero if R(t) < r). Expand each elementary polynomial tensor

u i (ε) ⊗ v i (ε) ⊗ w i (ε) into a polynomial in ε with coefficients in U ⊗ V ⊗ W .
Keep only the terms of degree q -1: there are at most as many as there are triples of nonnegative integers (i, j, k) whose sum is q -1 (because each monomial is built from the product of three monomials in ε), which is

q(q+1) 2
. The tensor obtained by removing terms of higher degree is ε q-1 t, and it can be written as the sum of at most r q(q+1) 2 elementary tensors, so that R(t) ≤ r q(q + 1) 2 .

We can now conclude the proof of Lemma 7: R ( m, n, p ) ≤ r so that R m, n, p ⊗N ≤ r N for N ∈ N. Thanks to Lemma 6, we get R m, n, p ⊗N ≤ r N • P (N, q) where P (N, q) is a bivariate polynomial in N and q. Now we can apply Lemma 7 and let N tend to infinity and we get our result.

The intuition behind this result is that finding the coefficient of ε q-1 is asymptotically smaller than the overhead cost of computing with polynomials, since the degree of the approximation grows linearly when the rank grows exponentially with N . This is in contrast with the case of low values of N , for which this overhead is less justifiable. Thus, although border rank constitutes a formidable mathematical jump forward and allowed for spectacular improvements on ω, it is often mentioned as one of the reasons for matrix product algorithms to be deemed impractical, because extracting the coefficients only becomes faster for extremely large matrices.

A Quick Overview of Some Matrix Product Algorithms

This chapter describes some of the landmark results which have improved the best lower bound on ω from the naive estimate ω < 3 (see Equation (4.4)) to the current record by Le Gall [START_REF] François | Powers of tensors and fast matrix multiplication[END_REF]: ω < 2.37287. It is not exhaustive; for this we suggest a few references [116][22][2] [START_REF] Williams | Multiplying matrices faster than Coppersmith-Winograd[END_REF]. Dumas and Pan recently published a review [START_REF] Dumas | Fast matrix multiplication and symbolic computation[END_REF] with an emphasis on practical matrix multiplication, which we can informally define as the multiplication of matrices of size n ≤ 10 6 .

Strassen's algorithm (1969)

Strassen's algorithm was the first to break the ω < 3 barrier, much like Karatsuba was the first to break the barrier of O(n 2 ) for multiplying polynomials.

Strassen [START_REF] Volker Strassen | Gaussian elimination is not optimal[END_REF] proves that R ( 2, 2, 2 ) ≤ 7

by exhibiting an algorithm to compute 2, 2, 2 in 7 multiplications. From there, invoking Lemma 4, we get ω ≤ log 2 (7) ≈ 2.81. Note that Winograd proposed an improvement on Strassen's algorithm which reduces the number of additions [START_REF] Winograd | On multiplication of 2× 2 matrices[END_REF], but also proved that 7 multiplications is optimal.

Bini's Approximate Algorithms (1979)

The idea to use the border rank, presented in Section 4.1.4 was first introduced by Bini [START_REF] Bini | O(n 2.7799 ) complexity for n*n approximate matrix multiplication[END_REF]. It was the first departure from the idea of finding "Strassenlike" exact matrix multiplication algorithms for a specific format with fewer multiplications than the naive algorithm. Bini introduced an approximative (i.e., using ε's) algorithm for the product

a 1,1 a 1,2 a 2,1 0 b 1,1 b 1,2 b 2,1 b 2,2 (4.7)
using 5 polynomial multiplications:

φ(ε) :=(a 1,2 + εa 1,1 ) ⊗ (b 1,2 + εb 2,2 ) ⊗ c 2,1 + (a 2,1 + εa 1,1 ) ⊗ b 1,1 ⊗ (c 1,1 + εc 1,2 )- a 1,2 ⊗ b 1,2 ⊗ (c 1,1 + c 2,1 + εc 2,2 ) -a 2,1 ⊗ (b 1,1 + b 1,2 + εb 2,1 ) ⊗ c 1,1 + (a 1,2 + a 2,1 ) ⊗ (b 1,2 + εb 2,1 ) ⊗ (c 1,1 + εc 2,2 )
From this, one can easily deduce an approximative algorithm for 3, 2, 2 by gluing two copies of φ (and changing indices appropriately). By using permutations, one can deduce approximative algorithms for 2, 3, 2 and 2, 2, 3 , and subsequently for 12, 12, 12 in 10 • 10 • 10 = 1000 multiplications.

Finally, this gives R ( 12 n , 12 n , 12 n ) ≤ 1000 n so that ω < 2.78 by Lemma 7.

Schönhage's τ -theorem (1981)

In 1981, Schönhage [START_REF] Schönhage | Partial and total matrix multiplication[END_REF] departed even further from Strassen's algorithm. While Bini had relaxed the goal of finding matrix product algorithms to that of finding approximate matrix product algorithms, Schönhage relaxed it even further to approximate, independent and parallel matrix product algorithms. The following theorem states that if we have a tensor t which is a degeneration of the sum of s independent matrix products of various sizes, we can deduce from it an asymptotic upper bound for the complexity of matrix multiplication of square matrices. This complexity depends on the sizes of these matrices and on the border rank of t.

Theorem 2 (Schönhage's τ -theorem). Suppose that for s ∈ N * and for positive integers

m i , n i , p i , 1 ≤ i ≤ s, we have R s i=1 m i , n i , p i ≤ r. Then, if β > 0 is such that s i=1 (m i n i p i ) β = r, we have ω ≤ 3 β.
An equivalent statement is that

s i=1 (m i n i p i ) ω ≤ r.
Schönhage then presents a tensor computing 4, 1, 4 ⊕ 1, 9, 1 with border rank 17, which immediately gives the corollary that ω < 2.55. This tensor is presented as an illustration of Pan tables in the next section.

Although this theorem has been a substantial tool to improve bounds on ω, the usual proof is non-constructive. By this, we mean two things. First, it relies on an implicit argument of existence of a family of algorithms whose complexity tends to the resulting value of ω rather than exhibiting one. Secondly, it does not give a way to build the algorithm concretely. This discourages one from building them and looking at the complexity of these algorithms, and encourages the idea that the resulting complexity is only interesting asymptotically.

Pan's Trilinear aggregation (1984)

In this section, we present an idea by Pan which consists in combining disjoint matrix product tensors to take advantage of Theorem 2. The tensors are presented as rows of three-column tables, and Pan's idea consists in adding them column by column. The tensor obtained this way then needs a correction term to remove terms of low degree: the hope is that the rank of this correction term is low. This idea is called trilinear aggregation and was extensively developed in a book [100] and is more clearly and concisely recalled in [START_REF] Dumas | Fast matrix multiplication and symbolic computation[END_REF]. Here is a simple example of such a table. This tensor is actually from Schönhage, but we adopt a "Pan table" view of it.

x i,0 y 0,k ε 2 z k,i εu 0,k,i εv k,i,0 w 0,0 (4.8)
Each line almost represents a different (and disjoint) matrix multiplication, and each column represents a component of the tensor. The almost is because the variables u 0,k,0 and v 0,i,0 are adjustment variables and do not correspond to actual matrix entries.

One should specify bounds, which are here i = 0..m -1 and k = 0..p -1. Thus the tensor represented by the table in Equation (4.8) reads:

t 1 := m-1 i=0 p-1 k=0 (x i,0 + εu 0,k,i ) ⊗ (y 0,k + εv k,i,0 ) ⊗ (ε 2 z k,i + w 0,0 ). (4.9)
One should also specify the value of the adjustment variables, under the form of linear constraints, like for instance

m-1 i=0 u 0,k,i = 0 (4.10) p-1 k=0 v k,i,0 = 0 (4.11)
The tensor t we are interested in obtaining is the coefficient of degree 2 in ε of the tensor in Equation (4.9). In order to have t 2 t 1 (see Definition 7), we need to get rid of terms of degree lower than 2. Thus we add the correction term This would be done in a naive way in (m

σ = - m-1 i=0 x i,0 ⊗ p-1 k=0 y 0,k ⊗ w 0,
• 1 • p) + (1 • (m -1)(p -1) • 1) = 2 m p-(m+p-1)

A Constructive Proof of the τ -Theorem and a New Family of Algorithms

In this section, we first present a constructive proof of the τ -theorem (Theorem 2). This proof then inspires a theorem about a new family of matrix multiplication algorithms.

Proof

Part of the following proof is inspired from that of [2,Chapter 7]. We assume that R (

s i=1 m i , n i , p i ) ≤ r and we take β ∈ R, β > 0, such that s i=1 (m i n i p i ) β = r.
Let us show that ω ≤ 3 β. We can assume that r > s, unless all products are 1, 1, 1 , which is not an interesting case: it would mean we are doing s independent scalar multiplications. We will make the assumption that we are not in that case.

Recall that if t is a tensor and s ∈ N, s t stands for

t ⊕ t• • • ⊕ t s times .
We should first remark that ⊗ and ⊕ behave just like multiplication and addition up to an isomorphism. Analogously to Theorem 1, we have for k ∈ N: In fact, Lemma 9 is an easier case of the τ -theorem when all disjoint matrix products are of the same size.

s i=1 m i , n i , p i ⊗k ∼ = µ µ1+•••+µs=k k µ s i=1 m µi i , s i=1 n µi i , s i=1 p µi i ( 4
Proof. The idea is to recursively compute m k , n k , p k . To this end, we see

m k , n k , p k as m k-1 , n k-1 , p k-1 ⊗ m, n, p .
We define the sequence (K k ) k≥1 by

K k = r s k-1 mnp + r k-1 i=1 s r i .
Notice that if k ≥ 1, we have

K k+1 = r s K k + r.
We now prove by induction on k ≥ 1 the assertion

H k := R( m k , n k , p k ) ≤ K k . • H 1 is true, since K 1 = mnp. • Suppose H k-1 is true for some k ≥ 2. Recall that m k , n k , p k ∼ = m k-1 , n k-1 , p k-1
⊗ m, n, p : this consists in cutting the matrix multiplication m k , n k , p k as a higher-level matrix multiplication of size m k-1 , n k-1 , p k-1 , but where the base operation is the lower-level matrix multiplication m, n, p . The higher level multiplication can be done by induction hypothesis in less than K k-1 operations, each being a m, n, p multiplication. Now, we cut this set of multiplications into K k-1 s stacks of s multiplications, padding the last stack with zeroes if necessary. We get the formula:

R( m k , n k , p k ) ≤ K k-1 s (R (s m, n, p ))) ≤ K k-1 s + 1 r = K k .
Our induction is proved. Remember that s r < 1 , so that the sum in K k can be bounded by

1 1-s r . Thus R( m k , n k , p k ) ≤ r s k-1 mnp + r 1 -s r .
Thus by Lemma 7,

kω 3 ≤ (k -1) log( r s ) log(mnp) + log mnp + r 1-s r log(mnp) .
Dividing on both sides by k and letting k tend to infinity, we get the expected result.

We can now use Lemma 9 to write from (4.16):

ω ≤ 3 log( r k S k ) log(M k N k P k )
. (4.17)

Consider the equality:

r k = s i=1 (m i n i p i ) β k = µ µ1+•••+µs=k k µ s i=1 (m i n i p i ) µi β .
(4.18)

The biggest term on the right is bigger than the average over the k+s-1 s-1 ≤ (k + 1) s-1 terms of the sum, so that we may choose the term S k M k , N k , P k in Equation (4.16) such that:

S k (M k N k P k ) β ≥ r k (k + 1) s-1 . (4.19)
However, this is not effective if we want to actually find S k . In practice, we may choose a central multinomial coefficient as in the following lemma.

Lemma 10

(Maximal Multinomial Coefficient). Let k, s ∈ N * . Write the euclidean division of k by s, k = qs + r. Then k! q! s-r (q+1)! r is maximal among all multinomial coefficients k ν1,•••,νs .
Proof. First, we prove that in a maximal tuple ν 1 ,• • • , ν s , all ν i necessarily satisfy q ≤ ν i ≤ q + 1. Notice first that for every element ν i < q there exists some ν j ≥ q + 1. Assume for convenience that i = 1 and j = 2:

k! ν 1 ! ν 2 ! • • • ν s ! < ν 2 ν 1 + 1 k! ν 1 ! ν 2 ! • • • ν s ! = k! (ν 1 + 1)! (ν 2 -1)! • • • ν s ! .
This contradicts our maximality assumption, which proves our first point. Assume there are t values of ν i which are equal to q, and s -t which are equal to q + 1. Then we have k = t q + (s -t) (q + 1) = k so that t = s -r, which proves our lemma.

In order to compute the asymptotics, we keep to Equation (4.19). Using Equation (4.17) and Equation (4.19) we get: and so by letting k tend to infinity, we have the desired result.

ω ≤ 3 log (k + 1) s-1 (M k N k P k ) β log(M k N k P k ) ≤ 3β + 3(s -1) log (k + 1) log(M k N k P k ) . Since • log(M k N k P k ) ≥ 1 β log(

Interpretation of this result

Using the τ -theorem, spectacular theoretical improvements can be made. Equation (4.14) is mentioned as being obvious in most papers and books on the topic, but it is not so straightforward to actually write this bijection. We will not write it down here, but we programmed it in Maple and later Ocaml. We wrote a procedure Extract(t, k, µ) which takes as input a tensor t such that

s i=1 m i , n i , p i t, an integer k ∈ N and a vector µ = (µ 1 ,• • • , µ s ) ∈ N s with µ 1 +• • • + µ s = k,
and returns a tensor (i.e., an algorithm) t such that

k µ s i=1 m µi i , s i=1 n µi i , s i=1 p µi i t .
Of course, if one uses a direct sum of naive matrix multiplications, one only gets k µ naive products of the type M k , N k , P k . Only tensors already constituting an improvement on the naive algorithm are susceptible to yield improvements.

Removing ε's and counting operations

As seen in Section 4.2.4, Pan proposes a method called trilinear aggregation to build degenerate tensors representing disjoint matrix products using as few operations as possible.

Combining this idea with the above algorithm Extract associated to Schönhage's τ -theorem proves to be a fruitful approach. The following are some new definitions and propositions.

Definition 8 (Homogeneity). Let t be a disjoint matrix product tensor of length h with the 3s variables (x 1,k ) k=1..s , (x 2,k ) k=1..s , (x 3,k ) k=1..s the 3s exponents (q 1,k ) k=1..s , (q 2,k ) k=1..s , (q 3,k ) k=1..s and the 3sh coefficients

(λ i,1,k ) i=1..h,k=1..s , (λ i,2,k ) i=1..h,k=1..s , (λ i,3,k ) i=1..h,k=1..s , such that t = h i=1 s k=1 λ i,1,k x 1,k ε q 1,k ⊗ s k=1 λ i,2,k x 2,k ε q 2,k ⊗ s k=1 λ i,3,k x 3,k ε q 3,k .
Then we say that t is s-homogeneous.

The key point is that each variable always appears with the same exponent. An example of a homogeneous tensor is presented in Equation (4.9). Theorem 3. Let t be a s-homogeneous tensor. Let µ

= (µ 1 ,• • • , µ s ), k such that s i=1 µ i = k and let t = Extract(t, k, µ). The k µ products k µ M k , N k , P k :=
t are of the form ε q t 1 , where q = 3 u=1 s i=1 q u,i µ i and t 1 is an ε-free tensor (i.e, q only depends on µ).

Proof. Let s l=1 λ ij ,1,l x 1,l ε q 1,l ⊗ s l=1 λ ij ,2,l x 2,l ε q 2,l ⊗ s l=1
λ ij ,3,l x 3,l ε q 3,l , j = 1..k be k terms (possibly repeated) of t. Their product appears in the expansion of t ⊗k :

k j=1 s l=1 λ ij ,1,l x 1,l ε q 1,l ⊗ s l=1 λ ij ,2,l x 2,l ε q 2,l ⊗ s l=1 λ ij ,3,l x 3,l ε q 3,l (4.20)
= 3 u=1   k j=1 s l=1 λ ij ,u,l x u,l ε q u,l   (4.21)
From the 3ks terms of

3 u=1   k j=1 s l=1 λ ij ,u,l x u,l ε q u,l
  , the only elementary tensors which appear in t are those of the shape x = x 1 ⊗ x 2 ⊗ x 3 with non-zero components of the form

x u ∼ =   s l=1   j∈J l λ ij ,u,l x u,l ε q u,l   ⊗µ l   (4.22) ∼ = ε qu,1µ1+•••+qu,sµs s l=1   j∈J l λ ij ,u,l x u,l   ⊗µ l . (4.23)
As all elementary tensors of t are of the shape (4.23), we can set q = 3 u=1 s i=1 q u,i µ i and factor t by ε q . The remaining term of this factorization is ε-free.

Theorem 3 means that we can benefit from the τ -theorem's improvements on Pan's degenerated tensors without the ε's . Indeed, in the proof of the τ -theorem, we were interested in upper bounding the border rank R (S k M k , N k , P k ). Then, we computed an asymptotic complexity for this upper bound as we knew that the tensor rank R (S k M k , N k , P k ) would be higher because of Lemma 8. However, in this case, our tensor is of the shape ε q t with t an ε-free tensor, so that by setting ε to 1, we can see that our bound on the border rank of t is also a bound on the tensor rank of

S k M k , N k , P k .
Remark 1. In Theorem 3, the terms of t all had the same degree in ε. If we look closely, there might be other terms of the multinomial expansion of Theorem 2 which happen to have the same degree. An improvement would be to keep, not only k µ M k , N k , P k but all terms corresponding to a k ν such that

3 u=1 s i=1 q u,i ν i = 3 u=1 s i=1 q u,i µ i .
One could then iterate the "τ -theorem algorithm" on the newly obtained disjoint matrix product. This might motivate a particular choice of the powers of ε in Pan's trilinear aggregating tables.

A natural question arises after this: how do we count the number r of elementary tensors in Extract(t, k, µ)? We will be able to compute k µ products of the type

s i=1 m µi i , s i=1 n µi i , s i=1 p µi i
(where the m i , n i , p i are the dimensions of each of the s disjoint matrix products) in parallel:

R 1, k µ , 1 ⊗ s i=1 m µi i , s i=1 n µi i , s i=1 p µi i ≤ r (4.24) so that ω ≤ log(r) log k µ s i=1 m µi i n µi i p µi i (4.25)
There is a way to compute r without computing t k .

Definition 9. Let t be a s-homogeneous tensor and let

t i = 3 u=1 s k=1 λ ij ,1,l x u,k ε λ u,k
be an elementary tensor of t. A 3s-tuple p(t i ) ∈ {0, 1} 3×s is called a pattern of

t i if ∀1 ≤ k ≤ s, 1 ≤ l ≤ 3, p(t i ) l,k = 0 if λ i,j,k = 0 1 otherwise Definition 10. The monomial M (p) associated with a pattern p ∈ {0, 1} 3×s is 3 u=1 s i=1 X p1,i i Y p2,i i Z p3,i i The polynomial P t1,•••,ts ∈ Z[X 1 ,• • • , X s , Y 1 ,• • • , Y s , Z 1 ,• • • , Z s ]
associated with a decomposition of an s-homogeneous tensor t as a sum of elementary tensors

t = s i=1 t i is defined as s i=1 M (p(t i )).
From now on we write P t for P t1,•••,ts as we fix a decomposition of t.

Theorem 4. Let t be a s-homogeneous tensor of length s. Let P = P t be its associated polynomial. Let k ∈ N and let µ be a s-tuple of sum k. Let Q be the polynomial obtained by keeping only the monomials of P k which are multiples

of R := s i=1 X µi i s i=1 Y µi i s i=1 Z µi i i.e., Q = P k -P k mod R .
Then the number of multiplications in the algorithm yielded by Theorem 3 for N and µ is equal to the sum of the coefficients of Q (the value of Q when all its variables are set to 1).

Proof. As there is one monomial for each term of t written as a homogeneous tensor, the products of monomials M 1 • • • M k in the expansion of P k t are in oneto-one correspondence with the tensor products of the elementary tensors in t ⊗k (see Equation (4.21)).

The operation P k -P k mod R keeps only the monomials of P k which are multiples of the monomial R. Thus, we need to show that these monomials are exactly the ones corresponding to the elementary tensors we kept in Equation (4.21). We define V

[1] = X, V [2] = Y and V [3] = Z.
We kept an elementary tensor of t ⊗k as in Equation (4.21) if and only if one of its sub-terms had the shape of Equation (4.22). This corresponds to keeping a monomial

M 1 • • • M k if and only if for each l, among the k monomials M 1 ,• • • , M k , for u ∈ {1, 2, 3}
, there are at least µ l which correspond to a tensor containing a non-zero term x u,l . The latter is equivalent to saying that for each l, there are at least

µ l monomials among M 1 ,• • • , M k which are divisible by V [u] l , which is equivalent to saying that the product M 1 • • • M k is divisible by R = ( s i=1 X µi i ) ( s i=1 Y µi i ) ( s i=1 Z µi i ).
Theorem 4 gives us an algorithm to compute the number of multiplications of the result of Extract(t, k, µ) for t a s-homogeneous tensor. We first decided to write an implementation of the Extract algorithm in the Maple symbolic computation software. Recall that this algorithm takes as input a tensor t such that

s i=1 m i , n i , p i t, a natural number k ∈ N, a vector µ = (µ 1 ,• • • , µ s ) such that s i=1 µ i = k,
and returns a tensor (i.e., an algorithm)

t such that k µ s i=1 m µi i , s i=1 n µi i , s i=1 p µi i t .
We illustrate the symbolic operations performed by this algorithm by going informally through an example with k = 2.

Input tensor: Suppose we start with the direct sum of two naive matrix multiplications, 2, 1, 2 ⊕ 1, 3, 1 :

t := 1 i=0 1 j=0 a i,0 ⊗ b 0,j ⊗ c i,j + 2 k=0 x 0,k ⊗ y k,0 ⊗ z 0,0
For this example we choose k = 2 and µ = (1, 1).

Computing t ⊗k : First, we compute t ⊗2 . The tensor product of variables is done in the following way. Indices are computed according to the rule already exposed just after Equation (4.6), which correspond to block matrix multiplication, and names are concatenated. For example when we take the product

1 i=0 a i,0 ⊗ 2 k=0
x 0,k , the variable a implicitly represents a matrix of size (2, 1) and x represents a matrix of size (1,[START_REF] Ahmed | Ahmed's integral: the maiden solution[END_REF]. Thus we create a variable ax representing a matrix of size (2 • 1, 1 • 3) = (2, 3) and the tensor product gives

1 i=0 2 k=0 ax i,k . Here is t ⊗2 : 1 i=0 1 j=0 1 i =0 1 j =0 aa (2•i+i ),0 ⊗ bb 0,(2•j+j ) ⊗ cc (2•i+i ),(2•j+j ) + 1 i=0 1 j=0 2 k =0 ax i,k ⊗ by k ,j ⊗ cz i,j + 2 k=0 1 i =0 1 j =0 xa i ,k ⊗ yb k,j ⊗ zc i ,j + 2 k=0 2 k =0 xx 0,(3•k+k ) ⊗ yy (3•k+k ),0 ⊗ zz 0,0
By grouping sums together, we can see the terms of the expansion in Equation (4.14) which we recall is:

s i=1 m i , n i , p i ⊗k ∼ = µ µ1+•••+µs=k k µ s i=1 m µi i , s i=1 n µi i , s i=1 p µi i
• The variable aa (respectively bb, cc) represents the vector µ = (2, 0); in other words the one in which all variables come from the first component of the direct sum 2, 1, 2 ⊕ 1, 3, 1 ;

• The variables ax and xa (respectively by and yb, cz and zc) represent the vector µ = (1, 1); in other words the one in which there is one variable (a) from the first component of the direct sum and one variable (x) from the second component of the direct sum;

• The variable xx (respectively yy, zz) represents the vector µ = (0, 2).

Extracting the term corresponding to µ and renaming variables to reflect parallel independent matrix multiplications: Suppose now that we want to extract the term corresponding to µ = (1, 1), which according to our expansion above (and in Equation (4.14)), should give us

2 1,1 = 2 independent copies of 2 • 1, 1 • 3, 2 • 1 = 2, 3, 2 .
We need to set to zero all variables which correspond to another term. We also rename the variables in a more convenient way, with a0, b0, c0 for the first copy of 2, 3, 2 and a1, b1, c1 for the second copy:

1 i=0 1 j=0 2 k0=0 a0 i,k0 ⊗ b0 k0,j ⊗ c0 i,j + 2 k=0 1 i0=0 1 j0=0 a1 i0,k ⊗ b1 k,j0 ⊗ c1 i0,j0
which is the output of our algorithm. Of course, the tensor t was rather uninteresting: it consisted of two copies of a naive matrix multiplication algorithm, and we can't expect the output to be better than several copies of a naive algorithm. The main goal of this example was to explain how variables are combined, renamed and extracted depending on µ.

Let us now consider a tensor t which is a degeneration of 2, 1, 2 ⊕ 1, 3, 1 , that is a tensor with ε's with strictly less elementary tensors than the naive algorithm, that is stricly less than 2

• 1 • 2 + 1 • 3 • 1 = 7 elementary tensors.
Let us apply the exact same transformations as above to t ⊗2 , and let us call the result t . The tensor t is a degeneration of 2 2, 3, 2 , i.e it computes 2 matrix products of the shape (2, 3, 2). The hope is to be able to pick sufficiently big values for k and µ, so that t features strictly less than 2 • (2 • 3 • 2) = 24 elementary tensors, or, algorithmically speaking, strictly less than 24 polynomial multiplications. Of course, in the case of our example, sticking to k = 2 makes finding such good values unlikely.

Symbolic computation systems are natural candidates for the implementation of such an algorithm: manipulating variables, subscripts and polynomials is one of the core use cases for such programs. We hence made a first attempt using Maple [START_REF] Michael | Maple 10 Programming Guide[END_REF]. In fact, the numerous commands available in Maple to extract information from data (normalizing polynomials to see if they are equal, looking at a specific coefficient, substituting in an algebraic way) were invaluable to make quick experiments which led to the observations of Section 4.3.2.

However, there were several reasons which ultimately led us to depart from a computer algebra system like Maple to a functional programming language like Ocaml.

The Maple implementation of the τ -theorem worked really well with concrete instances of tensors, i.e. when the number of elementary tensors was static. When we started to do substitutions under a binder " k i=1 " for tensors represented as sums of a variable number of elementary tensors, it became hard to tame the behavior of Maple, although this may very well be due to inexperience with the latter on our part. In particular, in Pan's tensors as presented in his book [100], there is often a clean presentation such as in Equation (4.9) which is then followed by a set of rewriting rules such as in Equation (4.10). Correctly rewriting these under a sum can prove tricky and we found it easier to write this ourself in Ocaml.

Related to the previous part is that we realized that some of what we wanted was similar to manipulations of an abstract syntax tree as it is done for instance in compilers; a statically typed functional programming language such as Ocaml is much more suited to that purpose than Maple. The abstract data structures used in Ocaml enables us to export our tensors in a modular way to a Symbolic Computation Software such as Maple or Sage, but also to C++ (in order to get actual runnable algorithms for parallel matrix computations) and to L A T E X.

Ocaml Implementation

We built a software tool in Ocaml enabling us to do symbolic manipulations on tensors.

Input

The input can be given this way. It consists of three parts: a list of tables, a set of constraints, and a list of the dimensions of the matrices. Each table is itself a collection of matrix multiplications (one on each row) representing Pan tables as described in Section 4.2.4. The constraints are a collection of rewriting rules of the form

x[i 1 ,• • • , i n ] := e where x is a variable with subscripts i 1 ,• • • , i n (usually n = 2 in
our examples, but we don't constrain it) and e is an arithmetic expression of variables. Subscripts are either static natural numbers such as 0 or 1, free variables such as k, l or they can be explicit substitution such as k := 0. The difference between x[0] = e and x[k := 0] = e is that in the first case, we will only substitute the syntactic expression x[0] by e but we won't substitute under the sum in 1 i=0 x[i]; whereas in the second case, we will substitute in both cases. When a variable has no subscript, it can also be written

x instead of x[].
Here is Schönhage's algorithm from Section 4.2.4 for instance. The table T1 corresponds to Equation (4.8) and the table T2 corresponds to

Table T1 : Line l11 : Seq([x[i,0];y[0,k];z[k,i] * eps^2],i = 0 .. (m-1),k = 0..(p-1)) Line l12 : Seq([u[0,k,i] * eps;v[k,i,0]*eps;w[0,0]],i = 0.. m-1,k = 0..p-1) Table T2 : Line l21 : [sum(x[dd,0],dd=0..(m-1));sum(y[0,ee],ee=0..(p-1));-w[0,0]] ;; Constraints : u[x:=0,k:=0,i] = 0 v[k,i:=0,x:=0] = 0 u[x:=0,k,i:=0] = -sum(u[0,k,dd],dd=1..m-1) v[k:=0,i,x:=0] = -sum(v[ee,i,0],ee=1..p-1) x[i,0] = a1[i,0] x[dd,0] = a1[dd,0] y[0,k] = b1[0,k] y[0,ee] = b1[0,ee] z[k,i] = c1[i,k] u[0,k,i] = a2[0,(k-1)*(m-1)+(i-1)] v[k,i,0] = b2[(k-1)*(m-1)+(i-1),0] w[0,0] = c2[0,0] eps = 1 ;; Spaces : a1 b1 c1 (m,1,p) a2 b2 c2 (1,(m-1)*(p-1),1) ;;
Notice that eps = 1 is part of the constraints, where eps is the variable ε. This is justified by our Theorem 3. This input format is not very restrictive, and in fact we can input more classical matrix multiplication algorithms. Here is Strassen's algorithm for instance: The first seven lines, shown as one-line Pan tables represent each of the seven elementary tensors in Strassen's algorithm:

Table T1 : Line l11 : [a[1,2] -a[2,2];b[2,1] + b[2,2]; z[1,1]] Table T2 : Line l21 : [a[2,1] -a[1,1];b[1,2] + b[1,1]; z[2,2]] Table T3 : Line l31 : [ a[1,1];b[1,2] -b[2,2];z[2,1] + z[2,2]] Table T4 : Line l41 : [a[2,2];b[2,1] -b[1,1];z[1,2] + z[1,1]] Table T5 : Line l51 : [a[2,1] + a[2,2];b[1,1];z[1,2] -z[2,2]] Table T6 : Line l61 : [a[1,2] + a[1,1];b[2,2];z[2,1] -z[1,1]] Table T7 : Line l71 : [a[1,1]+a[2,2];b[1,1]+b[2,2];z[1,1]+z[2,
(a 1,2 -a 2,2 ) ⊗ (b 2,1 + b 2,2 ) ⊗ z 1,1 + (a 2,1 -a 1,1 ) ⊗ (b 1,2 + b 1,1 ) ⊗ z 2,2 + a 1,1 ⊗ (b 1,2 -b 2,2 ) ⊗ (z 2,1 + z 2,2 )+ a 2,2 ⊗ (b 2,1 -b 1,1 ) ⊗ (z 1,2 + z 1,1 )+ (a 2,1 + a 2,2 ) ⊗ b 1,1 ⊗ (z 1,2 -z 2,2 )+ (a 1,2 + a 1,1 ) ⊗ b 2,2 ⊗ (z 2,1 -z 1,1 )+ (a 1,1 + a 2,2 ) ⊗ (b 1,1 + b 2,2 ) ⊗ (z 1,1 + z 2,2 )
Finally, the "Spaces" line expresses that a is a matrix variable of size 2, 2, as well as b and c.

Here are the options currently permitted by our tool, which are showed when it is called with the option --help:

Tool to parse tensors in the shape of Pan's tables -f Name of file to parse -pow Specifies the power at which we want to compute the tensor -output {plain|tautheorem} Type of output to compute: plain tensor (nothing set to zero) or tau-theorem -style {latex|cpp|maple|sage} Style of output -raw Set to true if your tensor is not a direct sum of matrix products -constraints Should constraints specified in <file.in> be applied? Defaults to true -laser Special option for the laser method, sets all parameters right --help Display this list of options

The -raw option decides whether we want to interpret our tensor in terms of matrix products, which is the case for all examples presented here, or if we just want to see it as a "raw" tensor, meaning a tensor not directly encoding a matrix product. This would be the case for Strassen's initial tensor for the "laser" method [START_REF] Volker Strassen | Relative bilinear complexity and matrix multiplication[END_REF], which does not represent a direct sum of matrix products. We chose not to detail here our support for the laser method, as it is still too rudimentary for now.

Style of Output

Maple: Any tensor can be outputted to a format which will accepted as the input of a Symbolic Computation Software format. For now, only the syntax of Maple is supported, but little effort would be needed to support others such as Sage.

C++: If a tensor represents a matrix product or a direct sum of matrix products in a way that is made explicit by the spaces section, and if it doesn't have any extra variables left (such as ε in our previous examples), it can be extracted to a C++ program computing the corresponding matrix multiplications. Of course, the latter are just necessary conditions: the correctness of the algorithm is contingent on the correctness of the tensor.

Latex: Tensors can also be exported to Latex. This feature was used to write down the examples in Section 4.4.1.

Power

One can compute powers of tensors using the -pow option. Composition will be done differently depending on whether the raw option is set to :

• true: To compose two variables x i1,..,i k and y j1,..,j l we concatenate both their names and indices into xy i1,..,i k ,j1,..j l ;

• false: we compose variables in a way that reflects block matrix product, as shown in Section 4.1.2. This assumes that each variable has exactly two indices and has an associated matrix size which has been specified to the program, in the "spaces" part of the input file.

Output

The -output option can be either

• plain: in which case the computed power of the tensor will be exported as is (with no terms set to zero, possibly with ε's remaining, and with no attempt to rename variables)

• tautheorem: this will only be valid if raw is false, and will apply the algorithm Extract described in Section 4.4.1 to the power of the input tensor with a µ which maximizes k µ .

Constraints

The constraints option defaults to true, but can be set to false if the user does not want to take constraints specified in the input into account.

Exact Complexity of Some Matrix Product Algorithms

In this section, we explore exactly how many operations (scalar multiplications and additions) it takes to compute matrix products using some of the algorithms explored.

Methodology

Contrary to the number of multiplications implied by an algorithm represented by a tensor, which is straightforward (it is exactly the number of elementary tensors, which our program outputs), the number of additions is more tricky to estimate, at least given a tensor depending on parameters. Given an elementary tensor (a

1,1 + a 1,2 ) ⊗ (b 2,1 + b 2,2 ) ⊗ (c 1,1 + c 2,2 )
we could say there is one addition to compute A := a 1,1 + a 1,2 , one for B := b 2,1 + b 2,2 and then one to add AB to c 1,1 and another to add it to c 2,2 . If we do this naively across all elementary tensors, we might miss common terms: in

a 1,1 ⊗ (b 2,1 + b 2,2 ) ⊗ c 1,1 + a 1,2 ⊗ (b 2,1 + b 2,2 ) ⊗ c 2,1
, we don't want to count the addition b 2,1 + b 2,2 twice as it could easily be done only once. We tried two approaches for counting additions:

• The first one was to use the optimize function from the codegen library in Maple. Given a set of algebraic expressions S to compute, it gives an "optimized computation sequence", meaning a sequence of steps which will obtain all the elements of S in a hopefully small number of steps by identifying linear common subexpressions. This works well on concrete tensors, as for example when parameters such as m and p are replaced with small numeric values in the case of Schönhage's tensor in Equation (4.9);

• The second one was to statically estimate in Ocaml the number of additions, memorizing each (symbolic) encountered expression and putting it in a hashtable so that whenever we saw it again, we didn't count it again.

The downside of the first method is that it only works on relatively small tensors, so we took the minimum of both values when we had access to both, and only the Ocaml estimate otherwise.

Remark: Both methods are actually slightly overestimating the number of additions, because they are counting each variable representing the output matrix as one addition, whereas the first value stored at the initialization of the variable should not count as an addition. For example, our method applied to the naive algorithm for m, n, p gives mnp additions instead of mn(p -1), and it gives 22 additions for Strassen's algorithm instead of 18.

Schönhage's tensor

Recall Schönhage's tensor, already presented in Equation (4.9) as an example for Pan's tables:

(m-1) i=0 (p-1) k=0 (u 0,k,i • ε + x i,0 ) • (v k,i,0 • ε + y 0,k ) • (w 0,0 + z k,i • ε 2 )+ (-1) • (m-1) dd=0 x dd,0 • (p-1) ee=0 y 0,ee • w 0,0
It can be seen to be a degeneration of m, 1, p ⊕ 1, (m -1)(p -1), 1 in mp + 1 (polynomial) multiplications. In Table A.1, we computed the 25 best values we get for ω with concrete values for m and p between 2 and 15, and by taking the tensor to a power between 2 and 4 and applying our method. For example, the last row of this table gives a "practical ω" of 2.774065. This means if one makes multiple copies of the 6 49, 1764, 64 that we obtained and packs them into a rectangular matrix product, one can build an algorithm with the corresponding ω for that rectangular matrix product. Alternatively, one can directly compute 6 matrix products of size 49, 1764, 64 and get a speedup of 6 * 49 * 1764 * 64 10308816 3.22 in terms of the number of multiplications.

A tensor from Pan

We wrote an input file for a tensor presented by Pan [100], which can be seen in Appendix A.2. This tensor is a degeneration of the following direct sum of matrix products:

m -1, 1, 2p -2 ⊕ 2, p -1, m -1 ⊕ p -1, 2m -2, 1 (4.26)
and it consists in 2m(p + 1) (polynomial) multiplications. We found a mistake in the description of this tensor and had to incorporate some correcting (scalar) factors for it to be a degeneration of the disjoint matrix multiplications described in Equation (4.26). We detected this thanks to our implementation by outputting a term to Maple and trying to match it with the claimed matrix products. Once we did the correction, on top of checking the terms in Maple again, we used a testing function in Ocaml to multiply random matrices of integers to make sure we had gotten it right. The result is the tensor presented in Listing A.1, and which also comes with our distributed code.

Conclusion

In this part, we have presented the result of our experiments with programming theoretical matrix multiplication algorithms in Maple and Ocaml. We wrote a constructive proof of Schönhage's τ -theorem. We discovered that starting with elements of a new class of s-homogeneous degenerate tensors, we could build exact (ε-free) algorithms computing several matrix products of reasonable sizes in parallel, with significantly less multiplications than the naive algorithm and in some cases than Strassen's algorithm. We built an Ocaml tool which allows us to manipulate matrix multiplication tensors symbolically, and output them in various formats. The data structure of Pan tables is general enough that we were able to encode every matrix multiplication tensor we could get our hands on into it.

It is our hope that this tool can be of help to discover more properties of matrix multiplication tensors in the future. We have started looking at Strassen's Laser method [START_REF] Volker Strassen | Relative bilinear complexity and matrix multiplication[END_REF] and the Coppersmith-Winograd algorithm [START_REF] Coppersmith | Matrix multiplication via arithmetic progressions[END_REF] which are both crucial for more recent improvements by Williams [START_REF] Williams | Multiplying matrices faster than Coppersmith-Winograd[END_REF] and Le Gall [START_REF] François | Powers of tensors and fast matrix multiplication[END_REF]. We do not yet have results on the last two, beyond the fact that they can be encoded as Pan tables and so they can be made into actual, concrete algorithms. To the best of our observations so far, it seems like the number of multiplications needed by these two types of algorithms for reasonable matrix sizes is not competitive, consistently with what is usually claimed. Nevertheless there could be a pedagogical utility to being able to manipulate such tensors.

It turns out that the ε-free independent parallel matrix product algorithms from Theorem 3 do not perform well (about 10 times worse than the naive multiplication algorithm), at least in the way we extract them to C++. Although this is disappointing, it need not be the end of the story. First, the way we output code to C++ is probably far from optimal. Secondly, we hope that our tool and findings can be used to experiment and build new ε-free matrix product algorithms, just like we used the Combine library to find new tilings in Chapter 3.

Chapter 5

The result of a computation can provide a varying degree of confidence in a mathematical fact. On the skeptical side, the output of a calculator is not enough to be convinced of a numerical fact: Ramanujan famously conjectured that exp(π When possible, one would wish to obtain as much confidence as possible in outputs of computations while still benefiting from the considerable engineering feats of scientific software. In particular, this would not require to reinvent the wheel and implement a suboptimal algorithm in order to obtain better guarantees that the result was obtained in a sound way.

Proof assistants are programs in which one can manually define mathematical objects, make mathematical statements about them and write proofs of such statements. They emphasize soundness: a great amount of care is taken so that only valid proofs are accepted. This often comes at the expense of efficiency: even though they usually support doing computations, the latter are typically slower than equivalent code in a powerful programming language.

This power of computer algebra software comes with different assumptions. Maple will gladly simplify x x = 1 without knowing what type of object x is or whether it can be 0 when that makes sense, and thus it is not clear what correctness might mean in most settings. Moreover, Section 12.3 will establish that regarding numerical integral computations, neither correctness (the outputted value is within some guaranteed precision of the actual value) nor termination (the program finished computing in finite time) are guaranteed.

One way to obtain more trustworthy computations would be to write programs inside a proof assistant, prove them correct, and execute them to get an indisputable answer. However, on top of the efficiency problem already mentioned, proving modern scientific software correct is a huge endeavor. A formal proof of the correctness of GMP's square root algorithm [START_REF] Bertot | A proof of gmp square root[END_REF] took 13, 000 lines of Coq, and this is a small routine in a huge piece of software. Not all problems have solutions which are as trivial to check as the correctness of a Sudoku board. However, many problems lie in a larger class: those for which checking a solution is easier than finding a solution. This distinction is natural if one thinks about complexity theory, whose most famous challenge consists in comparing the class P of problems for which there exists a deterministic algorithm which can find a solution in polynomial time to the larger class N P of problems for which there exists a deterministic algorithm which can check a solution in polynomial time.

Analogously, a hybrid solution for problems which lie in the right class of easily checkable solutions is to divide the work between an oracle and a verifier. A tongue-in-cheek example of this paradigm is the "twittersort" algorithm1 for sorting arrays. When the user inputs an array of numbers to sort, the twittersort program publishes the array on the Twitter social network, and waits for someone to respond with a sorted array. If and when someone does, it checks that the numbers in the array are sorted and correspond to the input array, and if it is correct it accepts it as the final result. The oracle produces a result and a certificate.

In this paradigm, the only part which needs to be carefully proved is the specification of the work of the verifier, i.e., the verifier accepts a result if and only if it is valid. If the verifier function is not too complex, this can be done in a proof assistant so as to obtain high confidence in our results. The oracle can be implemented however one wishes, in any programming language, without any requirement of termination, completeness or correctness. This laxness is fortunate in the case of scientific software for which none of these is guaranteed.

Harrison and Théry [START_REF] Harrison | A skeptic's approach to combining HOL and Maple[END_REF] develop the idea of a skeptic's approach in the particular context of Computer Algebra Software (CAS). They implement an interface between the Maple CAS and the HOL proof assistant which can call the simplify and factor Maple commands.

This separation of concerns in the paradigm of certificate-based proofs pro-vides an ideal combination between the power and speed of scientific computing software and the uncompromising mathematical rigor of formal proof systems. One uses each system for what it does best and avoids it for what it does worst. The technique of creative telescoping is a class of algorithms which take as input recurrence relations satisfied by a summand, and produce certificates from which one can deduce recurrence relations satisfied by a sum of this summand [START_REF] Chyzak | The ABC of Creative Telescoping-Algorithms, Bounds, Complexity[END_REF]. They output a formal identity which is often easy to check, and then getting the recurrence boils down to summing a telescoping sum, that is, a sum of the shape k (u k+1 -u k ). The easy and mechanical aspect of checking these outputted certificates may lead one to think that they constitute proofs of recurrence relations. A natural step would be to make these proofs into formal proofs, by using a proof assistant to fully verify them. In a Maple sheet [START_REF] Salvy | An Algolib-aided version of Apéry's proof of the irrationality of ζ(3)[END_REF], Salvy uses creative telescoping algorithms and in particular their implementation in the Algolib [START_REF] Chyzak | Algolib[END_REF] Maple library to present a creative-telescoping-based proof of the following theorem in Number Theory: Theorem 5 (Apéry, 1978).

The constant ζ(3) := ∞ i=1 1 i 3 is irrational.
This result was the first dent in the problem of the irrationality of the evaluation of the Riemann ζ function at odd positive integers. As of today, this problem remains a long-standing challenge of number theory. Although Rivoal [START_REF] Rivoal | Propriétés diophantiennes de la fonction zêta de Riemann aux entiers impairs[END_REF] and Zudilin [START_REF] Vadim | One of the numbers ζ(5), ζ(7), ζ(9), ζ(11) is irrational[END_REF] showed that at least one of the numbers ζ( 5 3) is the only one known to be irrational.

In the present part, we describe a formal proof of Theorem 5 inside the Coq proof assistant, using the Mathematical Components libraries [1]. This formalization follows the structure of Apéry's original proof. However, we replace the manual verification of the recurrence relations proposed by Apéry with an automatic discovery of these equations, using the symbolic computations proposed by Salvy. For this purpose, we use Maple packages to perform calculations outside the proof assistant, and we verify a posteriori the resulting claims inside Coq. By combining these verified results with additional formal developments, we obtain a complete, constructive formal proof of Theorem 5. We report on the implementation of this cooperation between a computer algebra system and a proof assistant.

We describe in detail the formalization of an upper bound on the asymptotic behavior of lcm(1, ..., n), the least common multiple of the integers from 1 to n, a part of the proof which was missing in a published paper [START_REF] Chyzak | A Computer-Algebra-Based Formal Proof of the Irrationality of ζ(3)[END_REF].

The rest of the present part is organized as follows. We first describe the background formal theories used in our development (Section 7). We then outline the proof of Theorem 5 (Chapter 6). We summarize the algorithms used in the Maple session, the data this session produces and the way this data can be used in formal proofs (Section 8.1). We then describe the proof of the consequences of Apéry's recurrence (Section 8.2). Finally, we present an elementary proof of the bound on the asymptotic behavior of the sequence lcm(1, ..., n), which is used in this irrationality proof (Section 8.3), before concluding (Chapter 9).

Chapter 6

History and outline of Apéry's theorem Van der Poorten [START_REF] Van Der Poorten | A proof that Euler missed: Apéry's proof of the irrationality of ζ(3)[END_REF] reports that Apéry's announcement of the result of Theorem 5 was at first met with wide skepticism. His obscure presentation featured "a sequence of unlikely assertions" without proofs, not the least of which was an enigmatic recurrence (Lemma 13) satisfied by two sequences (a n ) and (b n ). It took two months of collaboration between Cohen [START_REF] Cohen | Démonstration de l'irrationalité de ζ (3)(d'apres apéry)[END_REF], Lenstra, and Van der Poorten, with the help of Zagier, to obtain a thorough proof of Theorem 5.

There exists other proofs of Apéry's theorem, such as the one by Beukers [START_REF] Beukers | A note on the irrationality of ζ(2) and ζ(3)[END_REF]. As explained for example in Fischler' survey [START_REF] Fischler | Irrationalité de valeurs de zêta (d'après Apéry, Rivoal[END_REF], all these proofs share a common structure. They rely on the asymptotic behavior of the sequence n , the least common multiple of integers between 1 and n, and they proceed by exhibiting two sequences of rational numbers a n and b n , which have the following properties:

1. For a sufficiently large n:

a n ∈ Z and 2 3 n b n ∈ Z; 2. The sequence δ n = a n ζ(3) -b n is such that: lim sup n→∞ |2δ n | 1 n ≤ ( √ 2 -1) 4 ;
3. For an infinite number of values n, δ n = 0.

Altogether, these properties entail the irrationality of ζ(3). Indeed, if we suppose that there exists p, q ∈ Z such that ζ(3) = p q , then 2q 3 n δ n is an integer when n is large enough. One variant of the Prime Number theorem [START_REF] Harrison | Formalizing an analytic proof of the Prime Number Theorem[END_REF] states that n = e n (1+o(1)) and since ( √ 2 -1) 4 e 3 < 1, the sequence 2q 3 n δ n has a zero limit, which contradicts the third property. Actually, the Prime Number theorem can be replaced by a weaker estimation of the asymptotic behavior of n , that can be obtained by more elementary means.

Lemma 11. Let n be the least common multiple of integers 1 . . . n, then

n = O(3 n ).
Since we still have ( √ 2 -1) 4 3 3 < 1, this observation [START_REF] Hanson | On the product of the primes[END_REF][START_REF] Feng | A simple elementary proof for the inequality d n < 3 n[END_REF] is enough to conclude. Section 8.3 discusses the formal proof of Lemma 11, an ingredient which was missing at the time of writing the previous report on this work [START_REF] Chyzak | A Computer-Algebra-Based Formal Proof of the Irrationality of ζ(3)[END_REF].

In our formal proof, we consider the pair of sequences proposed by Apéry in his proof [START_REF] Apéry | Irrationalité de ζ(2) et ζ(3)[END_REF][START_REF] Van Der Poorten | A proof that Euler missed: Apéry's proof of the irrationality of ζ(3)[END_REF]:

a n = n k=0 n k 2 n+k k 2 , b n = a n z n + n k=1 k m=1 (-1) m+1 n k 2 n+k k 2 2m 3 n m n+m m (6.1)
where z n denotes n m=1 1 m 3 , as already used in Proposition 1. By definition, a n is a positive integer for any n ∈ N. The integrality of 2 3 n b n is not as straightforward, but rather easy to see as well: each summand in the double sum defining b n has a denominator that divides 2 3 n . Indeed, after a suitable re-organization in the expression of the summand, using standard properties of binomial coefficients, this follows easily from the following slightly less standard property: Lemma 12. For any integers i, j, n such that 1 ≤ j ≤ i ≤ n, j i j divides n .

Proof. Using Lemma 14, observe that for any prime p, the p-valuation of j i j is smaller than the one of n .

The rest of the proof is a study of the sequence δ n = a n ζ(3) -b n . It is not difficult to see that δ n tends to zero, from the formulas defining the sequences a and b, but we also need to prove that it does so fast enough to compensate for 3 n , while being positive. In his original proof, Apéry derived the latter facts by combining the definitions of the sequences a and b with the study of the mysterious recurrence relation (6.2). Indeed, he made the surprising claim that Lemma 13 holds: Lemma 13. For n ≥ 0, the sequences (a n ) n∈N and (b n ) n∈N satisfy the same second-order recurrence:

(n + 2) 3 y n+2 -(17n 2 + 51n + 39)(2n + 3)y n+1 + (n + 1) 3 y n = 0. (6.2) Equation 6.2 is a typical example of a linear recurrence equation with polynomial coefficients and standard techniques [START_REF] Salvy | An Algolib-aided version of Apéry's proof of the irrationality of ζ(3)[END_REF][START_REF] Van Der Poorten | A proof that Euler missed: Apéry's proof of the irrationality of ζ(3)[END_REF] can be used to study the asymptotic behavior of its solutions. Using this recurrence and the initial conditions satisfied by a and b, one can thus obtain the last two properties of our criterion, and conclude with the irrationality of ζ(3). Our formal proof relies on an elementary version of this asymptotic study, mostly based on variations on the presentation of van der Poorten [START_REF] Van Der Poorten | A proof that Euler missed: Apéry's proof of the irrationality of ζ(3)[END_REF]. We detail this part of the proof in Section 8.2.

Using only Equation 6.2, and sufficiently many initial conditions, it would not be easy to obtain the first property of our criterion, about the integrality of a n and b n for a large enough n. In fact, it would also be difficult to prove that the sequence δ tends to zero: we would only know that it has a finite limit, and how fast the convergence is. By contrast, it is fairly easy to obtain these facts from the explicit Formulas 6.1.

The proof of Lemma 13 was by far the most difficult part in Apéry's original exposition. In his report [START_REF] Van Der Poorten | A proof that Euler missed: Apéry's proof of the irrationality of ζ(3)[END_REF], van der Poorten describes how he, with other colleagues, devoted significant efforts to this verification after having attended the talk in which Apéry exposed his result for the first time. The proof of Lemma 13 boils down to a routine calculation using the two auxiliary sequences U n,k and V n,k , themselves defined in terms of

λ n,k = n k 2 n+k k 2 (with λ n,k = 0 if k < 0 or k > n): U n,k = 4(2n + 1)(k(2k + 1) -(2n + 1) 2 )λ n,k , V n,k = U n,k n m=1 1 m 3 + k m=1 (-1) m-1 2m 3 n m n+m m + 5(2n + 1)k(-1) k-1 n(n + 1) n k n + k k
The key idea is to compute telescoping sums for U and V . For instance, we have:

U n,k -U n,k-1 = (n + 1) 3 λ n+1,k -(34n 3 + 51n 2 + 27n + 5)λ n,k + n 3 λ n-1,k (6.3)
Summing Equation 6.3 on k shows that the sequence a satisfies the recurrence relation of Lemma 13. A similar calculation proves the analogue for b, using telescoping sums of the sequence V . Not only is the statement of Formula 6.2 difficult to discover: even when this recurrence is given, finding the suitable auxiliary sequences U and V by hand is a difficult task. Moreover, there is no other known way of proving Lemma 13 than by exhibiting this nature of certificates. Fortunately, the sequences a and b belong in fact to a class of objects well known to combinatorialists and computer-algebraists, called ∂-finite sequences. Following seminal work of Zeilberger's [START_REF] Zeilberger | A holonomic systems approach to special functions identities[END_REF], algorithms have been designed and implemented in computer-algebra systems, which are able to obtain linear recurrences for these sequences. For instance the Maple package Mgfun (distributed as part of the Algolib [START_REF] Chyzak | Algolib[END_REF] library) implements these algorithms, among others. Salvy's worksheet [START_REF] Salvy | An Algolib-aided version of Apéry's proof of the irrationality of ζ(3)[END_REF] is based on this implementation and follows Apéry's original method but interlaces Maple calculations with human-written parts. In particular, this worksheet illustrates how parts of this proof, including the discovery of Apéry's mysterious recurrence, can be performed by symbolic computations. Our formal proof of Lemma 13 follows an approach similar to the one of Salvy. It is based on calculations performed using the Algolib [START_REF] Chyzak | Algolib[END_REF] library, and certified a posteriori. This part of the formal proof is discussed in Section 8.1.

Chapter 7 Preliminaries

This section provides some hints about the representation of the different natures of numbers at stake in this proof in the libraries backing the formal development. It also describes a few extensions devised for these libraries and fixes some notations used throughout the present part. Most of the material presented here is related to the Mathematical Components libraries [1] [52].

Integers

In Coq, the set N of natural numbers is usually represented by the type nat:

Inductive nat := O | S : nat -> nat.
This type is defined in a prelude1 library, which is automatically imported by any Coq session. It models the elements of N using a unary representation: Coq's parser reads the number 2 as the term S (S O). This inductive type comes with the usual recurrence scheme on natural numbers. This is convenient for defining elementary functions on natural numbers, like comparison or arithmetical operations, and for developing their associated theory. However, the resulting programs are usually very naive and inefficient implementations, which should only be evaluated for the purpose of small scale computations.

The set Z of integers can be represented by gluing together two copies of type nat, which provides a signed unary representation of integers:

Inductive int : Set := Posz of nat | Negz of nat.
If the term n : nat represents the natural number n ∈ N, then the term (Posz n) : int represents the integer n ∈ Z and the term (Negz n) : int represents the integer -(n + 1) ∈ Z. In particular, the constructor Posz : nat -> int implements the embedding of type nat into type int, which is invisible on paper because it is just the inclusion N ⊂ Z. In order to mimic the mathematical practice, the constant Posz is declared as a coercion, which means in particular that unless otherwise specified, this function is hidden from the terms displayed by Coq to the user (in the current goal, in answers to search queries,. . . ) and automatically inserted.

The Mathematical Components libraries provide formal definitions of a few concepts and results from elementary number theory, defined on the type nat. For instance, they provide the theory of Euclidean division, a boolean primality test, the elementary properties of the factorial function, of binomial coefficients,. . . In the rest of the present part, we use the standard mathematical notations n! and n m for the corresponding formal definition of factorial and binomial coefficients. These libraries also define the p-valuation v p (n) of a number n: if p is a prime number, v p (n) corresponds to the exponent of p in the prime decomposition of n. However, we had to extend the available formal theory with a few extra standard results like the formula giving the p-valuation of factorials: Lemma 14 (Legendre's Formula). For any n ∈ N and for any prime number p:

v p (n!) = log p n i=1 n p i .
Incidentally, the formal version of this formula is a typical example of the slight variations one may introduce in a mathematical statement, in order to come up with a formal sentence which is not only correct and faithful to the original mathematical result, but also a tool which is easy to use in subsequent formal proofs. In the formal library, Lemma 14 is in fact stated as: For any prime p and any j, n ∈ N, with n < p j+1 , v p (n!) = Adding an extra variable to generalize the upper bound of the sum is a better option because it will ease unification when this formula is applied or used for rewriting. For example, suppose that we want to prove that v p (n!) = log p n +1 i=1 n p i . Before applying Lemma 14, we would first need to break down the sum into two parts, and then prove that the second term n p log p n +1 is zero. However, we can directly apply the lemma in Equation (7.1): Coq will unify j with log p n + 1 and ask us to prove n < p ( log p n +1)+1 which is easy.

Moreover, we do not really need to introduce logarithms to state this lemma: indeed, log p n is used to denote the largest power of p smaller than n. For this purpose, we could use the function trunc_log : nat -> nat -> nat provided by the Mathematical Components libraries, which computes the greatest exponent α such that n α ≤ m, in other words log n m . Better yet, since the summand is zero when the index i exceeds this value, we can simplify the side condition on the extra variable and require only that n < p j+1 . Finally, although the fraction in the original statement of Lemma 14 may suggest that rational numbers play a role here, n m is in fact exactly the quotient of the Euclidean division of n by m. In the present part, for n, m ∈ N and m non-zero, we thus write n m for the quotient of the Euclidean division of n by m.

The basic theory of binomial coefficients present in the Mathematical

Components libraries describes their role in elementary enumerative combinatorics. However, when viewing binomial coefficients as a sequence which is a certain solution of a recurrence system, it becomes natural to extend their domain of definition to integers: we thus developed a small library about these generalized binomial coefficients. We also needed to extend these libraries with properties of multinomial coefficients. We used the formula

l i=1 k1+•••+ki ki from Equation (1.
2) in our formal definition, as it provides for free the fact that multinomial coefficients are non-negative integers, and proved the other characterizations, including the generalized Newton formula of Theorem 1 describing the expansion of (x

1 +• • • + x l ) n .

Sharing theories and notations

The Mathematical Components libraries feature a hierarchy of algebraic structures [START_REF] Garillot | Packaging Mathematical Structures[END_REF], which organizes a corpus of theories and notations shared by all the instances of the structures. This hierarchy implements inheritance, so that for example Z seen as a ring inherits from the properties of Z as an additive group. The hierarchy also implements sharing, so that the properties of the + and * operations of rings can be used in any instance of ring. Inheritance and sharing are made possible in Coq by mechanisms called coercions and canonical structures [START_REF] Mahboubi | Canonical Structures for the working Coq user[END_REF]. Each structure in the hierarchy is modeled by a structure called a dependent record, which packages a carrier type with some operations on this type and with requirements on these operations. For example, these structures are all discrete, which means that they require a boolean equality test. Part of this hierarchy deals with ordered structures [START_REF] Cohen | Formal proofs in real algebraic geometry: from ordered fields to quantifier elimination[END_REF], which means that they require an additional binary boolean predicate. This predicate has an infix notation ≤ and it is used to model an order relation which is (possibly partial and) compatible with the algebraic laws of the structure. Some structures in the hierarchy feature operations that make sense only on a subset of the elements of the carrier type: for instance, by definition, only the units of a ring have an inverse. In the dependent type theory implemented by Coq, it would be possible to use a dependent pair in order to model the source type of such an inverse operation. Instead, as a rule of thumb, the signature of a given structure avoids using rich types as the source types of their operations but rather "curry" the specification. For instance, the source type of the inverse operation in the structure of ring with units is the carrier type itself, but the signature of this structure also has a boolean predicate, which selects the units in this carrier type. The inverse operation has a default behavior outside units and the equations of the theory that involve inverses are typically guarded with invertibility conditions. Hence although the expression x^-1 * x makes sense for any term x of an instance of ring with units, it can be rewritten to 1 only when x is known to be invertible.

In order to equip a given type with a certain structure, one should package this type with enough operations and properties, following the signature prescribed by the structure. For example, the type int introduced in Section 7.1 is an instance of ordered integral domain. Unlike the case of Posz, the canonical embedding of type int in any instance of ring cannot be declared as a coercion, and eluded in formal statements. Its formal definition hence comes with a generic postfix notation _%:~R, used to cast an integer as an element of another ring. As for rational numbers, they are represented using a dependent pair of a pair (p, q) of integers, together with a proof that the corresponding fraction is normalized: p and q integers have to be coprime and q should be positive. By Hedberg's theorem [START_REF] Hedberg | A coherence theorem for Martin-Löf's type theory[END_REF], this proof can be made irrelevant by using a boolean predicate in order to express this normal form condition. Thus it is possible to compare rational numbers by comparing only the first components of the dependent pairs, and therefore to implement boolean operations of comparison. As a result, the type rat modeling the set Q of rational numbers is equipped with a structure of ordered field.

The Mathematical Components libraries also include a construction of the algebraic closure Q of the field of rational numbers (more in Section 7.3), resulting in the definition of a type algC which is an instance of a (partially) ordered, algebraically closed field [START_REF] Gonthier | A Machine-Checked Proof of the Odd Order Theorem[END_REF]. As a result, the statement 0 < x * y makes sense in Coq whether x and y are of type int, rat or even algC, because all these types share the notations of the signature of ordered rings. Moreover, in all these cases, this statement can be turned into 0 < y * x using the same lemma:

Lemma mulrC (R : comRingType) (r1 r2 : R) : r1 * r2 = r2 * r1.

because all these types are instances of the structure of commutative ring.

Algebraic numbers, real numbers

Almost all the irrational numbers involved in the proof are real algebraic numbers, and more precisely, they are of the form r 1 n for some rational number r (and with n ∈ N). These numbers are involved in inequalities expressing signs and estimations. It might come as a surprise that we use the type algC of algebraic (complex) numbers mentioned in Section 7.2 to cast these quantities, although we do not actually need complex numbers. But this choice proved convenient due to the fact that the type algC features both a definition of n-th roots, and a clever choice of partial order. Indeed, although Q cannot be ordered as a field, it is equipped with a binary relation, denoted ≤, which coincides with their real order relation on Q ∩ R:

∀x, y ∈ Q, x ≤ y ⇔ y -x ∈ R + .
In particular, for any z ∈ Q:

0 ≤ z ⇔ z ∈ R + and z ≤ 0 ⇔ z ∈ R -.
Moreover, the type algC comes with a function n.-root : algC -> algC, for any (n : nat), which computes the n-th (complex) root of its input which has a minimal non-negative argument. Crucially, when (z : algC) represents a nonnegative real number, (n.-root z) coincides with the definition of the real n-th root:

Lemma rootC_ge0 (n : nat) (x : algC) : n > 0 -> (0 <= n.-root x) = (0 <= x).

The shape of Lemma rootC_ge0 is typical of the style pervasive in the Mathematical Components libraries, where equivalences between decidable statements are stated as boolean equalities. It expresses that for an algebraic number x ∈ Q, x

The one notable place at which we need to resort to a larger set of real numbers is the definition of the number ζ(3), if only because as of today, it is not even known whether ζ(3) is algebraic or transcendant. This number is defined as the limit of the sums z n = n m=1 1 m 3 , so we start our formal study by establishing the existence of this limit. In fact we formalize the part of this proof concerned with asymptotic properties using Cauchy sequences. A Cauchy real is a sequence (x n ) n∈N ∈ Q N together with a modulus of convergence m x such that if ε ∈ Q + * , any two elements of index greater than m x (ε) are separated at most by ε. Two Cauchy reals x and y are equal, written x = y, if eventually |x n -y n | < ε, for any ε > 0. The Cauchy real x is smaller than y, written x < y, if there is an ε > 0 such that eventually x n + ε ≤ y n . There is no effective way to compare two Cauchy reals, so these numbers cannot be equipped with the structures described in Section 7.2. We resort to the construction of Cauchy reals provided by Cohen [START_REF] Cohen | Construction of real algebraic numbers in Coq[END_REF] [START_REF] Cohen | Formalized algebraic numbers: construction and first-order theory[END_REF]. His libraries provides a type for Cauchy reals, and implements field operations for these numbers. It also provides a tactic called bigenough, which eases construction of the effective moduli of convergence required in the proofs that a certain property on Cauchy reals is eventually true.

Computations

In this section, we consider some aspects of computing inside Coq which are relevant to the present part, both manual and automatic.

Using the unary representation of integers described in Section 7.1, the command:

Compute 100*1000.

which asks Coq to evaluate this product, triggers a stack overflow. For the purpose of running computations inside Coq's logic, on integers of a medium size, an alternate datastructure is required, together with less naive implementations of the arithmetical operations. The present formal proof requires this nature of computations at several places, for instance in order to evaluate sequences defined by a recurrence relation at a few particular values. For these computations, we used the binary representation of integers provided by the ZArith library included in the standard distribution of Coq, together with the fast reduction mechanism embarked inside Coq's kernel [START_REF] Grégoire | A compiled implementation of strong reduction[END_REF].

These two ingredients are also used behind the scene by tactics implementing decision procedures. For instance, the development makes extensive use of proof commands dedicated to the normalization of algebraic expressions like the field and ring tactics [START_REF] Mahboubi | Proving equalities in a commutative ring done right in Coq[END_REF]. The field tactic generates proof obligations describing sufficient conditions for the simplifications it made. For instance, the command ring might be used to automatically prove a goal of the shape (x -1) (x -2) = x 2 -3 x + 2, for x in a ring. On the goal (x -1) / (x -1) = 1, the command field would generate a side condition x <> 1 as a new goal, and prove in the background the implication that x <> 1 -> (x -1) / (x -1) = 1.

In our case, such side conditions in turn are automatically solved using the lia decision procedure for linear integer arithmetic [START_REF] Besson | Fast reflexive arithmetic tactics the linear case and beyond[END_REF] so that in practice, calling field and then lia completely solves the goal if these side conditions are implied by the current proof context.

In the case of automated tactics, the conversion of the formulas in the goal into data structures suitable for larger scale computation is not visible to the user: it is performed by extra-logical code which is part of the internal implementation of these tactics. The situation is different when a computational step in a proof requires the evaluation of a formula at a given argument, and when both the formula and the argument are described using proof-oriented, inefficient representations. In that case, for instance for evaluating terms in a given sequence, the CoqEAL library [START_REF] Cohen | Refinements for Free! In Certified Programs and Proofs[END_REF] was used, which provides an infrastructure automating the conversion between different datastructures and algorithms used to model the same mathematical objects, like different representations of integers or different implementations of a matrix product.

a n = n k=0 n k 2 n+k k 2 , b n = a n z n + n k=1 k m=1 (-1) m+1 n k 2 n+k k 2 2m 3 n m n+m m
.

Lemma 13 is the bottleneck in Apéry's proof; recall that it states that for n ≥ 0, (a n ) n∈N and (b n ) n∈N both satisfy the second-order recurrence:

(n + 2) 3 y n+2 -(17n 2 + 51n + 39)(2n + 3)y n+1 + (n + 1) 3 y n = 0.

Both sums a n and b n are instances of parametrised summation: they follow the pattern F n = β(n) k=α(n) f n,k in which the summand f n,k , potentially the bounds, and thus the sum, depend on a parameter n. This makes it appealing to resort to the algorithmic paradigm of creative telescoping, which was developed for this situation in computer algebra.

To demonstrate how this works, let us look at a very simple example. Suppose we want to find a recurrence satisfied by the sequence F n = n k=0 n k . We start by adopting a characterization of f n,k := n k through an implicit representation by two recurrences, which are natural if we think of the usual definition

n k = n! k!(n-k)! : n + 1 k = n + 1 n + 1 -k n k , n k + 1 = n -k k + 1 n k . (8.1)
We put aside for now the question of the initial conditions which would be necessary to fully characterize them. Creative telescoping derives the following formal relation, with finite difference with respect to k on the right-hand side:

n + 1 k -2 n k = n k + 1 - n + 1 k + 1 - n k - n + 1 k (8.2)
Let us check this formal identity, which is the output of an algorithm. If we expand the first part of the right-side of Equation (8.2) using the rules (8.1), we get

n k + 1 - n + 1 k + 1 - n k - n + 1 k = (8.3) 1 - n + 1 n -k n k + 1 - n k - n + 1 k = (8.4) - k + 1 n -k n -k k + 1 n k - n k - n + 1 k = (8.5) n + 1 k -2 n k . (8.6)
Notice that in order to prove this, we only used the equations (8.1) so that this property is general to any bivariate sequence satisfying them.

The idea is then to sum Equation (8.2) over k from 0 to n + 1, which gives:

F n+1 -2F n = n n + 2 - n + 1 n + 2 - (8.7) n 0 - n + 1 0 (8.8) =0 -(1 -1) = 0 (8.9)
since the right-hand side of Equation (8.2) was of the shape v k+1 -v k . From this we get that F n = 2 n F 0 , so that with the right initial conditions for f n,k , namely f n,0 = 1 and f n,n = 1, we get the usual identity n k=0 n k = 2 n . This kind of proof is a typical albeit very elementary example of what the book A=B [START_REF] Petkovšek | [END_REF] calls automating the discovery and proof of identities.

Given both the magical nature of Apéry's recurrence in Lemma 13 and the painful efforts of Van der Porten, Cohen and Zagier [START_REF] Van Der Poorten | A proof that Euler missed: Apéry's proof of the irrationality of ζ(3)[END_REF] to prove it, one could try to apply creative telescoping techniques to re-discover both the recurrence and its proof automatically. This was done by Salvy in a Maple sheet [START_REF] Salvy | An Algolib-aided version of Apéry's proof of the irrationality of ζ(3)[END_REF] using the Algolib [START_REF] Chyzak | Algolib[END_REF] library. From the formal proof point-of-view, the prospect of only checking what are essentially pre-baked proofs or rather, proof certificates, seems enticing enough.

Recurrences as a data structure for sequences

As already encountered in Equation (8.1), a fruitful idea from the realm of Computer Algebra is to represent sequences not explicitly, such as the univariate (n!) n or the bivariate ( n k ) n,k , but by a system of linear recurrences with coefficients which are polynomials in the variables on which this sequence depends, accompanied with sufficient initial conditions. We borrow the following description from Chyzak [START_REF] Chyzak | The ABC of Creative Telescoping-Algorithms, Bounds, Complexity[END_REF].

A sequence

(u n1,•••,nr ) n1,•••,nr∈N is called hypergeometric if for each i ∈ {1,• • • , r}, there exists a rational function R i (n 1 ,• • • , n r ) such that u n1,•••,ni+1,•••,nr = R i (n 1 ,• • • , n r ) • u n1,•••,nr ,
in other words, if u is cancelled by the operators R i -S ni where S ni is the shift operator which maps u n1,•••,nr to u n1,•••,ni+1,•••,nr . For example, n k is hypergeometric: the recurrences (8.1) once rewritten as equalities to zero can be represented as P • f = 0 for P = S n -n+1 n+1-k and P = S k -n-k k+1 , respectively. This notion can be generalized. First, consider a skew polynomial algebra

A = Q(n 1 ,• • • , n r ) < S n1 ,• • • , S nr > with the commutation rules S i1 n1 • • • S ir nr c(n 1 ,• • • , n r ) = c(n 1 + i 1 ,• • • , n r + i r )S i1 n1 • • • S ir nr for any rational fraction c(n 1 ,• • • , n r ). A polynomial P = (i1,•••,ir) p(n 1 ,• • • , n r )S i1 n1 • • • S ir nr ∈ A
acts on a r-variable sequence (u n1,•••,nr ) by

(P • u) n1,•••,nr = (i1,•••,ir) p(n 1 ,• • • , n r )u n1+i1,•••,nr+ir ,
where subscripts denote evaluation.

To any sequence f , one associates the set of skew polynomials that annihilate it. This set, {P ∈ A : P • f = 0} is a left ideal of A, named the annihilating ideal of f , and denoted ann f . A hypergeometric sequence can be redefined as a function f whose quotient module A/ ann f is a vector space of dimension 1 over the rational function field

Q(n 1 ,• • • , n r ). A ∂-finite sequence is a function whose quotient module A/ ann f has finite dimension ≥ 1.
A non-commutative extension of the usual Gröbner-basis theory is available, together with algorithmic analogues. In this setting, a good representation of a ∂-finite sequence is obtained as a Gröbner basis of ann f for a suitable ordering on the monomials in S n1 ,• • • , S nr . For the example of f n,k = n k , a Gröbner basis consists of both already-mentioned skew polynomials encoding (8.1). In general, a Gröbner basis provides us with a (vectorial) basis of the quotient module A/ ann f , which is isomorphic to Af . This basis can be explicitly written in the form B = {f n1+i1,•••,nr+ir } (i1,•••,ir)∈U , where the finite set U of indices is given as the part under the classical stair shape of the Gröbner-basis theory. Given a Gröbner basis GB for ann f , the normal form NF(p, GB) is unique for any p ∈ A. Again in the binomial example, the finite set is U = {(0, 0)}, and normal forms are rational functions. This is the basis of algorithms for a number of operations under which the ∂-finite class is stable, which all operate by looking for enough dependencies between normal forms: application of an operator, addition, product. Although the following applies in the general case, we now restrict ourselves to the bivariate case which interests us more specifically from now on.

The case of summing a sequence (f n,k ) into a parametrized sum F n = n k=0 f n,k is more involved: it follows the method of creative telescoping [START_REF] Zeilberger | The method of creative telescoping[END_REF], in two stages. First, an algorithmic step determines pairs (P, Q) satisfying

P • f = (S k -1)Q • f (8.10) with P ∈ Q(n)[S n
] and Q ∈ A. To continue with our example f n,k = n k , this corresponds to Equation 8.2 with P = S n -2 and Q = 1 -S n . Second, a systematic step follows: summing (8.10) for k between 0 and n + deg P yields

(P • F ) n = (Q • f ) k=n+deg P +1 -(Q • f ) k=0 .
(8.11) Continuing with our binomial example, this corresponds to the summation leading to Equation (8.7).

Apéry's sequences are ∂-finite constructions

The sequences a and b in (6.1) are ∂-finite: they have been announced to be solutions of (6.2). But more precisely, they can be viewed as constructed from "atomic" sequences by operations under which the class of ∂-finite sequences is stable. This is summarized in Table 8.1.

step explicit form GB operation input(s) 1 c n,k = n k 2 n+k k 2 C direct 2 a n = n k=1 c n,k A creative telescoping C 3 d n,m = (-1) m+1 2m 3 ( n m )( n+m m ) D direct 4 s n,k = k m=1 d n,m S creative telescoping D 5 z n = n m=1 1 m 3 Z direct 6 u n,k = z n + s n,k U addition Z and S 7 v n,k = c n,k u n,k V product C and U 8 b n = n k=1 v n,k B creative telescoping V Table 8
.1: Construction of a n and b n : At each step, the Gröbner basis named in column GB, which annihilates the sequence given in explicit form, is obtained by the corresponding on the inputs given in the last column.

In this table, Gröbner bases are systems of recurrence operators: at each line in the table, the sequence given in explicit form is a solution of the system of recurrences described by the operators in the Gröbner basis column. Note that in fact none of these results rely on the specific sequences in the explicit form: at each step, a new Gröbner basis is obtained from known ones, the ones that are cited in the input column, just like we obtained the operator S n+1 -2S n in Equation (8.7) from S n -n+1 n+1-k and S k -n-k k+1 . The table can also be read bottom-up for the purpose of verification: the Gröbner basis obtained at a given step can be verified using only the Gröbner bases obtained at some previous steps, all the way down to C and D. These operators describe a more general class of sequences than just the explicit sequences used in this table, thus initial conditions are needed to describe a precise sequence.

Provisos and sound creative telescoping

Let us go back to our proof of the binomial identity in Equations (8.3) to (8.6). This proof in fact only holds if n = k which makes taking the sum for k from 0 to n + 1 an illegitimate operation. In another example, one can "almost prove" Pascal's triangle rule using only the recurrences (8.1):

n + 1 k + 1 - n k + 1 - n k = n + 1 n -k n -k k + 1 - n -k k + 1 -1 n k = 0 × n k = 0,
but this requires k = -1 and k = n. Therefore, this does not prove Pascal's rule for all n and k. This incomplete modelling of sequences by algebraic objects may cast doubt on the output of creative-telescoping algorithms. By contrast, in our formal proofs, we augmented the recurrences with provisos that restrict their applicability. In this setting, we validate a candidate identity like the Pascal triangle rule by a normalization modulo the elements of a Gröbner basis plus a verification that this normalization only involves legal instances of the recurrences. In the case of creative telescoping, Equation (8.10) takes the form:

(n, k) / ∈ ∆ ⇒ (P • f ,k ) n = (Q • f ) n,k+1 -(Q • f ) n,k , (8.12)
where ∆ ⊂ Z 2 guards the relation and where f ,j denotes the univariate sequence obtained by specializing the second argument of f to j. Thus our formal analogue of Eq. (8.11) takes this restriction into account and has the shape

(P • F ) n = (Q • f ) n,n+β+1 -(Q • f ) n,α + r i=1 i j=1 p i (n) f n+i,n+β+j + α≤k≤n+β ∧ (n,k)∈∆ (P • f ,k ) n -(Q • f ) n,k+1 + (Q • f ) n,k , (8.13)
for F the sequence with general term F n = n+β k=α f n,k and where P = r i=0 p i (n)S i n , with the p i polynomials with rational coefficients.

The proof of identity (8.13) is a straightforward reordering of the terms of the left-hand side, (P •F ) n = r i=0 p(n) F n+i , after unfolding the definition of F and applying relation (8.12) everywhere allowed in the interval α ≤ k ≤ n + β. The first part of the right-hand side is the usual difference of border terms, already present in Equation (8.11). The middle part is a collection of terms that arise from the fact that the upper bound of the sum defining F n depends linearly on n and that we do not assume any nullity of the summand outside the summation domain. The last part, which we will call the singular part, witnesses the possible partial domain of validity of relation (8.12). The operator P is a valid recurrence for the sequence F if the right-hand side of Equation (8.13) normalizes to zero, at least outside of an algebraic locus ∆ that will guard the recurrence.

The Coq counterpart to Equation (8.13) is the following

Lemma sound_telescoping (cf0 : int -> R := fun _ => 0) (r : nat := size Pseq) (P : (int -> R) -> int -> R := horner_seqop Pseq) (U : int -> R := fun n => \sum_(a <= k < n + b :> int) u n k) (n : int) (range_correct : a <= n + b) (PeqDQ : ∀ (n k : int), not_D n k -> P (u ^~k) n = Q u n (k + 1) -Q u n k) : P U n = Q u n (n + b) -Q u n a + \sum_(a <= k < n + b :> int | ~~not_D n k) (P (u ^~k) n -(Q u n (k + 1) -Q u n k)) + \sum_(0 <= i < degree P) \sum_(0 <= j < i) P i * u (n + i) (n + j + b).
Listing 8.1: Our sound telescoping lemma in Coq

The notation u ^~k stands for fun n => u n k, P i denotes the i-th coefficient of the polynomial P. The outside of the locus ∆ is represented by the predicate not_D.

Let us now apply the theorem of Equation (8.13) to our example of the binomial identity. The forbidden domain ∆ is the singleton {(n, n)}, because that is the only point at which Equation (8.2) is not valid. We get

(S n -2)F n = F n+1 -2F n = ((f n,n+1 -f n+1,n+1 ) + (f n,0 -f n+1,0 )) + f n+1,n+1 + (((S n -2)f ,n ) n -(f n,n+1 -f n+1,n+1 ) + (f n,n -f n+1,n )) =(-1) + 1+ (n + 1 -2) -(0 -1) + (1 -(n + 1)) =0.
The automatic nature of the process of creative telescoping takes a hit here, both because ∆ has to be built manually and because it may depend on the order in which recurrences are rewritten in the proof of identity (8.12). Furthermore, we are aware of no guarantees that given a rewriting order, the various terms of Formula (8.13) will cancel as nicely as in the case of the binomial identity above.

Generated operators, hand-written provisos, and formal proofs

For each step in Table 8.1, we make use of the data computed by the Maple session in a systematic way. Figure 8.1 illustrates this pattern on the example of step 7. As mentioned in Section 8.1.3, we annotate each operator produced by the computer-algebra program with provisos (see below) and turn it this way into a conditional recurrence predicate on sequences. To each sequence in the program corresponds a file defining the corresponding conditional recurrences, for instance annotated_recs_c, annotated_recs_u, and annotated_recs_v for c, u, and v, respectively. More precisely these files contain all the operators obtained by the Maple script for a given sequence, not only the Gröbner basis. We use rounded boxes to depict the files that store the definitions of these predicates. These are generated by the Maple script which pretty-prints its output in Coq syntax, with the exception of the definition of provisos. Throughout this section, a maple leaf tags the files that are generated by our Maple script. Yet automating these annotations is currently out of reach.

In our formal proof, each step in Table 8.1 consists in proving that some conditional recurrences on a composed sequence can be proved from some conditional recurrences known for the arguments of the operation. We use square boxes to depict the files that store these formal proofs. The statement of the theorems proved in these files are composed from the predicates defined in the round boxes: a dashed line points to (predicates used to state) conclusions and a labelled solid line points to (predicates used to state) hypotheses. The definition of these recurrences is named to reflect the term it rewrites, e.g., the left-hand sides in (8.1): these names are the result of pretty-printing the (skew) monomial that encodes these left-hand sides, the prefix S standing for "shift". For example Sn is the name of a recurrence defining c n+1,k , while SnSk would be for c n+1,k+1 . Rewriting a given term with such an equation makes the term decrease for the order associated with the Gröbner basis. Another part of the file defines the recurrences obtained from a creative-telescoping pair (P, Q) generated for the purpose of the summation defining the sequence a.

Listing 8.2: Maple-generated part of annotated rec c

Observe that these generated definitions feature named provisos that are in fact placeholders. In the preamble of the file, displayed on Listing 8.3, we provide by a manual annotation a concrete definition for the proviso of each recurrence defined in the generated part. Observe however that part of these definitions can be inferred from the coefficients of the recurrences. For example the k = n + 1 condition in precond.Sn, the proviso of recurrence Sn, is due to the denominator (-n -1 + k) 2 of the coefficient (Sn_cf0_0 n k). The proof of the goal of Listing 8.6, proceeds in three steps. The first step is to inspect the terms in singular_part and to chase ill-formed denominators, like n -n. These can arise from the specialisations, like k = n, induced when unrolling the definition of (the negation of) precond.CT_premise. In our formalization, a division by zero is represented by a conventional value: we check that these terms vanish by natural compensations, independently of the convention, and we keep only the terms in singular_part that represent genuine rational numbers. The second step consists in using the annihilator ann_c of the summand to reduce the resulting expression under the stairs of the Gröbner basis. In fact, this latter expression features several collections of terms, that will be reduced to as many independent copies of the stairs. In the present example, we observe two such collections: (i) terms that are around the lower bound (n, 0) of the sum, of the form c n,0 , . . . c n,s ; (ii) terms that are around the upper bound (n, n) of the summation, of the form c n,n , . . . , c n,n+s for a constant s. The border terms induce two such collections but there might be more, depending in particular on the shape of the precond.CT_premise proviso. For example, the sum n k=0 (-1) k n k 3k n = (-3) n leads to a proviso involving n = 3k + 1 and similar terms: an additional category of terms around (n, n/3) drifts away from both (n, 0) and (n, n) when n grows. The collection of terms around the upper bound in our running example is displayed on Listing 8.7. The script of Listing 8.8 reduces this collection under the stairs of ann_c, producing the expression displayed on Listing 8.9. The premise of each rule in this basis being an integer linear arithmetic expression, we check its satisfiability using our front-end intlia to the lia proof command [START_REF] Besson | Fast reflexive arithmetic tactics the linear case and beyond[END_REF], which automates the formal proof of first-order formulae of linear arithmetics. The third and last step consists in checking that the rational-function coefficient of every remaining evaluation of c is zero. For this purpose, we start by unfolding the definitions of the coefficients P_cf2, Sn_cf0. Previous steps kept them carefully folded as these values play no role in the previous normalizations but can lead to large expressions if expanded, significantly slowing down any other proof command. The resulting rational function is proved to be zero by a combination of the field [START_REF] Delahaye | Dealing with algebraic expressions over a field in Coq using Maple[END_REF] and lia [START_REF] Besson | Fast reflexive arithmetic tactics the linear case and beyond[END_REF] proof commands. The former reduces the rational equation into a polynomial one between the cross product of two rational functions. This equation is then solved by the ring proof command [START_REF] Mahboubi | Proving equalities in a commutative ring done right in Coq[END_REF]. The algebraic manipulations performed by field produce a set of non-nullity conditions for the denominators. These are solved by the lia proof command. To this end, our Maple script generates rational fractions with factored denominators, that happen to feature only linear factors in these examples. formally in file algo_closures that each sequence involved in the construction of a n and b n is a solution of the corresponding Gröbner system of annotated recurrence, starting from c n , d n , and z n and applying the lemmas proved in the ops_for_* files all the way to the the final conclusions of ops_for_a and ops_for_b. This proves that a n is a solution of the recurrence (6.2) but provides only a recurrence of order four for b n ; let us call this recurrence R 4 .

Composing closures and reducing the order of B

In file reduce_order, we prove that b as well satisfies the recurrence (6.

Consequences of Apéry's recurrence

In this section, we detail the elementary proofs of the properties obtained as corollaries of Lemma 13. We recall from Chapter 6 that these properties describe the asymptotic behavior of the sequence δ n = a n ζ(3) -b n , with:

a n = n k=0 n k 2 n+k k 2 , b n = a n z n + n k=1 k m=1 (-1) m+1 n k 2 n+k k 2 2m 3 n m n+m m (8.14)
In this section, we use the vocabulary and notations of Cauchy real numbers (see Section 7.3). For instance, we have the equality between the two Cauchy reals:

Lemma 15. ζ(3) = b a
Proof. Easy from the definition of the sequences a and b.

However, the Mathematical Components libraries do not cover any topic of analysis, and even the most basic definitions of transcendental functions like the exponential or the logarithm are not available. Fortunately, it is possible to obtain the required properties of the sequence δ by very elementary means, and almost all these elementary proofs can be inferred from a careful reading and a combination of Salvy's proof [START_REF] Salvy | An Algolib-aided version of Apéry's proof of the irrationality of ζ(3)[END_REF] and of van der Poorten's description [START_REF] Van Der Poorten | A proof that Euler missed: Apéry's proof of the irrationality of ζ(3)[END_REF].

Following van der Poorten, we introduce an auxiliary sequence (w n ) ∈ Q n , defined as:

w n = b n+1 a n+1 b n a n = b n+1 a n -a n+1 b n .
The sequence w is called a Casoratian: as a and b are solutions of a same linear recurrence relation (6.2) of order 2, this can be seen as a discrete analogue of the Wronskian for linear differential systems. For example, w satisfies a recurrence relation of order 1, which provides a closed form for w:

Lemma 16. For n ≥ 2, w n = 6 (n+1) 3 .
Proof. Since a and b satisfy the recurrence relation (6.2), we have for n ∈ N:

b n+2 a n+2 b n+1 a n+1 = (17n 2 +51n+39)(2n+3) (n+2) 2 -(n+1) 3 (n+2) 3 1 0 • b n+1 a n+1 b n a n .
Taking determinants on both sides, we get

w n+1 = (n + 1) 3 (n + 2) 3 w n .
The result then follows from the computation of w 0 .

From this formula, we can obtain the positivity of the sequence δ, and an evaluation of its asymptotic behavior in terms of the sequence a. ).

Proof. Since ζ(3) = b a , it is sufficient to show that there exists a constant K, such that for any

k < l, b l a l -b k a k ≤ K a 2 k
. But since a is an increasing sequence, Lemma 16

proves that for any k < l, b l a l -b k a k ≤ l-1 i=k wi aiai+1 ≤ l-1 i=k 6 (i+1) 3 a 2 k ≤ K a 2 k
, for any K greater than 6ζ(3).

The last remaining step of the proof is to show that the sequence a grows fast enough.

Lemma 17. 33 n = O(a n ).
Proof. Consider the auxiliary sequence ρ n = an+1 an . Since ρ 51 is greater than 33, we only need to show that the sequence ρ is increasing. For the sake of readability, we denote µ n and ν n the fractions coefficients of the recurrence satisfied by a, obtained from Equation (6.2) after division by its leading coefficient. Thus a satisfies the recurrence relation:

a n+2 -µ n a n+1 + ν n = 0.
For n ∈ N, we also introduce the function h n (x) = µ n + νn x , so that ρ n+1 = h n (ρ n ). The polynomial P n (x) = x 2 -µ n x + ν n has two distinct roots x n < x n , and the formula describing the roots of polynomials of degree 2 show that 0 < x n < 1 < x n and that the sequence x n is increasing. But since h n (x) -x = -Pn(x)

x , for 1 < x < x n , we have h n (x) > x. A direct recurrence shows that for any n ≥ 2, ρ n ∈ [1, x n ], which concludes the proof.

In the formal proof of Lemma 17, the computation of ρ 51 was made possible by using the CoqEAL library, with the help of Cyril Cohen. This proof also requires a few symbolic computations that are a bit tedious to perform by hand: in these cases, we used Maple as an oracle to massage algebraic expressions, before formally proving the correctness of the simplification. This was especially useful to study the roots x n and x n of P n .

We can now conclude with the limit of the sequence 3 n δ n , under the assumption that n = O(3 n ).

Corollary 3. lim n→∞ ( 3 n δ n ) = 0.
Proof. Immediate, since δ n = O( 1 an ) by Corollary 2, and 3 n = O((3 3 ) n ), and 3 3 < 33.

In the next Section, we describe the proof of the last remaining assumption, about the asymptotic behavior of n .

Asymptotics of lcm(1, ..., n)

The asymptotic behavior of the sequence lcm(1, ..., n) is an elementary corollary of the Prime Number Theorem. The proof of this corollary is based on a simple remark about the p-valuations of lcm(1, ..., n).

Remark 2. For any integer n, let n denote the least common multiple of the integers no greater than n, denoted lcm(1, ..., n). For any prime number p, p vp( n ) is the highest power of p not exceeding n, so that:

v p ( n ) = log p (n) .
Proof. Noticing that v p (lcm(a, b)) = max(v p (a), v p (b)), we see by induction on n that v p ( n ) = n max i=1 v p (i). Recall from Section 7.1 that log p (n) is a notation for the greatest integer α such that p α ≤ n. Since α = v p (p α ), we have α ≤ v p ( n ). Now suppose that v p ( n ) = v p (i) for some i ∈ {1,• • • , n}. Then i = p vp(i) q with gcd(p, q) = 1 so that p vp( n) = p vp(i) ≤ i ≤ n and thus v p ( n ) ≤ α. This proves that v p ( n ) = α.

By Remark 2, n can hence be written as p≤n p log p (n) and therefore:

ln( n ) = p≤n log p (n) ln(p).
If π(n) is the number of prime numbers no greater than n, we hence have:

ln( n ) ∼ π(n) ln(n).
The Prime Number theorem states that π(n) ∼ n ln(n) ; we can thus conclude:

n = e n (1+o(1)) .

J. Avigad and his co-authors provided the first machine-checked proof of the Prime Number theorem [START_REF] Avigad | A formally verified proof of the prime number theorem[END_REF], which was considered at the time as a formalization tour de force. Their formalization is based on a proof by A. Selberg and P. Erdös. Although the standard proofs of this theorem use tools from complex analysis like contour integrals, their choice was guided by the corpus of formalized mathematics available for the Isabelle proof assistant, or the limits thereof. Although less direct, the proof by A. Selberg and P. Erdös is indeed more elementary and avoids complex analysis completely.

Notations and Outline

In order to prove Corollary 3 in Section 8.2, we already mentioned in Chapter 6 that a weaker description of the asymptotic behavior of ( n ) is sufficient, stated in Lemma 11:

n = O(3 n ). (8.15)
This sort of weaker description was in fact known before the first proofs of the Prime Number theorem but our formal proof is a variation on an elementary proof proposed by Hanson [START_REF] Hanson | On the product of the primes[END_REF].

The idea of the proof is to replace the study of n by that of another sequence C(n). The latter is defined as a multinomial coefficient depending on elements of a fast-growing sequence α. The fact that

n i=1 α 1/αi i < 3 independently of n then allows to show that C(n) = O(3 n ).

Proof

Define the sequence (α n ) n∈N by α 1 = 2, and

α n+1 = α 1 α 2 • • • α n + 1 for n ≥ 1.
By an induction on n, this is equivalent to

α n+1 = α 2 n -α n + 1. For n, k ∈ N, let C(n, k) = n! n/α 1 ! n/α 2 !• • • n/α k ! .
As soon as α k ≥ n, C(n, k) is independent of k and we denote C(n) = C(n, k) for all such k. Hanson directly works with C(n), but we found this to be inconvenient to manipulate in the proof. Moreover, most inequalities stated on C(n) actually hold for C(n, k) with little or no more hypotheses. Notice that

C(n, k) = n n/α1 , n/α2 ,•••, n/α k • n - k i=0 n αi ! . In particular, C(n, k) ∈ N.
The goal is now to show that n ≤ C(n) < K • 3 n for some K.

Lemma 18. For k ∈ N, k i=1 1 α i = α k+1 -2 α k+1 -1 < 1 and thus for x ∈ Q with x ≥ 1, x > k i=1 x α i .
The proof is done by induction and relies on the fact that if a ∈ Q and m ∈ N + , we have a m = a m . In the following, for n, k ∈ N and p prime, we denote β p (n, k) for val p (C(n, k)).

Lemma 19. For all n, k ∈ N and p prime,

β p (n, k) ≥ log p (n) = val p ( n ). Therefore C(n, k) ≥ n .
Proof. The proof uses Lemma 14.

β p (n, k) =val p (n!) - k i=1 val p ( n/α i !) = log p (n) i=1 n/p i - k i=1 log p ( n α i ) j=1 n α i p j = log p (n) i=1   n/p i - k j=1 n p i α j   ≥ log p (n) i=1 1 (thanks to the fact that 1 α i < 1 from lemma 18). Since n = p≤n p log p (n) from remark 2, we get n ≤ C(n, k) = p≤n p βp(n,k) . Lemma 20. For i ≥ 1 and n ≥ α i , n αi n α i n αi n α i < 10 n α i α i -1 α i . Proof. If n = α i , we have 1 < √ 10 ≤ 10 α i -1
α i , hence the result. Otherwise n > α i and then, if we write n = bα i + r, with 0 ≤ r < α i , we have n-αi+1 αi

-n αi = -1 + 1 αi < 1 2 < 0 so that n-αi+1 αi < n αi .
It is easy to deduce that n-αi+1 αi

n-α i +1 α i < n αi n α i , hence we have n αi n α i n αi n α i < n αi n α i n-αi+1 αi n-α i +1 α i = 1 + α i -1 n -α i + 1 n-α i +1 α i -1 α i -1 α i • n α i α i -1 α i
.

The first operand in the last expression is of the shape 1 + 1

x x α i -1 α i
, where x is rational. We showed using only elementary properties of algebraic numbers that for x ∈ Q, 0 < x, 1 + 1 x x < 10, hence the result. Note that we only needed that there exist a constant K > 0 such that 1 + 1

x x < K.
Of course, using elementary real analysis allows for the tighter bound e for 1 + 1 x x , which was used in Hanson's paper, but this bound is irrelevant for the final result. Using the previous facts, we obtain the following bound:

Lemma 21. For k ≥ 1 and n ≥ 2, C(n, k) < n n n α1 n α 1 • • • n α k n α k . Proof. First observe that if m = m 1 +• • • + m k
where m and the m i are (not all zero) nonnegative integers, we have because of Theorem 1:

(m 1 +• • • + m k ) m ≥ m m 1 ,• • • , m k m m1 1 • • • m m k k . (8.16)
Let k ≥ 1, and define:

t = k i=1 n α i .
Then t < n by lemma 18. We have

C(n, k) =n • (n -1)• • • (t + 1) t n/α 1 , n/α 2 ,• • • , n/α k . (8.17)
Because of equation ( 8.16), we know that

t n/α 1 , n/α 2 ,• • • , n/α k ≤ t t n/α 1 n/α1 n/α 2 n/α2 • • • n/α k n/α k . (8.18)
From equations (8.17) and (8.18) we deduce that

C(n, k) < n n n/α 1 n/α1 n/α 2 n/α2 • • • n/α k n/α k . Lemma 22. Let k ≥ 3, n ∈ N. If α k ≤ n then k < log 2 log 2 n + 2.
Proof. First observe by a simple induction that for all k ≥ 3,

α k > 2 2 k-2 + 1 so that k -2 < log 2 ( log 2 α k ) ≤ log 2 ( log 2 n ) . Lemma 23. Let k ≥ 1, n ∈ N. If α k ≤ n, C(n, k) < n n ( 10 n α1 ) α 1 -1 α 1 ( 10 n α2 ) α 2 -1 α 2 • • • ( 10 n α k ) α k -1 α k
10 n α1

10 n α 1

10 n α1

10 n α 1 • • • 10 n α k 10 n α k .
Proof. The result is straightforward by combining lemmas 20 and 21.

Lemma 24. Let w k = k i=1 α 1 α i i , k ≥ 1.
Then w k is increasing and there exists w ∈ R, with w < 2.98, such that w k < w.

Proof. The sequence w k is increasing because α i

1 α i > 1 (because α i > 1). Since α 2 i > α i+1 > (α i -1) 2 , one can see that for i ≥ 3, α 1 α i+1 i+1 < α 1 α i
i , so that for all k ≥ 1 and l ≥ 0,

w k+l ≤ k i=1 α 1 α i i • α 1 α k+1 l i=0 1 2 i k+1 ≤ w k • α 2 α k+1 k+1 .
We establish by an elementary external computation verified in Coq that

α 1 -1 α 1 + α 2 -1 α 2 +• • • + α k -1 α k = k -1 + 1 α k+1 -1 .
Proof. It is a direct consequence of Lemma 18.

Note that the statement of Remark 3 actually corrects a typo in the original paper.

Theorem 6. If α k ≤ n < α k+1 , C(n, k) = C(n) < (10 n) k-1 2 w n+1 .
Proof. From lemma 23, recall that we have

C(n, k) < n n ( 10 n α1 ) α 1 -1 α 1 ( 10 n α2 ) α 2 -1 α 2 • • • ( 10 n α k ) α k -1 α k n α1 n α 1 n α2 n α 2 • • • n α k n α k = n n (10 n) k i=1 α i -1 α i k i=1 α 1 α i i n n n k i=1 1 α i k i=1 α α i -1 α i i .
It can be seen using lemma 18 that:

n n 1-k i=1 1 α i ≤ n. Thus C(n, k) <n (10 n) k-1+ 1 α k+1 -1 w n k k i=1 α α i -1 α i i thanks to remark 3 ≤n (10 n) k-1+ 1 α k+1 -1 w n k i=1 α α i -1 α i i because w k ≤ w. Since n < α k+1 = 1 + k i=1 α i , n ≤ k i=1
α i , and we have

k i=1 α α i -1 α i i = k i=1 α i w k ≥ n w k . Thus C(n, k) <(10 n) k-1+ 1 α k+1 -1 w n w k ≤(10 n) k-1 2 w n+1 as α k+1 ≥ 3 and w k ≤ w.
We can now prove lemma 11:

Proof. n /3 n ≤ C(n, k)/3 n = (10 n) k-1 2 w 3 n+1 .
Remembering that k < log 2 log 2 n + 2 and w < 3, it is elementary to show that the quantity on the right is eventually decreasing to 0 and therefore bounded, which proves the result. We once again make use of the fact that 1 + 1 x x is bounded in the course of this elementary proof.

Chapter 9

Conclusion 9.1 Proving Inequalities in Coq

In the proof of Section [START_REF] Hanson | On the product of the primes[END_REF], we repeatedly had to prove inequalities between expressions in an ordered ring or field. This was quite tedious for two reasons.

We present these reasons and suggest ideas for tools which could relieve the user. Such tools would replace something humans do without thinking about it.

Chained Inequalities: Proofs of statements of the style e 1 ≤ e n which are done using a chain of inequalities e 1 ≤ e 2 ≤ • • • ≤ e n are currently wearisome, especially in the case where the e i are big expressions. The way one typically proceeds is by first proving that e 1 ≤ e 2 , which involves writing down the expressions of e 1 and e 2 in Coq. Once this is done, one invokes transitivity of inequality so that one is left with the goal e 2 ≤ e n . One then states e 2 ≤ e 3 , which again involves writing down both e 2 and e 3 , and so on. This incurs a lot of manual copy-pasting, which makes the code inelegant and the work of doing proofs tedious. One could hope for something similar as the "calc" environment in the Lean theorem prover [START_REF] De Moura | The Lean Theorem Prover (System Description)[END_REF]. For example, we could have a tool allowing something in the style where the vernacular # <= X on a goal of the shape a ≤ b would apply transitivity and create the two new goals a ≤ X and X ≤ b.

Side Conditions: We lack a tool for the automation of side-condition proving in inequalities. Suppose we want to prove inequalities in an ordered ring R between products

x 1 • y 1 • z 1 ≤ a 1 • b 1 • c 1 ,
where x 1 , z 1 , a 1 and c 1 are a (possibly empty) product of nonnegative expressions, and y 1 and b 1 are nonnegative elements of R.

Typically, this is done by inlining inequalities on y 1 and b 1 alone. For example, we may know that y 1 ≤ y 2 , and we can use this fact so that the initial goal

x 1 • y 1 • z 1 ≤ a 1 • b 1 • c 1 becomes x 1 • y 2 • z 1 ≤ a 1 • b 1 • c 1 .
Thus, we do a series of 'rewritings' of the shape

x 1 • y 1 • z 1 ≤ a 1 • b 1 • c 1 x 1 • y 2 • z 1 ≤ a 1 • b 1 • c 1 ,
justified by proving that y 1 ≤ y 2 . To do so, the only lemma at our disposal is typically of the shape:

∀a, b, c, d, 0 ≤ a ⇒ 0 ≤ b ⇒ a ≤ b ⇒ c ≤ d ⇒ a • c ≤ b • d.
This means that we will have to re-prove several times that each quantity appearing in the product is nonnegative in the course of the rewritings. Moreover, this lemma only works directly for products of two elements of R; when we apply it to the product of k elements, we need to apply it k -1 times. .

A tool to tackle this problem would have to be able to deal with products of arbitrary arity, to propose a sensible deterministic mechanism to describe where and how we want to rewrite inequalities, and to keep track of already proven goals (typically the non-negativity of terms in the product).

Casts

As mentioned in Section 7.2, inclusions between sets in the Mathematical Components hierarchy are represented by injective functions between datatypes. Sometimes, coercions and notations can make two different representations indistinguishable to the eye in a Coq session. One consequence is that the user will fail to be able to rewrite or apply a theorem on what seems to be a valid instance. For example, if n is of type nat, Posz (n + 1) and (Posz n) + 1 can both appear as n+1, because the coercion Posz is hidden by default. Rewriting n + 1 to S n (successor of n) will succeed in the first case and fail in the second.

For a similar reason, proof commands like ring will fail on what seems like a valid equation, and only careful inspection of the Coq terms in the goal will reveal the source of failure. In fact, we had to rewrite tactics which started by normalizing expressions to a standard form in order to be able to apply the ring and field proof commands to the type rat representing the set Q.

However, the formalization of nontrivial mathematics would not be possible if it were not for the abundant notations and coercions at our disposal to be able to "forget" the precise nature of every object manipulated, just like we do when we write proofs on paper. A balance must be struck between an abundance of tools which may lead to not understanding what object we are manipulating, or which may lead to the inexplicable failure of proof automation tools, and a lack of such tools which would prevent formal proofs from being done.

Related Work

Number theory is an area of mathematics where very simple statements often require extremely difficult proofs. We are not aware of a comprehensive formal proof library on the topic, although there exist formal proofs of a few emblematic results. The elementary fact that √ 2 is irrational was used as an example problem in a comparative study of the styles of various theorem provers [START_REF] Wiedijk | The Seventeen Provers of the World: Foreword by[END_REF],

including Coq. The Prime Number theorem was proved formally for the first time by Avigad et al. [START_REF] Avigad | A formally verified proof of the prime number theorem[END_REF] and later by Harrison [START_REF] Harrison | Formalizing an analytic proof of the Prime Number Theorem[END_REF]. The transcendence of e was formalized by Bingham [START_REF] Bingham | Formalizing a proof that e is transcendental[END_REF] in HOL-Light. More recently, the transcendence of both π and e was formalized in Coq [START_REF] Bernard | Formal proofs of transcendence for e and π as an application of multivariate and symmetric polynomials[END_REF]. Some of the ingredients needed in the present proof are however not specific to number theory. For instance, we here use a very basic infrastructure to represent asymptotic behaviors, but "big Oh" notations have been discussed by Avigad [START_REF] Avigad | A formally verified proof of the prime number theorem[END_REF] in his formalization of the Prime Number Theorem, and by Boldo et al. [START_REF] Boldo | Formal Proof of a Wave Equation Resolution Scheme: the Method Error[END_REF] in a continuous context. Another example of such a secondary topic is the theory of multinomial coefficients, which is also relevant to combinatorics, and which is also defined by Hivert in his Coq library Coq-Combi [START_REF] Hivert | Coq-combi[END_REF]. However, up to our knowledge this library does not feature a proof of the generalized Newton identity.

We used Cohen's Cauchy reals to define ζ(3), and algebraic complex numbers in our simplification of Hanson's proof. We considered using instead the C-CoRN library of constructive analysis [START_REF] Cruz-Filipe | C-CoRN, the constructive Coq repository at Nijmegen[END_REF], but we found out that it would have been too much work compared to what we eventually did, in part because of the necessity to write a "bridge" to link mathematical structures from the Mathematical Components libraries to those of C-CoRN.

Harrison [START_REF] Harrison | Formal proofs of hypergeometric sums (dedicated to the memory of andrzej trybulec)[END_REF] recently presented a way to produce rigorous proofs from certificates produced by the Wilf-Zeilberger certificates, by seeing sequences as limits of complex functions. His method applies to the sequence a, which satisfies the recurrence Equation (6.2). However, this method does not allow for a proof that b satisfies Equation (6.2), because the summand is itself a sum but is not hypergeometric. Up to our knowledge, there is no way to justify the output of the efficient algorithms of creative telescoping used here without handling a trace of provisos.

As already mentioned, the idea to use computer algebra software (CAS) as an oracle outputting a certificate to be checked by a theorem prover, dubbed a skeptic's approach, was described in several papers [8] [7] [65]. This technique takes the best of both worlds to produce reliable proofs requiring large scale computations. In the case of Coq, this viewpoint is especially fruitful since the kernel of the proof assistant includes efficient evaluation mechanisms for the functional programs written inside the logic [START_REF] Grégoire | A compiled implementation of strong reduction[END_REF]. Notable successes based on this idea include the use of Pocklington certificates to check primality inside Coq [START_REF] Grégoire | A Computational Approach to Pocklington Certificates in Type Theory[END_REF] or external computations of commutative Groebner bases, with applications for instance in geometry [START_REF] Pottier | Connecting Gröbner bases programs with Coq to do proofs in algebra, geometry and arithmetics[END_REF]. Delahaye and Mayero proposed [START_REF] Delahaye | Dealing with algebraic expressions over a field in Coq using Maple[END_REF] to use CAS to help experimenting with algebraic expressions inside a proof assistant, before deciding what to prove and how to prove it. Unfortunately, their tool was not usable in our case, where algebraic expressions are made with operations that come from a hierarchy of structures.

Organizing the cooperation of a CAS and a proof assistant sheds light on their respective differences and drawbacks. The initial motivation of this work was to study the algorithms used for the automatic discovery and proof of recurrences. Our hope was to be able to craft an automated tool providing formal proofs of recurrences, by using the output of these algorithms, in a skeptical way. This plan did not work and Section 8.1 illustrates the impact of confusing the rational fractions manipulated by symbolic computations with their evaluations, which should be guarded by conditions on the denominators. On the other hand, proof assistants need to be handled carefully if one wants to manipulate the large expressions imported during the cooperation, even those which are of a small to moderate size for the standards of computer algebra systems. For instance, we have highlighted in our previous report [START_REF] Chyzak | A Computer-Algebra-Based Formal Proof of the Irrationality of ζ(3)[END_REF] the necessity to combine two distinct natures of data-structures in our libraries: one devoted to formal proofs, which may use computation inside the logic to ease bureaucratic steps in formal proofs, and one devoted to large scale computations, which provides a fine-grained control on the complexity of operations.

On several occasions in this work, we wrote more elementary versions of the proofs than what we had found in the texts we were formalizing. We agree with Avigad [START_REF] Avigad | A formally verified proof of the prime number theorem[END_REF] when he says that this can be both frustrating and enjoyable: on one hand, it can illustrate the lack of mathematical libraries for theorems which mathematicians would find simple, such as elementary analysis for studying the asymptotics of sequences as in Section 8.2. Ten years later, "the need for elementary workarounds" is still present, despite his fear that it would "gradually fade, and with it, alas, one good reason for investing time in such exercises" [START_REF] Avigad | A formally verified proof of the prime number theorem[END_REF]. On the other hand, this need gives an opportunity to better understand the minimal scope of mathematical theories used in a proof, with the help of a computer. For instance, it was not clear to us that we could manage to completely avoid the need to define transcendental functions, or even the constant e, to formalize Hanson's paper [START_REF] Hanson | On the product of the primes[END_REF], although it sometimes comes at the price of some slightly pedestrian calculations.

Chapter 10 Introduction

Sometimes validating external computations is not an option. For some problems, there does not seem to exist a certificate for which verification is cheaper than computation. In those cases, what matters is that there exist a formal specification of the way the result is obtained, so that one can prove that this result is correct by construction.

The problem described in this part of the thesis falls in the latter category. The goal is to estimate integrals numerically: given two bounds u and v and a function f : R → R, find A and B such that

A ≤ v u f (t)dt ≤ B (10.1)
and with the interval [A; B] preferably tight. This problem was already posed long before computers, for example when trying to estimate π by enclosing the unit disk between two n-gons. Numerical integration distinguishes itself from symbolic integration which aims at finding an equivalent exact mathematical formula for an integral, for example:

1 0 cos(t)dt = sin(1).
However, this is not always possible using elementary functions, as in the example, for x ∈ R, of

x 0 e t 2 dt.
Harnessing the power of computers to obtain a precise numerical value for previously unattainable integrals is becoming an issue for mathematicians, as estimates of numerous values of integrals become part of some modern mathematical proofs, such as that of the Ternary Goldbach Conjecture by Helfgott [START_REF] Helfgott | Major arcs for Goldbach's problem[END_REF] or the Double Bubble Conjecture by Hass and Schlafly [START_REF] Hass | Double bubbles minimize[END_REF]. The problem of numerical integration is typically presented as an approximation problem. The goal is to find a natural number n, as well as n points

u ≤ x 1 ≤ x 2 ≤ • • • ≤ x n ≤ v and n weights w 1 , w 2 ,• • • , w n ∈ R such
that that the integral of f is well approximated by a linear combination of evaluations of f :

v u f w 1 f (x 1 ) + w 2 f (x 2 ) +• • • + w n f (x n ).
(10.2)

Many numerical integration methods exist with different choices of weights. We borrow from Fousse [START_REF] Fousse | Intégration numérique avec erreur bornée en précision arbitraire[END_REF] the following list: the family of Newton-Cotes methods, the Romberg method, the Gaussian quadrature methods, the Gauss-Kronrod quadrature, the Clenshaw-Curtis method and the double-exponential method. With enough hypotheses on f and its derivatives when they exist, each method comes with an error bound for the expression on the right of Equation (10.2).

As an example, the 2-point trapezoidal rule, which is an instance of the Newton-Cotes method, uses n = 2, x 1 = u, x 2 = v, and

w 1 = w 2 = v-u 2 . Then if f ∈ C 2 ([u; v]), v u f (t) - v -u 2 (f (x 1 ) + f (x 2 )) = - (v -u) 3 f (ξ) 12 (10.3) for some ξ ∈ [u; v] 1 .
The first problem with this approach is that this error bound analysis is done assuming that all the performed operations (e.g. additions, multiplications, divisions) are exact. Most of the actual implementations of quadrature methods are made using floating-point number operations, and are thus prone to roundoff errors.

Another problem is that an error bound such as in Equation ( 10.3) is often not used systematically by mathematical software to provide an error interval at the end of the computation, but simply as a heuristic to both justify that the method does eventually converge to the value of the integral, and to know how to adjust the parameters to reach the user-requested error.

There exists rigorous software which takes into account both these pitfalls. By "rigorous", we mean that the correctness of the method used by the software is proved on paper. The INTLAB [START_REF] Rump | Verification methods: Rigorous results using floatingpoint arithmetic[END_REF] library by Rump uses a rigorous Romberg method [108], which is a Newton-Cotes quadrature formula, to integrate four times differentiable functions using interval arithmetic for both the value and the error bound. The software library Correctly Rounded Quadrature (CRQ) by Fousse [START_REF] Fousse | Intégration numérique avec erreur bornée en précision arbitraire[END_REF] implements integration algorithms with rigorous error bounds for the Newton-Cotes and Gauss-Legendre methods.

At the end of his PhD thesis, Fousse calls for the validation of his algorithms using formal proofs. To our knowledge, formally verified computations currently provide the strongest assurance of correctness. This is appropriate in order to obtain the kind of guarantee required by large mathematical proofs. The main contribution of this part is the description of a tool to automatically verify integral enclosures inside the Coq proof assistant.

In Section 11, we start by giving an introduction to the efficient and rigorous numeric real computation methods which we use to obtain our enclosures, along with the Coq libraries we rely on.

In Section 12, we describe how we obtain formally proved enclosures of both proper and improper integrals inside the Coq proof assistant. We present benchmarks on various integrals found in the literature. These benchmarks include two instances where we found previously undetected mistakes.

Finally, in Section 13, we compare these results with a custom implementation of Simpson's method; on the way, we present several iterations of our attempts to implement automatic differentiation of order n, for fun and profit.

Notations and definitions

In the following, an interval is a closed connected subset of the set of real numbers. We use I to denote the set of intervals:

I := {[a; b] | a, b ∈ R ∪ {±∞}}.
A point interval is an interval of the shape [a; a] where a ∈ R. The set I can be equipped with an addition

+ I :I × I → I ([a; b], [c; d]) → [a + c; b + d]
with neutral element [0; 0] which makes I a commutative monoid. The set I can also be equipped with a multiplication

• I :I × I → I ([a; b], [c; d]) → [min(a • c, a • d, b • c, b • d); max(a • c, a • d, b • c, b • d)]
with neutral element [1; 1]. When there will be no ambiguity, + I (resp. • I ) will be simply noted + (resp. •).

Any interval variable will be denoted using a bold font. For any interval x ∈ I, inf x (resp. sup x) denotes its left (resp. right) bound, with inf

x ∈ R ∪ {-∞} (resp. sup x ∈ R ∪ {+∞}). An enclosure of x ∈ R is an interval x ∈ I such that x ∈ x.
Definition 11. The closed convex hull of a set A ⊆ R is the smallest interval containing A, denoted here hull(A). Moreover, the interval hull(a, b) denotes the convex hull of (the union of ) two intervals a and b. Finally, hull(a, +∞) designates the interval [inf a; +∞).

Throughout the present part, we manipulate real numbers which are the evaluation of expressions. These expressions belong to a set E built from π, n ∈ Z, a countable set of variables x, y, z, • • •, the addition (+), the product (•), division (/), the inverse (• -1 ), the absolute value (|•|), the integer power (• n ), the exponential function (exp), the natural logarithm function (ln), and the trigonometric functions cos, sin, tan, and arctan.

Fix an enumeration of the variables. We define a function eval R in the following way. Let s ∈ N. The function eval R takes as input an expression e ∈ E, a vector x = (x 1 ,• • • , x s ) ∈ R s , and it returns the evaluation eval R (e, x) of the expression by following inductively the structure of e, applying the corresponding real operation, and replacing the i-th variable with x i . Whenever an operation is not defined or a variable of index greater than s appears, the function eval R returns an error ⊥.

When there is no ambiguity, we identify a variable-free expression with its evaluation. When we write "e(x) ∈ E", this means that e(x) is an expression whose only variable is x. If there is no ambiguity, we write t → e(t) for the function t → eval R (e, t), and if t ∈ R, we write e(t) for eval R (e, t).

Interval Arithmetic

The paradigm of interval arithmetic As the name indicates, interval arithmetic is a paradigm of computation in which, instead of manipulating real values, one handles intervals of real values. For example, we might represent the number π by the interval: One of the goals of interval arithmetic is, given a real number x ∈ R presented as an expression e ∈ E, to build an interval x such that x ∈ x.

(11.1)

Naive interval arithmetic

There are several ways to implement interval arithmetic. In the simplest one called naive interval arithmetic, usual unary and binary operations are overloaded to the realm of intervals and an enclosure of an expression e is built by composing these overloaded operators following the inductive structure of e. The following inclusion properties ensure that equation (11.1) holds by an induction on the structure of e ∈ E.

• Given a unary operator on real numbers, naive interval arithmetic provides a unary operator ♦ on intervals such that

∀x ∈ R, ∀x ∈ I, x ∈ x ⇒ (x) ∈ ♦(x). (11.2)
♦ is said to be an interval extension of .

An example would be the exponential function: if = exp, define

Exp : I → I [a, b] → [exp(a), exp(b)].
Since exp is an increasing function, ♦ := Exp satisfies (11.2).

• Given a binary operator on real numbers, naive interval arithmetic provides a binary operator ♦ on intervals such that

∀x, y ∈ R, ∀x, y ∈ I, x ∈ x ∧ y ∈ y ⇒ x y ∈ x♦y. (11.3)
In this binary case, we also say that ♦ is an interval extension of .

If we take to be the multiplication operation, then • I defined in Section 11.1.1 satisfies (11.3).

Informative Extensions

An interval extension can be more or less informative. In the unary case of a function f : R → R, the most informative extension F is the convex hull of the image of f :

F : I → I x → hull({f (t) | t ∈ x})
Thanks to (11.2), for any interval extension G of f and any x, F (x) ⊆ G(x).

In the other direction, the least informative interval extension of a total function f : R → R is certainly

F : I → I x → (-∞, +∞).

Correlation

The approach of naive interval arithmetic cannot keep track of correlations between subexpressions and might compute overestimated enclosures which are thus useless for proving some goals. This is true even in the case when each operator in the definition of E is overloaded with the most informative interval extension.

For instance, assume that x ∈ [3; 4], so -x ∈ [-4; -3] using the (optimal) interval extension of the negation, so x + (-x) ∈ [3 + (-4); 4 + (-3)] using the (optimal) interval extension of the addition. If the goal was to prove that x -x is always 0, the interval [-1; 1] obtained by naive interval arithmetic is useless.

Interval Arithmetic and CoqInterval

The CoqInterval library CoqInterval is a Coq library for computing numerical enclosures of real-valued expressions [START_REF] Martin | Proving tight bounds on univariate expressions with elementary functions in Coq[END_REF]. It also provides a tactic interval to automatically deduce certain goals from these enclosures.

The tactic typically takes a goal A ≤ e ≤ B where e is an expression in E, and A and B are constants. Using the paradigm of interval arithmetic described in Section 11.1.2, it builds a set e such that e ∈ e holds by construction and such that e evaluates to an interval [inf e; sup e] by computation. Then it checks that A ≤ inf e and sup e ≤ B, again by computation, from which it proves A ≤ e ≤ B. All the computations on interval bounds are performed using a rigorous yet efficient formalization of multi-precision floating-point arithmetic.

Because of the problem of correlation described in Section 11.1.2, the Co-qInterval library also comes with refinements of naive interval arithmetic, such as automatic differentiation and rigorous polynomial approximations, so as to reduce this loss of correlations.

Data structures for reals and intervals

CoqInterval defines R as the set R∪{⊥ R } of extended reals, that is the set of real numbers completed with the extra point ⊥ R . This is an option type: an element of R is either a real number r, or it is a "default value" ⊥ R . The value ⊥ R is typically outputted when an operation is applied outside of its usual domain of definition. For instance, the result of dividing one by zero in R is ⊥ R , while it is unspecified in R. This ⊥ R element is then propagated along the subsequent operations.

Similarly, CoqInterval defines a type I = I ∪ {⊥ I } of extended intervals. An interval operator produces the value ⊥ I whenever the input intervals are not fully included in the definition domain of the corresponding real operator. In other words, an interval operator produces ⊥ I whenever the corresponding operator on R would have produced ⊥ R for at least one value in one of the input intervals. Just like ⊥ R , ⊥ I propagates along computations. The • ∈ • membership relation is hence extended: Definition 12 (membership in R and I). The elements ⊥ R and ⊥ I have the following properties with regard to membership:

• ⊥ R / ∈ [u, v] for all u, v ∈ R ∪ {-∞, +∞}; • x ∈ ⊥ I for all x ∈ R.
In particular, a crucial remark is that ⊥ I is the only interval containing ⊥ R .

The implementation of interval arithmetic on which we build uses pairs of floating-point numbers [inf x; sup x]. Unless otherwise specified, we will be identifying I and I from now on.

Representation of functions, efficient computations: straight-line programs

The CoqInterval library represents expressions in E as straight-line programs. This allows for explicit sharing of common subexpressions. Such a program is just a list of statements indicating what the operation is and where its inputs can be found. The place where the output is stored is left implicit: the result of an operation is always put at the top of the evaluation stack. The evaluation model is simple: the stack grows linearly with the size of the expression since no element of the stack is ever removed. The stack is initially filled with values corresponding to the constants of the program. The result of evaluating a straight-line program is at the top of the stack.

Below is an example of a straight-line program corresponding to the expression (t + π) √ t -(t + π). It is a list containing the operations to be performed. Each list item first indicates the arity of the operation, then the operation itself, and finally the relative depth at which the inputs of the operation can be found in the evaluation stack. Note that, in this example, t and π are seen as constants, so the initial stack contains values that correspond to these subterms 1 . The comments in the term below indicate the content of the evaluation stack before evaluating each statement. The evaluation of a straight-line program depends on the interpretation of the arithmetic operations and on the values stored in the initial stack. For instance, if the arithmetic operations are the operations from the Reals library (e.g. Rplus) and if the stack contains the symbolic value of the constants, then the result is the actual expression over real numbers. If u is a unary operator, we will denote by u R the real interpretation of u.

So if u is the Sqrt operator, u R = √ •.
Similarly, if b is a binary operator, then b R will denote its real interpretation, so that if b is Add, b R = +. We will use a similar notation u A and b A for other evaluation schemes into a set A.

Let us denote p R ( x) the result of evaluating the straight-line program p with operators from Reals over an initial stack x of real numbers. Similarly, p I ( x) denotes the result of evaluating p with interval operations over a stack of intervals. Then, thanks to the inclusion property of interval arithmetic, we can prove the following formula once and for all:

∀p, ∀ x ∈ R n , ∀ x ∈ I n , (∀i ≤ n, x i ∈ x i ) ⇒ p R ( x) ∈ p I ( x). (11.4)
In particular, Formula (11.4) implies that if a function f can be reified into

t → p R (t, x), then t → p I (t, x) is an interval extension of f if ∀i, x i ∈ x i .
Formula (11.4) is the basic block used by the interval tactic for proving enclosures of expressions [START_REF] Martin | Proving tight bounds on univariate expressions with elementary functions in Coq[END_REF]. Given a goal A ≤ e ≤ B, the tactic first looks for a program p and a stack x of real numbers such that p R ( x) = e. This process of finding a program p is called reification. It is done by a tactic outside the logic of Coq, so that the fact that p R ( x) = e holds is checked by Coq. More precisely, the goal A ≤ e ≤ B is convertible to p R ( x) ∈ [A; B] if A and B are floating-point numbers and if the tactic successfully reified the term.

The tactic then looks in the context for hypotheses of the form A i ≤ x i ≤ B i so that it can build a stack x of intervals such that ∀i, x i ∈ x i . If there is no such hypothesis, the tactic just uses (-∞; +∞) for x i . The tactic can now apply Formula (11.4) to replace the goal by p I ( x) ⊆ [A; B]. It then attempts to prove this new goal entirely by computation. Note that even if the original goal holds, this attempt may fail due to loss of correlation inherent to interval arithmetic as seen in Section 11.1.2. Although naive interval arithmetic is useful for approximating real expressions, we have seen in Section 11.1.2 that it suffers from loss of correlation. A less crude way of approximating real numbers on an interval is the use of uniform approximations of functions by polynomials:

Rigorous Polynomial Approximations

Definition 13 (Rigorous Polynomial Approximations). Let x ∈ I and let f :

x → R be a function which is n+1 times differentiable.

A rigorous polynomial approximation (RPA) for a f is a pair (P, ∆) where

P ∈ R n [X] and ∆ ∈ I, such that ∀x ∈ x, f (x) -P (x) ∈ ∆.
RPAs as defined in Definition 13 are agnostic as to the way in which they are obtained. However, the usual privileged method of producing them is by truncating a Taylor series expansion of the function f around a point x 0 ∈ x. This method uses an interval version of the Taylor-Lagrange formula to estimate the remainder:

∀x ∈ x, f (x) - n i=0 f (i) (x 0 ) i! (x -x 0 ) i ∈ f (n+1) (x -x 0 ) (n + 1)! (x -x 0 ) n+1 (11.5)
Such a RPA is called a Taylor Model [START_REF] Makino | Taylor models and other validated functional inclusion methods[END_REF]. The CoqApprox [START_REF] Martin | Proving tight bounds on univariate expressions with elementary functions in Coq[END_REF] library is an extension to CoqInterval which provides a generic theory of RPAs. It also comes with an formalized effective implementation which can build Taylor-Model-based RPAs for univariate expressions in the set E described in Section 11.1.1. The formalization is based on the algorithms developed in the PhD thesis of Joldes [START_REF] Mioara | Rigorous Polynomial Approximations and Applications[END_REF].

Since this formalization contains an effective version of RPAs using interval arithmetic, the actual data types used during computations are different from those in the idealized context of Definition 13: polynomials have interval coefficients, f takes its values in R to accommodate for its domain, x 0 becomes an interval X 0 . The CoqApprox library defines a predicate specifying what it means for such a Rigorous Polynomial Approximation to be correct with respect to a function it approximates: Definition 14 (Validity predicate for RPAs in CoqApprox). Let P := n i=0 α i X i be a degree-n polynomial with interval coefficients and let ∆ ∈ I.

Then (P, ∆) is a valid RPA for f : R → R on I around X 0 if 1. f is defined on I, that is, ∀x ∈ I, (f (x) = ⊥ R =⇒ ∆ = ⊥ I ) ; 2. I = ⊥ I =⇒ ∆ = ⊥ I ; 3. 0 ∈ ∆; 4. X 0 ⊆ I; 5. ∀x 0 ∈ X 0 , ∃Q ∈ R n [X], Q ∈ P ∧ ∀x ∈ I, f (x) -Q(x -x 0 ) ∈ ∆.
where

Q ∈ P means that if Q = n i=0 a k X k , ∀i, 0 ≤ i ≤ n =⇒ a i ∈ α i .

Primitives of RPAs

The following work is a contribution and was integrated into CoqInterval.

We place ourselves within the context of Definition 13 in order not to make the exposition cumbersome with interval coefficient polynomials. If T = (P, I) is a valid RPA for f : X → R with respect to X 0 ⊂ X, we want to build a valid RPA U = (Q, J) for x → x x0 f (s)ds with respect to the same X 0 ,X. We write P = n k=0 a k X k , and for all c ∈ R, we define prim(P, c) by

prim(P, c) = c + n k=0 a k k + 1 X k+1
In other words, prim(P, c) is the only primitive of P which evaluates to c at 0. Consider R = prim(P, 0). If

x 0 , x 1 ∈ X 0 , then ∀x ∈ X, x x0 f (s) ds -R(x -x 1 ) = x x0 f (s) ds -[R(s -x 1 )] x x1 = x x0 f (s) ds - x x1 P (s -x 1 ) ds = x x1 (f (s) -P (s -x 1 )) ds + x1 x0 f (s) ds but for s ∈ [x 0 , x 1 ] 2 , thanks to Definition 13, f (s) -P (s -x 0 ) ∈ I so that f (s) ∈ P (s -x 0 ) + I and then x1 x0 f (s) ∈ x1 x0 P (s -x 0 ) ds + (X 0 -X 0 )I and thus x1 x0 f (s) ds ∈ R(X 0 -X 0 ) + (X 0 -X 0 )I
which brings us to the inclusion:

∀x ∈ X, x x0 f (s) ds -R(x -x 1 ) ∈ (x -x 0 )I + R(X 0 -X 0 ) + (X 0 -X 0 )I
and finally in all generality ∀x ∈ X,

x x0 f (s) ds -R(x -x 1 ) ∈ (X -X 0 )I + R(X 0 -X 0 ) + (X 0 -X 0 )I so that (R, (X-X 0 )I +R(X 0 -X 0 )+(X 0 -X 0 )I) is an RPA for x → x x0 f (s) ds on (X 0 , X).
Notice that if X 0 is thin (which is a reasonable scenario, as we are developing around X 0 ), the error reduces to (X -X 0 )I which is simply the scaling of the error by the size of the interval.

Polynomial Approximations Are Not Always Better

Polynomial approximations usually give tighter enclosures of an integral v u f than the naive approach. In fact, the naive approach is an instance of the polynomial one: let y be the midpoint of F(hull(u, v)) and let ∆ := F(hull(u, v))-y. Then (y, ∆) is a polynomial approximation of degree 0 for f since for all t ∈ [u; v], f (t) ∈ F(hull(u, v)) and thus f (t) -y ∈ ∆.

However, due to the way polynomial approximations are obtained in Coq-Approx, it can happen that a polynomial approximation is less accurate than the use of naive interval arithmetic.

One first instance could be around a singularity preventing the Taylor series of f to exist at all points of [u; v], as in the example of x → √ x on [0; ε], with ε > 0. In this case, CoqApprox returns a pair (P, ⊥ I ) which yields the uninformative enclosure ⊥ I for ε 0

√

x dx, whereas the naive estimate would have given (ε -0)

• [0; √ ε] = [0; ε √ ε].
Another instance would be if we wanted to integrate a strongly oscillating function over a large interval. In Section 12.3, we will be looking at the integral

8 0 sin(x + exp (x)) dx
Indeed, let us use the software Sollya [START_REF] Sylvain Chevillard | Sollya: An environment for the development of numerical codes[END_REF], which implements the same algorithm as CoqApprox, to find a Taylor Model for the integrand on the whole interval [0; 8]. In Figure 12.1, we use the command taylorform which takes as input the function f , the approximating degree n, the point around which to develop the Taylor series x 0 , and the interval I on which to find the approximation. It returns a Taylor model (P, ∆). We only look at ∆ for our purposes here, in the cases of n = 10 and n = 50. In the case of n = 10, the error ∆ has width greater than 2 • 10 30 ; when we try to increase the degree to 50, it is greater than 2 • 10 111 . The reason is that the n-th derivative of this function is big and the error interval obtained from Equation (11.5) is not tight. Lemma 27 is going to yield a useless estimate. Now if we use a naive estimate, we get the much better enclosure

v u f = x1 x0 f + . . . + xn xn-1 f with x i = u + i n (v -u).
Let w(x) = sup x -inf x denote the width of an interval. The smaller w(x) is, the more accurately any real x ∈ x is approximated by x. Any sensible interval arithmetic respects w(x + y) w(x) + w(y) and w(k

• x) k • w(x).
We consider the case of the naive approach first. We assume that F is an optimal interval extension of f and that f has a Lipschitz-constant equal to k 0 , that is, w(F (x))

k 0 • w(x). Since w(naive([x i ; x i+1 ])) (x i+1 -x i ) • w(F ([x i ; x i+1 ]
)), we get the following accuracy when computing the integral:

w i naive([x i ; x i+1 ]) i (x i+1 -x i ) • k 0 • w(([x i ; x i+1 ])) (12.
3)

k 0 • (v -u)/n • i w(([x i ; x i+1 ])) (12.4)
k 0 • (v -u) 2 /n. (12.5)
Suppose we want to gain one bit of accuracy by going from n 1 integrals in the above sum to some number n 2 of integrals. This means that we want

k 0 •(v-u) 2 /n 2 ≤ 1 2 •k 0 •(v-u) 2 /n 1 , so that n 2 ≥ 2 n 1
, so that we have to multiply the computation time by at least two, hence the complexity is exponential. Now for the polynomial enclosure. Let us assume we can compute a polynomial approximation of f on any interval x with an error ∆(x). Because of Formula (11.5), we can expect this error to satisfy w(∆(x)) k d • w(x) d+1 with d the degree of the polynomial approximation and k d depending on f . Since

w(poly([x i ; x i+1 ])) (x i+1 -x i ) • w(∆([x i , x i+1 ])),
developing as in (12.3), the accuracy is now

w i poly([x i ; x i+1 ]) k d • (v -u) d+2 /n d+1 . For a fixed d, k d • (v -u) d+2 /n d+1 2 ≤ 1 2 • k d • (v -u) d+2 /n d+1 1 implies that n 2 ≥ 2 1
d+1 n 1 : one still has to increase n exponentially with respect to the target accuracy. The power coefficient, however, is much smaller than for the naive method. By doubling the computation time, one gets d + 1 additional bits of accuracy.

In order to improve the accuracy of the result, one can increase d instead of n. If f behaves similarly to exp or sin, Taylor-Lagrange formula tells us that k d decreases as fast as (d!) -1 . Moreover, the time complexity of computing a polynomial approximation usually grows like d 3 . If e.g. n v -u, doubling the computation time is thus done by increasing d to d by a factor of about 2

1 3
1.25. This gives about 25% more bits of accuracy. As can be seen from the considerations above, striking the proper balance between n and d for reaching a target accuracy in a minimal amount of time is difficult, so we have made the decision of letting the user control d (see Section 12.1.4) while the implementation adaptively splits the integration interval.

Automatic Continuity and Integrability

When computing the enclosure of an integral, the tactic should first obtain a formal proof that the integrand is integrable on the integration domain, as this is a prerequisite to all the theorems in Section 12. In fact we can be more clever: we prove that, if we succeed in numerically computing an informative enclosure of the integral, the function was actually integrable. This way, the tactic does not have to prove anything beforehand about the integrand.

For this, we use the data type for extended reals R = R ∪ {⊥ R } defined in Section 11.1.3. In order to gain intuition of what we are trying to do, let us look at a simple example first.

Suppose we have a function ln : R → R such that

∀x ∈ R, x > 0 =⇒ ln x = ln x; (12.6) ∀x ∈ R, x ≤ 0 =⇒ ln x = ⊥ R ; (12.7) ln(⊥ R ) = ⊥ R . (12.8)
Suppose that ln : I → I is an extension of ln and ln(

[u; v]) = ⊥ I , that is ln([u; v]) = [u ; v ] for some u , v ∈ R. Since ⊥ I [u ; v ], ⊥ R / ∈ [u ; v ]
because of Definition 12. Furthermore, we can deduce from Equations (12.6), (12.7) and (12.8) that [u; v] ⊆ (0; +∞). Hence ln is defined on [u; v], and since ln is continuous on its domain of definition, we can deduce that ln is continuous on [u; v]. We just deduced from computation that a function was continuous. As any continuous function on a compact interval is also integrable on that interval, we will be using this principle to automatically deduce integrability.

Let us make this more formal. We already saw two evaluation schemes p R and p I in Section 11.1.4. The alternate scheme p R produces the value ⊥ R as soon as an operation is applied to inputs that are outside the usual definition domain of the operator (this would be a way to define ln in the example above). This ⊥ R element is then propagated along the subsequent operations. Thus, the following equality holds, using the trivial embedding from R into R:

∀p, ∀ x ∈ R n , p R ( x) = ⊥ R ⇒ p R ( x) = p R ( x).
(12.9)

Taking advantage of the properties of ⊥ R and ⊥ I in Definition 12, we can also prove a variant of Formula (11.4):

∀p, ∀ x ∈ R n , ∀ x ∈ I n , (∀i ≤ n, x i ∈ x i ) ⇒ p R ( x) ∈ p I ( x). (12.10)
In CoqInterval, Formula (11.4) is actually just a consequence of both Formulas (12.9) and (12.10). This is due to two properties of ⊥ I . First, (-∞; +∞) ⊆ ⊥ I holds, so the conclusion of Formula (12.10) trivially holds whenever p I ( x) evaluates to ⊥ I . Second, ⊥ I is the only interval containing ⊥ R . As a consequence, whenever p I ( x) does not evaluate to ⊥ I the premise of Formula (12.9) holds.

Let us go back to the issue of proving integrability. By definition, whenever p R ( x) does not evaluate to ⊥ R the inputs x are part of the definition domain of the expression represented by p. But we can actually prove a stronger property: not only is x part of the definition domain, it is also part of the continuity

Catalog of supported integrable functions

In order to be able to use Lemma 30, we need to be able to find a suitable extension I g : I → I for the improper part of the integral of g. We thus look at two scales of well-known integrable functions for which we can build I g : Bertrand functions (Section 12.2.2) and exponential functions (Section 12.2.2).

Bertrand integrals

We consider functions g(x) = x α ln β x with α ∈ R, β ∈ R, and which from now on we will call Bertrand functions. These functions are of constant positive sign on [1; +∞). They are integrable at +∞ only when α < -1, or when α = -1 and β < -1. Now we focus on how to compute them. If α < -1, β = 0 and u > 0,

+∞ u x α dx = - u α+1 α + 1 . ( 12.13)
When β ≥ 1, integrating by parts shows that

+∞ u x α ln β x dx = - u α+1 ln β u α + 1 - β α + 1 +∞ u x α ln β-1 x dx. (12.14)
Note that in order to prove this identity, we had to extend Coquelicot with a proof of the general formula for integration by parts (more details in Section 12.2.5). When α < -1 and β < 0, there is no closed form, but by moving ln β x into the bounded part of Lemma 29, we can nevertheless compute bounds on the integral. If f is bounded on [u; ∞], then

+∞ u f • x α ln β x dx = +∞ u f ln β x bounded •x α dx
and then Formula (12.13) can be used in Lemma 29.

When α = -1 and β < -1, we have a closed form:

+∞ u ln β x x dx = - ln β+1 u β + 1 .
When β ∈ N, using Equations (12.13) and (12.14), we can recursively obtain a closed form to evaluate Bertrand integrals with α < -1 and β ≥ 0. For instance, using (12.14) then (12.13), we get:

+∞ 1 ln x x 2 dx = - 1 -1 ln 1 -1 - 1 -1 +∞ 1 dx x 2 = 0 + (-) 1 -1 -1 = 1.

Exponential

We also handle the case of the positive function g(x) = e γx with γ < 0, using the fact that +∞ u e γx dx = -e γu γ .

Case of 0 +

When the singular bound is 0 + instead of +∞, we use a variant of Lemma 29.

Lemma 31. Let f, g : R → R. Suppose that, on (0; v], f is bounded, f and g are continuous, and g has a constant sign. Moreover, suppose that v 0 + g exists. Then v 0 + f g exists, and

v 0 + f g ∈ hull{f (t) | 0 ≤ t ≤ v} • v 0 + g.
As in the case of +∞, we provide a catalog of supported functions. Consider g(t) = t α (-ln t) β with α ∈ R, β ∈ R. This function is of constant sign on (0; v], where v < 1. Observe that using the substitution t = 1

x , we get

v 0 + t α (-ln(t)) β dt = ∞ 1/v
x -2-α ln β x dx.

The right-hand-side integral has the shape treated in Section 12.2.2, so we have a way to bound the left-hand-side integral. To do so, we added a proof of the substitution lemma to Coquelicot (see Section 12.2.5).

Note that the functions recognized by the CoqInterval library are actually of the shape t α ln β t, with α and β integers.

Integrability

Just like in the case of proper integrals in Section 12.1.6, we deduce integrability from computation.

As mentioned above, for improper integrals v u f g, the enclosure of the function f has to be bounded. This is not always very easy to obtain in the case of an integration domain extending to +∞. Indeed, the input domain x is no longer bounded in that case, which means that RPAs become useless. Let us consider the following integral to illustrate the issue (u > 0):

+∞ u x + 1 x + 2 e -x dx.
The quotient is bounded on [u; +∞). Yet using naive interval arithmetic gives [u+1; +∞)/[u+2; +∞) = [0; +∞), which is not bounded. Thus the tactic is unable to prove integrability and to compute an enclosure of the integral. To circumvent this issue, the user has to massage the bounded part of the integrand into a form suitable for naive interval arithmetic, e.g. (1 + 1/x)/(1 + 2/x). This time, the tactic obtains [1/(1 + 2/u); 1 + 1/u], which is bounded. Note that this kind of transformation of the integrand is not always possible. Consider for example the integral

∞ 1 ln x x e -x dx
One would hope to apply Lemma 29 by proving that f (x) := ln x x is bounded on [1; ∞] and using e -x as the function g. However, given the primitives of E there is no way to massage the expression ln x

x so that naive interval arithmetic sees is as bounded. In fact, we can not compute an enclosure of this integral automatically as of now.

We conclude this part about integrability with a remark about the fact that the only functions integrable at ∞ that we handle are bounded, whereas of course there exist unbounded integrable functions. One might be tempted to think that the only functions in the class E which are integrable at +∞ are necessarily bounded. This is actually false, as the following example suggested by Philippe Dumas shows:

∞ 0 t 1 + sin 2 t 2 t 5
dt < ∞.

Extending Coquelicot

In order to be able to handle improper integrals, a nontrivial amount of work was needed to improve the existing results provided by the Coquelicot library. This section dives into more technical details about implementations of concepts in Coq. Although some parts of this section may seem at times to be overly focused on proof-engineering details, these reflect the kind of challenges which arise when one wants to extend a library of proofs and tactics. First, we will look at the way the Coq standard library defines real numbers and the basic notions of Real Analysis. Then, we will examine the Coquelicot library which extends the standard library and improves it usability. Finally, we will explain the work which was necessary to extend the tactic to improper integrals.

The Coq standard library of real numbers

Coq's standard library Reals1 axiomatizes real arithmetic with a classical flavor [START_REF] Mayero | Formalisation et automatisation de preuves en analyses réelle et numérique[END_REF]. This means that a type R of real numbers is declared as an axiom (as opposed to being constructed) and that the elements 0, 1, the operators +, * ,/,• -1 are assumed as well as a list of axioms. For example, the comparison operator < is assumed to be total:

Axiom total_order_T : ∀ r1 r2:R, {r1 < r2} + {r1 = r2} + {r1 > r2}.

This means that given two reals x and y, it is decidable to compare them, i.e. there exists a Coq function whose result determines either x < y, x = y, or x > y. This is a very strong assumption since in practice, for any given concrete representation of real numbers, such a function can't be built. The library also contains the following axiom: The completeness lemma states that any bounded non-empty subset of R has a computable least upper-bound. Note that this is both very strong (it postulates a Coq program which takes as input any bounded non-empty subset of R and returns its upper-bound u) and weaker than what is typically assumed in textbook mathematics: this axiom is not enough to obtain a sequence (u n ) n∈N of elements of E converging to u.

The standard library comes with a formalization of the properties of usual mathematical functions like sin, cos, exp, and so on, defined as series. It also provides some notions of elementary real analysis, including the definition of continuity, differentiability, and Riemann integrability:

Definition RiemannInt (f:R -> R) (u v:R) (pr:Riemann_integrable f u v) : R.

This definition uses a dependent type: one needs a proof pr that f is Riemann-integrable between u and v in order to write the object v u f . If one wants to add two integrals v u f and v u g, one first needs to come up with a proof that f + g is integrable between u and v before one can even write down the integral v u f + g.

Coquelicot

This section describes the Coquelicot library as it stood before the present work. This library is a conservative extension of the theory of real numbers from the Coq standard library, meaning that it uses the same axiomatization of real numbers and does not add any axiom of its own.

Coquelicot introduces the notion of filters in order to manipulate limits in a uniform (e.g. between sequences, functions and functionals) and elegant (with less δ/ε manipulations) way; this idea is borrowed from a similar scheme in Isabelle/HOL [START_REF] Hölzl | Type Classes and Filters for Mathematical Analysis in Isabelle/HOL[END_REF]; the notion dates back to Cartan [START_REF] Cartan | Filtres et ultrafiltres[END_REF] [START_REF] Bourbaki | Eléments de mathématique : Livre III. Topologie générale[END_REF]. A filter is a family of sets with properties related to inclusion (see Definition 16). Filters allow to treat in similar ways various kinds of limits such as sequence limits, limits of functions at one point or on one side of a point, and asymptotic limits. They can in fact be thought of as a way to encode the familiar notion of neighborhoods in analysis.

Definition 16 (filter). Let E be a set. F ∈ P(E) is a filter of E if:

• E ∈ F ; • For all A, B ⊂ E, if A ∈ F and B ∈ F then A ∩ B ∈ F ; • If A ∈ F and A ⊂ B, then B ∈ F .
Here is an example of filter containing the neighbourhoods of ∞:

Definition 17. Let f ∞ := {E ⊂ R | ∃M ∈ R, ∀x ≥ M, x ∈ E}. Then f ∞ is a filter of R.
Coquelicot has an inductive type Rbar encoding reals extended with ∞ and -∞: A filter encoding a neighborhood of x ∈ Rbar can be obtained by evaluating Using such statements would have made it quite cumbersome for us to extend our tactic to generalized integrals.

New definitions

In order to be able to state and prove in Coq the theorems of Sections 12.2.1, 12.2.2, 12.2.3 and 12.2.4 as well as to extend the tactic to improper integrals, we needed a total operator computing such integrals and some more lemmas of Real Analysis on them.

A new definition of is_RInt_gen: We redefined is_RInt_gen in a more natural way in order to simplify proofs. If a ∈ R and f : R → R, a natural way to define

∞ a f is ∞ a f = lim b→∞ b a f,
where this equality is meant as: "there exists a neighborhood U of ∞ such that when b ∈ U , b a f exists and lim b∈U,b→∞ b a f exists". This is so cumbersome to express because the function (f, a, b) → b a f itself is partial. Since the generalized integral is in some sense the "limit" of this relation, this motivated the introduction of the notion of limit of a partial function defined implicitly by a predicate:

Definition filterlimi {T U : Type} (R : T -> U -> Prop) F G : Prop.
For any relation R on elements x of type T and y of type U, and for any filters F and G, the predicate filterlimi R F G states that for any neighborhood Q ∈ G, the set {x ∈ T |∃y ∈ Q, R x y} of "antecedents" of Q by R is itself a neighborhood of F . Now we are ready to read the definition of is_RInt_gen:

Definition is_RInt_gen (f : R -> V) (Fa Fb : (R -> Prop) -> Prop) (l : V) := filterlimi (fun (a,b) => is_RInt f a b) (filter_prod Fa Fb) (locally l).
This is saying that for any function f : R → V, for any real numbers a and b, for any neighborhood W of l, for any w ∈ W, the set of pairs (a,b) such that is_RInt f a b w holds is in the filter product of Fa and Fb.

A Total Operator for Generalized Integrals: The following non-computable iota function is useful for building a total function from a predicate on a complete space T: Definition iota (P : T -> Prop) : T := lim (fun A => (∀ x, P x -> A x)).

The function lim is here given as input a filter N formed of all sets A containing P. Assuming 1) that N does not contain the empty set and 2) that for all ε, there is some x such that the ball B(x, ε) ⊆ N, lim N returns an element x ∈ T such that all open balls centered on x are in N. In the particular case above, if one can prove that any two points in P are equal (which is the case for is_RInt or is_RInt_gen, for instance), iota is effectively a choice function for the input property. If the property P is true for exactly one value x of the complete space T, then iota P returns this value; otherwise, it returns a default value.

Using iota, we defined the function Just like for RInt, the output of this function only makes sense if f is integrable between the bounds a and b. Notice that RInt_gen is more general than RInt: it takes as argument a function whose values are in a complete normed R-module. We redefined RInt in a similar way (defined on more general functions):

Definition RInt (f : R -> V) (a b : R) := iota (is_RInt f a b).

Unanticipated Consequences of the New Definitions

As we saw above, the new definition of RInt using iota accepts functions with a more general type. This allows for example to define proper integrals of complexvalued functions. As we defined it in this more general way, we realized that some elementary theorems about RInt had been proven using arguments which are specific to R. For instance, and that F (a) = 0, so that ∀x ∈ [a, b], F (x) = 0 hence the result. However, it turns out that the following lemma, implicitly used in this usual proof, was not proved in Coquelicot with the necessary generality: This lemma states that if the derivative of a function f is uniformly zero in [a; b], then f is a constant function. It is in fact not so trivial as it may seem at first, because the first simple proof which comes to mind relies on the theorem is_RInt_derive which we are trying to prove. Our proof relies on the completeness axiom of the Reals library presented in Section 12.2.5 which states the existence of a least upper bound for bounded and non-empty subsets of R. Here is our proof of eq_is_derive:

Proof. Let ε > 0. We want to prove that for all

x ∈ [a, b], |f (x) -f (a)| ≤ ε. Let A = {x ≤ b | ∀t, a ≤ t ≤ x =⇒ |f (t) -f (a)| ≤ ε(t -a)}.
A is non-empty since a ∈ A, and A is certainly bounded by b, so that A has an upper bound u. Let c ∈ A be such that c < b.

Since f (c) = 0, we can find δ > 0 so that c + δ ≤ b and In order to finish the proof of is_RInt_derive, we also needed a proof of the fact that the derivative of

∀x ∈ [c -δ; c + δ], f (x) -f (c) ≤ (x -c)ε If x ∈ [c; c + δ], since c ∈ A, |f (x) -f (a)| ≤ |f (x) -f (c)| + |f (c) -f (a)| ≤ (x -c)ε + (c -a)ε = (x -a)ε.

Listing 12.1: Derivative of an integral in Coquelicot

The first hypothesis Hloc is very strong. Applied to our specific case, and assuming a ≤ b, it mandates that x a f (x)dx exist in a neighborhood of b, that is a little to the left of b (that we can do with, especially when a < b) but also to the right of b. This has no chance to be true in our case, as the following counterexample shows: Take a = -1, b = 0, and f (x) = x on [-∞, 0], and

f (x) = xχ Q (x) on [0, ∞],
where χ Q is the characteristic function of Q on R, i.e χ Q returns 1 on rational numbers and 0 on irrational numbers. Then f satisfies all the preconditions of is_RInt_derive, as it is C 1 at all points of [a, b] = [-1, 0]:

• f is C 1 trivially at any point in [-1, 0); • f is continuous at any point of [-1, 0].
• f is defined and is continuous to the left of 0, with f left (0) = 1;

• For x ≥ 0, f (x) -f (0) -1 • x = x(χ Q (x) -1) → 0 when x → 0, hence f right (0) = 1 and thus f (0) = 1.
However it can't be that f (x) =

x a df (x)dx in a neighborhood of 0, as the quantity on the right is continuous on that neighborhood while f is not, at least to the right of 0.

Therefore, how can we use is_derive_RInt in the proof of is_RInt_derive? We want ∀x ∈ [a; b],

x a df (x)dx = df (x) (12.18) However, in the hypotheses of is_RInt_derive, the function f is continuously differentiable on [a; b], so that df is continuous, which is weaker than what the Hloc hypothesis of is_derive_RInt mandates. Let g be a C 1 extension of f on [a; b], meaning that f and g coincide on [a; b], and g is affine on (-∞; a] and on [b; +∞), and C 1 on R. Then since g is continuous everywhere, x → x a g (x)dx is defined everywhere so that Hloc is true for x → But g and df coincide on [a; b], which entails (12.18).

Substitution lemma

We also (re-)proved the substitution lemma for proper integrals for the same reason as is_RInt_derive: the proof was using the Intermediate Value Theorem which is specific to real-valued functions and was broken by the generalization of the type of is_RInt. This lemma was crucial for handling Bertrand integrals at 0 + (see Section 12.2.3). Just like in Section 12.2.5, this proof uses a continuously differentiable extension of the function g.

Conclusion

The changes described in this section are now part of Coquelicot 3.0.0 and later versions.

Benchmarks

This section presents the behavior of the tactic on several integration problems, each given as a symbolic integral, its value (approximate if no closed form exists), and a set of absolute error bounds that must be reached by the tactic. Each problem is translated into a set of Coq scripts as follows, one for each bound:

Goal Rabs (RInt[_gen] function domainvalue) <= error . interval with options. Qed.

The tactic options have been set using the following experimental protocol. First, the target relative accuracy is computed from the error bound and the initial estimation of an integral. The floating-point precision is then set at about 10 more bits than the target accuracy, so that round-off errors do not make interval enclosures too large. The maximal depth is originally set to a large enough value. Then, various degrees of RPAs are tested and the one that leads to the fastest execution is kept. Finally, the maximal depth is reduced as long as the tactic succeeds in proving the bounds, so that we get an idea of how deep splitting has to be performed to compute an accurate enclosure of the integral. Note that reducing the maximal depth might improve timings in case the adaptive algorithm had been overly conservative and did too much domain splitting. Reducing the target relative accuracy could also improve timings (again by preventing some domain splitting), but this was not done.

The tables below indicate, for each error bound, the time needed and the tactic settings. Timings are in seconds and are obtained on a standard-grade laptop. All timings are obtained with a version of the tactic that uses the vm compute machinery to perform computations. Modern versions of Coq also have a native compute machinery [START_REF] Boespflug | Full reduction at full throttle[END_REF], which can improve the longest computation times by an order of magnitude, at the expense of a long initialization time which makes it unsuitable for the present benchmarks.

Proper integrals

For each proper integral, we also ran several quadrature methods from Octave [START_REF] Eaton | GNU Octave version 3.8.1 manual: a high-level interactive language for numerical computations[END_REF]: quad, quadv, quadgk, quadl, quadcc. We also used IntLab [START_REF] Rump | Verification methods: Rigorous results using floatingpoint arithmetic[END_REF]; it provides verifyquad, an interval arithmetic procedure that computes integral enclosures using a verified Romberg method. For each method, we ask for an absolute accuracy of 10 -15 . We only comment when the answer is off, or when the execution time exceeds 1 second. Finally, we also tested VNODE-LP [START_REF] Nedialko | Interval tools for ODEs and DAEs[END_REF] on each example by representing the integral as the value of the solution of a differential equation.

The first problem is the integral of the derivative of arctan, a highly regular function. As expected, the tactic behaves well on it, since it takes about 3 seconds to compute 18 decimal digits of π by integration. Note that the time needed for reifying the goal and performing the initial computations is incompressible, so there is not much difference between 10 -3 and 10 -6 . The second problem is Ahmed's integral [START_REF] Ahmed | Ahmed's integral: the maiden solution[END_REF]. It is a bit less regular and uses more operators than the previous problem, but the tactic still behaves well enough: adding ten bits of accuracy doubles the computation time. The third problem involves a function that is harder to approximate using RPAs, so the tactic performs more domain splitting, degrading performances.

π 0 x sin x 1 + cos 2 x dx = π 2 4
Error Time Accy Degree Depth Prec 10 The fourth problem is an example from Helfgott [START_REF] Helfgott | Rigorous numerical integration[END_REF] in the spirit of [START_REF] Helfgott | Major arcs for Goldbach's problem[END_REF]. The polynomial part crosses zero, so there is a point where the integrand is not differentiable because of the absolute value. Thus only degenerate Taylor models can be computed around that point. Although the tactic has to perform a lot of domain splitting to isolate that point, it still computes an enclosure of the integral quickly. Note that the approximate value of the integral was computed using the interval_intro tactic. The last two problems are inherently hard to numerically integrate. The first one is the 12-th coefficient of a Chebyshev expansion. Note that the initial estimation of the integral is completely off, which explains why the relative accuracy has to be set about 30 bits higher than one would expect. As with the previous problem, there are some points where no RPAs can be computed. The approximate value was again computed using the interval_intro tactic.

1 -1 2048x 12 -6144x 10 + 6912x 8 -3584x 6 + 840x 4 -72x 2 + 1 exp -x -3 4 2 1 -x 2 dx -3.2555895745 • 10 -6
The quad, quadl, and quadcc procedures give completely off but consistent answers without warning; quadv gives an answer which is off the mark as well, but it gives a warning "maximum iteration count reached"; verifyquad works only for functions that are four times differentiable, hence its failure here; quadgk gives yet another off answer with no warning. Finally, VNODE-LP fails here because of computational errors such as divisions by 0. The last problem is an example taken from Tucker's book [START_REF] Tucker | Validated Numerics: A Short Introduction to Rigorous Computations[END_REF] and originally suggested by Rump in [109, page 372]. This integral is often incorrectly approximated by computer algebra systems, because of the large number of oscillations (about 950 sign changes) and the large value of the n-th derivatives of the function. While the maximal depth is not too large, the tactic reaches it for numerous subdomains, hence the large computation time.

Error Time Accuracy

The quad, quadcc, and quadgk procedures give off values without any warning; quadv gives an off value with a warning; verifyquad takes 1.7 seconds to give a correct answer; quadl takes 9 seconds to return a correct answer. The second example is similar to the integral from Tucker's book, in the sense that the high number of oscillations of the integrand makes it hard to give an accurate approximation of the integral. For instance, Maple 18 forfeits after 10 seconds of computations. The tactic does not perform much better since it is not able to compute more than two digits in a reasonable amount of time. This is partly due to the adaptive splitting algorithm, which is built upon the assumption that splitting an integration domain into two parts eventually improves the accuracy by more than one bit on each part; this is not the case for the remainder of this example. The tactic cannot handle this integral fully automatically since the integrand is not syntactically a product with a term x α ln β x. It is up to the user to split the integral into two parts: one proper part between -100,000 and 100,000 (as was done in the original paper) and one improper part between 100,000 and +∞ (counted twice, since the integrand is an even function). The proper part is handled in the same way as all the previous examples. It takes about 30 seconds to get the relative accuracy of 10 -6 needed by the original paper. For the improper part, the integrand first has to be transformed into the following form, which was proved to be equal to the original integrand in a few lines of Coq: +∞ 100 000 Chapter 13

1 + 0.5•ln(1+2.
Quadrature and Automatic Differentiation of order n 13.1 Quadrature Methods: A quick introduction

As mentioned in Chapter 10, quadrature methods based on the evaluation of a linear combination of evaluations of the integrand f can be made rigorous under certain assumptions on the regularity of f . In this case, an error bound depending on some iterated derivative of the integrand on the integration interval is typically provided, as in the case of Formula (10.3) for the trapezoidal rule. Another example of a popular and efficient quadrature method is the so-called Simpson's method. Rabinowitz is quoted [START_REF] Davis | Chapter 4 -error analysis[END_REF] claiming that "95% of all pratical work in numerical analysis boiled down to applications of Simpson's rule and linear interpolation":

Theorem 7 (Simpson's method). Let f ∈ C 4 ([u, v]). Then there exists ξ ∈ [u; v] such that v u f (x) dx - b -a 6 f (u) + 4f u + v 2 + f (v) = - (v -u) 5 2880 f (4) (ξ).
In practice, the compound rule is often used:

Theorem 8 (Compound rule version of Simpson's method). Let f ∈ C 4 ([u, v]) and n ∈ N, with n ≥ 1 and N = 2n. Let x 0 = u < x 1 < • • • < x N -1 < x N = v be equally spaced points, with h = v-u N = x i+1 -x i for all i. For 0 ≤ k ≤ N , let f k denote f (x k ). There exists ξ ∈ [u; v] such that v u f (x) dx - h 3 f 0 + 4 n k=1 f 2k-1 + 2 n-1 k=1 f 2k + f 2n In(f ) = - (v -u) 5 180N 4 f (4) (ξ) en(f ) .
As claimed above, most popular integration software uses Simpson's rule or other Newton-Cotes-based methods to compute approximations of integrals.

They usually only make use of I n (f ) and compute it using floating-point approximations (as opposed to interval methods) and ignore the error term e n (f ).

In order to make Theorem 8 effective and rigorous, we need two things. The first requirement is to correctly evaluate f at the points x k , which thanks to the methods exposed in Sections 11.1 and 11.2 we know how to do for at least a class of functions f by building an interval extension F of f which is correct by construction. The second requirement is of a nature which has not been mentioned yet: we need to produce a correct enclosure of f (4) at the point ξ. Since ξ could be any point in [u; v], what we really want is an interval evaluation of f (4) on [u; v]. All this suggests an interval version of the theorem:

Theorem 9 (Interval compound rule version of Simpson's method). Let f ∈ C 4 ([u, v]) and n ∈ N, with n ≥ 1 and N = 2n. Let u, v ∈ I with u ∈ u, v ∈ v.
Let F be an interval extension of f , and let G be an interval extension of f (4) on

[u; v]. For 0 ≤ k ≤ 2n, pose h := v-u 2n and x k := u + k • h and F k := F(x k ). Then v u f (x) dx ∈ h 3 F 0 + 4 n k=1 F 2k-1 + 2 n-1 k=1 F 2k + F 2n In(F) - (v -u) 5 180N 4 G([u; v]) En(G)
.

In order to use Theorem 9 in an implementation, we need a way to estimate the k-th derivative of a function in E for k ≥ 2. The usual technique to obtain numerical approximations of k-th derivatives is called Automatic Differentiation (AD): AD of order k is a numerical method to compute an enclosure of the k-th derivative of a function f . AD of order 1 is already implemented in CoqInterval, but not for higher orders.

Once we obtained results using naive interval arithmetic and RPAs, we compared how far off we are from more traditional quadrature methods. For this, we made several attempts at implementing AD of arbitrary order, which are recounted (Section 13.1.1). Then, a few benchmarks about computation times are presented (Section 13.1.2).

AD of order n

Automatic Differentiation or Algorithmic Differentiation (AD) is a process for computing the numerical derivative of a function. We restrict ourselves to the case of a univariate differentiable real function f : R → R and a real number t 0 , and the goal is to simultaneously compute an interval approximation of f (t 0 ) and f (t 0 ). AD distinguishes itself from two other natural ways of approximating derivatives: the first one is to symbolically compute the derivative of the function when it is presented as an expression made of additions, multiplications, divisions, composition and usual mathematical functions; the second one is to try to approximate f (x) = f (x+h)-f (x) h for some well-chosen h. There exists a large literature on the topic of AD. The combination of AD with interval arithmetic is covered in Moore [START_REF] Moore | Introduction to Interval Analysis[END_REF] and Tucker [START_REF] Tucker | Validated Numerics: A Short Introduction to Rigorous Computations[END_REF] among others. Interval AD of order 1, which is described below, was formalized in CoqInterval [START_REF] Martin | Proving tight bounds on univariate expressions with elementary functions in Coq[END_REF]. To our knowledge, the only other implementation in a formal proof system is in a paper focused on differentiating FORTRAN programs using only basic arithmetic operations [START_REF] Mayero | Using theorem proving for numerical analysis[END_REF].

Here we restrict ourselves to a specific flavor of univariate AD which uses operator overloading. The functions f : R → R that we consider belong to the class E presented in Section 11.1.

AD of order 1

Let us first describe automatic differentiation of order 1. Let e(x) ∈ E be an expression, and let f be the function defined by t → e(t). Suppose we want to compute f (t 0 ) and f (t 0 ) for some t 0 ∈ R such that these quantities are well-defined.

The idea is to build the pair p e (t 0 ) = (f (t 0 ), f (t 0 )) ∈ R 2 recursively on the structure of e(x) ∈ E. The operations are done assuming that f is both defined and differentiable at t 0 : when f (t 0 ) or f (t 0 ) is ill-defined, it is replaced by the ⊥ R value which then propagates. The following rules are used:

• if e(x) = c where c is a constant of E, then p e (t 0 ) = (c, 0);

• if e(x) = x, then p e (t 0 ) = (x, 1);

• if e(x) = e g (x) + e h (x), if p eg (t 0 ) = (v g , v g ) and p e h (t 0 ) = (v h , v h ) then p e (t 0 ) = (v g + v h , v g + v h ); • if e(x) = e g (x) • e h (x), if p eg (t 0 ) = (v g , v g ) and p e h (t 0 ) = (v h , v h ) then p e (t 0 ) = (v g • v h , v g • v h + v h • v g ); • if e(x) = eg(x) e h (x) , then if p eg (t 0 ) = (v g , v g ) and p e h (t 0 ) = (v h , v h ) then p e (t 0 ) = ( vg v h , vg•v h -v h •v g v 2 h
);

• if e(x) = (e g (x)) -1 and p eg (t 0 ) = (v g , v g ), then p e (t 0 ) = (v -1 g , -v g v g -1 2 );

• If e(x) = u f (g(x)), where u f is a usual function (such as • n , exp, log, sin, cos, tan, or any unary function of E) and g

(x) ∈ E, then if p eg (t 0 ) = (v g , v g ), p e (t 0 ) = (u f (v g ), v g • u f (v g )); • If e(x) = |e g (x)| and p eg (t 0 ) = (v g , v g ), then p e (t 0 ) =      (v g , v g ) if v g > 0 (v g , -v g ) if v g < 0 ⊥ R otherwise
The following lemma describes the correctness of this process for constructing p e (t 0 ): Lemma 33. For any expression e(x) ∈ E, for all

t 0 ∈ R if f := t → e(t) is differentiable at t 0 , then p e (t 0 ) = (f (t 0 ), f (t 0 )).
Proof. The proof is by induction on the structure of e(x) ∈ E:

• It is straightforward for the two base cases of c and x;

• For the three binary operators +, •, and /, it simply amounts to the wellknown derivation rules for these operators;

• For unary operators corresponding to a usual function u f , correctness is deduces from the formula for differentiating u f and the chain rule;

• For the absolute value, it is a simple case analysis.

The above description is not yet an effective algorithm, because the data type of the pair p e (x) is R 2 . In order to make an implementation which is able to provide correct numerical enclosures of first-order derivatives, we replicate all the above steps with the usual interval operations and, from an input interval t 0 , build a pair of intervals p e (t 0 ) with the following property:

Lemma 34. For any expression e(x) ∈ E, for all t 0 ∈ I and t 0 ∈ R such that t 0 ∈ t 0 , p e (t 0 ) ∈ p e (t 0 ).

Proof. Since p e (t 0 ) is built the same way as p e (t 0 ), and since every interval operation used is an interval extension of the corresponding real operation, the result follows by induction on the structure of e.

The pair p e (t 0 ) := (i 0 , i 1 ) also has the nice property that if i 1 = ⊥ I , then f := t → e(t) is differentiable at t 0 . This property is similar to the way integrability is inferred in Section 12.1.6.

This process of interval automatic differentiation of order 1 is not a contribution: it is already implemented in CoqInterval and is used for example to refine some interval enclosures by looking at the sign of the derivative.

AD of Order n

It is natural to want to extend the principle of the previous subsection to the nth derivative by computing, still by induction on the structure of the expression e(x) defining f , the vector p n e (x) := (f (x), f (x),• • • , f (n) (x)) ∈ R n+1 . When implementing order 1 automatic differentiation, the only variation in the derivation formulas which might have an impact on numerical accuracy is the associativity and distributivity of operators.

A typical presentation of AD of order n establishes a direct formula for the n-th derivative of each operator. Up to an integer factor, the elements of p n e (t 0 ) are really the n + 1 first coefficients of the Taylor series development of f about t 0 . Formulas for each operation in the definition of E can be inferred from this fact, as explained in Tucker [START_REF] Tucker | Validated Numerics: A Short Introduction to Rigorous Computations[END_REF] for example.

However, in the context of an implementation in a formal proof system such as Coq which we expect to eventually verify, we would like to minimize the effort in proving each and every specific formula and benefit from the elegance of functional programming to obtain this n-th derivative recursively. We would also like for the proof of correctness to only involve facts about derivatives of order 1 when that is possible.

In this subsection, we explore two approaches meeting these criteria. The first one is presented mostly as an oddity which takes advantage of the abstractions of functional programming; however, it is extremely inefficient to implement in practice. The second one is our final implementation; although it has not been proved correct yet, it is both functionally elegant and amenable to formal proof. Expressions in E are left implicit in this part to keep formulas legible, but the functions mentioned should still be understood as elements of E univariate in the same variable x.

A peculiar kind of recursion Let us first describe a surprisingly simple way of doing AD of order n.

Looking back at the previous subsection, what we did was overload common mathematical operators. For example,

• : R × R → R (x, y) → x • y became • 1 : R 2 × R 2 → R 2 ((x 1 , x 2 ), (y 1 , y 2 )) → (x 1 • y 1 , x 1 • y 2 + x 2 • y 1 ).
Let us rename multiplication • as • 0 . We had the nice property that for all x ∈ R, for all f and g differentiable at x,

(f (x), f (x)) • 1 (g(x), g (x)) = (f (x) • 0 g(x), (f (x) • 0 g(x)) )
We can operate this renaming for all operators described in the previous subsection, for instance (x, y) + 1 (z, t) = (x + z, y + t). Now if we pose

• 2 : (R 2 ) 2 × (R 2 ) 2 → (R 2 ) 2 ((x 1 , x 2 ), (y 1 , y 2 )) → (x 1 • 1 y 1 , x 1 • 1 y 2 + 1 x 2 • 1 y 1 )
we can see the following: Let x ∈ R, f, g : R → R twice differentiable at x. For all functions h, let us denote h for h(x), h for h (x). We have:

((f, f ), (f , f )) • 2 ((g, g ), (g , g )) = ((f, f ) • 1 (g, g ), (f, f ) • 1 (g , g ) + 1 (f , f ) • 1 (g, g )) = ((f • 0 g, (f • 0 g) ), ((f • 0 g) , f • 0 g + f • 0 g + f • 0 g + f • 0 g )) = ((f • 0 g, (f • 0 g) ), ((f • 0 g) , (f • 0 g) ))
so that although we have a redundancy in the middle, we now have the derivatives of order 0, 1 and 2 of f • g at x.

In the same way, we could recursively define

• n+1 : (R 2 n ) 2 × (R 2 n ) 2 → (R 2 n ) 2 ((x 1 , x 2 ), (y 1 , y 2 )) → (x 1 • n y 1 , x 1 • n y 2 + n x 2 • n y 1 )
as well as for all the other operators of E. An induction on the index n of the operators would show that it is possible by iterating this construct to obtain a vector of exponential size in n containing all of f (x), f (x),• • • , f (n) (x). We have iterated the formal construction of automatic differentiation from base operations by applying it recursively to itself. Although this is a pretty result, it is not usable in practice, which is why we move on to a more reasonable way of implementing AD or order n.

Building AD of order n more realistically From now on the notation op n introduced in subsection 13.1.1 is superseded by the following. Let us define

• n : R n+1 × R n+1 → R n+1 x, ȳ → (x 0 • y 0 ,• • • , k i=0 k i x i x k-i ,• • • , n i=0 n i x i x n-i ) so that if f , g are two n-times differentiable functions at x ∈ R, we have (f (x), f (x),• • • , f (n) (x)) • n (g(x), g (x),• • • , g (n) (x)) = ((f • g)(x), (f • g) (x),• • • , (f • g) (n) (x)).
In the same way for op ∈ {+, -, /} we can define:

op n : R n+1 × R n+1 → R n+1 so that if f , g are two n-times differentiable functions at x ∈ R, (f (x), f (x),• • • , f (n) (x)) op n (g(x), g (x),• • • , g (n) (x)) = ((f op g)(x), (f op g) (x),• • • , (f op g) (n) (x)).
and similarly for unary operations op ∈ {exp, ln, cos, sin, tan, arctan}:

op n : R n+1 → R n+1 so that f is a n-times differentiable function at x ∈ R, op n (f (x), f (x),• • • , f (n) (x)) = ((op f )(x), (op f ) (x),• • • , (op f ) (n) (x).
If we want to inductively define each operation as in Section 13.1.1, we will need a formula for the k-th derivative of each operation for 1 ≤ k ≤ n. A first idea is to compute these formulas through iterated symbolic derivation: for example, the formula for the 2 first derivatives of multiplication given by this process is:

(x 0 , x 1 , x 2 ) • 2 (y 0 , y 1 , y 2 ) = (x 0 y 0 , x 0 y 1 + x 1 y 0 , (x 0 y 2 + x 1 y 1 ) + (x 1 y 1 + x 2 y 0 ))
This example exposes a major problem with using iterated derivation: expressions can blow up exponentially with n and when evaluating the same computations will be done several times. Furthermore, adding the same expression to itself several times will cause a loss of precision which may make the n-th derivative unusable at best, if it is computed in a rigorous way (through interval arithmetic for example), or completely incorrect at worst if it is computed using floating-point numbers.

An interval version of the symbolic derivation method was programmed in Coq using interval arithmetic and, as expected, it quickly becomes uninformative as n grows. The approach we chose uses a "hardcoded formula" for multiplication in order to avoid this explosion, but tries to be more generic for most operations to facilitate (future) formal proofs. Here is our construction.

• (x 0 ,• • • , x n ) + n (y 0 ,• • • , y n ) := (x 0 + y 0 ,• • • , x n + y n )

Test of Simpson's Method in Coq

In Section 12.3, we compared our implementation against some integration software which uses quadrature methods. This was fruitful for correctness and accuracy, but the speed of interpreted Coq code is still much slower than a well-tuned C implementation. It makes more sense to compare speeds with a Coq implementation of a quadrature method.

Our implementation uses the same strategy of dichotomy as the one described in Section 12.1.4, but uses our implementation of Simpson's method instead of Taylor Models. There are two variants: one which uses the simple version of Simpson's method as in Theorem 7, and one which uses the compound rules with n points as in Theorem 9.

First, we looked at the group of the three first integrals from Section 12.3. They have in common that the integrand is four times differentiable on the whole integration interval. For each of them, we started with the parameters from the corresponding array in Section 12.3 for the line with the best accuracy. The simple Simpson rule being always slower than Taylor Models in practice, we then tried to vary the depth and number of points to reach the fastest possible execution reaching comparable accuracy as the Taylor Model based one. We indicate the target accuracy, depth of the dichotomy, and precision of intermediary computations used to obtain the new result. The column "Points" gives the number of points used in Simpson's compound rule. Results are comparable with an edge in favor of Simpson's method. The second group of tests was more problematic as the integrands are not differentiable on the whole domain. Just like in the case of Rigorous Polynomial Approximations, when this was the case on a subinterval, we reverted to the naive method. For lack of space in the table below, we redefine here In both cases, the input depth has to be severely constrained compared to Section 12.3 in order to get a result in a reasonable time. The column "Accy" corresponds to the accuracy in bits which was asked for and is the same as with Taylor Models.

For I 1 , the computation took 33 seconds instead of 2.7, using the same (asked) accuracy of 2 -54 , using 50 points in the Simpson compound rule but only allowing a depth of 10 instead of 28 in the dichotomy. We can only conclude that in both cases, Rigorous Polynomial Approximations give more accurate results earlier in the binary search; thanks to this property, giving a high depth parameter like 28 or 35 does not penalize the search since this is only used around the singularities. However, Simpson's method doesn't seem so effective in reducing the breadth of search and we are forced to constrain the depth more, which in turns gives much less accurate results.

These results should be taken with caution as they are computed using an unproved implementation which could very well have flaws and inefficiencies. However, it seems to indicate that a quadrature method such as Simpson's method has the same order of magnitude as our methods based on Rigorous Polynomial Approximations. Our observations seem to be supported by Rall and Corliss [START_REF] Corliss | Adaptive, self-validating numerical quadrature[END_REF] whose FORTRAN implementation in the 1980s of a similar scheme was competitive with routines from the then state-of-the-art QUAD-PACK numerical integration library also written in FORTRAN.

Chapter 14

Conclusion and Perspectives

In this part, we have presented a method for computing and formally verifying numerical enclosures of univariate definite integrals inside the Coq proof assistant. This method deals with both proper and improper integrals. It has been integrated into the interval tactic of CoqInterval. In the same computation step, it proves that integrals exist and provides formally verified enclosures of their value. These proofs rely on the formal theory of generalized Riemann integrals provided by our extension to the Coquelicot library. In the proper case, the enclosing method only requires that there exist a Rigorous Polynomial Approximation of the integrand, so that its only limit is the underlying CoqInterval library. Any new functions added to the set of expressions E handled by the library would be almost immediately supported.

The current treatment of improper integrals is less automated. In particular, the syntactic expression of the integrand has to make explicit the element of the appropriate scale which models its asymptotic behavior near the singularity. The tactic currently supports two scales, the exponential e γx and x α ln β x. We could provide more scales to users, or at least merge these two into the more common e γx x α ln β x scale which has the advantage of being totally ordered with respect to the "small o" relation o(•) at e.g. +∞.

One possible direction to automate the discovery of asymptotics could then be to build yet another interpretation for elements in E in terms of their asymptotic behavior. Suppose for example that u > 0 and we define a function

• [u;+∞) : E → I×Z 3 with the specification that for all real expressions e(x) ∈ E, if e [u;+∞) = (x, (α, β, γ)) then there exists a continuous function f such that ∀x ∈ [u; +∞), f (x) ∈ x and e(x) R = f (x) • x α ln β (x) e γx .

(14.1)

We would need to explain how this interpretation propagates along the operators of E. For example, multiplication is straightforward: if e 1 [u;+∞) = (x 1 , (α 1 , β 1 , γ 1 )) and e 2 [u;+∞) = (x 2 , (α 2 , β 2 , γ 2 )), then e 1 • e 2 [u;+∞) could be defined as (x 1 • x 2 , (α 1 + α 2 , β 1 + β 2 , γ 1 + γ 2 )).

A more difficult example is addition, with two cases. In the easy case, if (α 1 , β 1 , γ 1 ) = (α 2 , β 2 , γ 2 ) then e 1 + e 2 [u;+∞) = (x 1 + x 2 , (α 1 , β 1 , γ 1 )) works. However, if x α1 ln β1 (x) e γ1x = o(x α2 ln β2 (x) e γ2x ) at +∞, then if f 1 and f 2 are the continuous functions from the specification (14.1) for e 1 and e 2 respectively, for x ∈ [u; ∞)]: e 1 (x) + e 2 (x) =f 1 (x) • x α1 ln β1 (x) e γ1x + f 2 (x) • x α 2 ln β2 (x) e γ2x = f 1 (x) • x α1-α2 ln β1-β2 (x) e (γ1-γ2)x + f 2 (x) • x α2 ln β2 (x) e γ2x . Since x α1-α2 ln β1-β2 (x) e (γ1-γ2)x is continuous and bounded, we would then be able to compute a new interval y such that e 1 (x) + e 2 (x) [u;+∞) = (y, (α 2 , β 2 , γ 2 )) while being correct with respect to the specification in Formula (14.1). This is possible because the scale constituted by the x α ln β (x) e γx is totally ordered with respect with the relation o(•).

After we would program a procedure which computes ∞ u x α ln β (x) e γx dx given α, β and γ, we would then treat previously inaccessible examples at the current stage of the code as seen in Section 12.2.4.

In the introduction of this part, we opposed specification to validation of computations. However, there need not be such a clear cut between the two notions, and parts of our computations in the context of integrals could be delegated to an external oracle. For example, the "depth tree" describing where to use dichotomy in an interval when computing a proper integral could be provided by a fast Ocaml program. Similarly, the initial estimate used to find the absolute error from the relative error provided by the user in the integral tactic need not be proved and we can replace it by the approximation part (without the error) of Simpson's method.

Before writing the programs and tactics described in Section 12, we wrote prototypes in Ocaml, as well as other routines for solving Ordinary Differential Equations using Taylor Models. It might seem at first that in the absence of a need for complex built-in features, prototyping in Ocaml before writing a Coq program would be a duplication of effort. However, certain realities of Coq development made this orders of magnitude faster in our case.

First, Coq is not such a great testing environment: computations are slow, there is no profiling, logging or printing possible1 and in particular no prettyprinting, for example of floats. These aspects of programming are crucial for rapid development.

Secondly, since the libraries we could have used for prototyping are geared towards their use in tactics, most of what they do is handled by Ltac which is untyped, unhelpful when it fails, and needs a lot of side-knowledge to be made to work. Alternatives like Mtac [START_REF] Ziliani | Mtac: A monad for typed tactic programming in coq[END_REF] may be steps in the right direction.

Nested integrals are not supported by our method. The naive enclosure approach could easily be adapted to support them, but performances would be even worse due to the curse of dimensionality. As there exists no general approach for integrating multivariate polynomials,2 being able to compute rigorous multivariate polynomial approximations would presumably not help.

While our adaptive bisection algorithm and our rigorous quadrature based on primitives of polynomials might seem crude, they proved effective in practice. They produce accurate approximations of non-pathological integrals in a few seconds, and thus they are usable in an interactive setting. Moreover, they are able to handle functions with unbounded second derivatives in a rigorous way, as well as unbounded integration domains. Another contribution is the way we are able to infer that a function is integrable from a successful computation of its integral.

Following Rall & Corliss [START_REF] Corliss | Adaptive, self-validating numerical quadrature[END_REF], we could have extended this adaptivity even more to include the order chosen for each subinterval.

We could also have tried a much more general method, that is, solving a differential equation built from the integrand, as we did with VNODE-LP. There has been some work done for Coq in the setting of exact real arithmetic [START_REF] Makarov | The Picard algorithm for ordinary differential equations in Coq[END_REF] [99], but the performances are not good enough in practice. Much closer to actual numerical methods is Immler's work in Isabelle/HOL [START_REF] Immler | Formally verified computation of enclosures of solutions of ordinary differential equations[END_REF], which uses an arithmetic on affine forms. This approach is akin to computing with degree-1 RPAs. Note that such computations in Isabelle are treated in a slightly different way from Coq: they are performed by an extracted Ocaml program whose result is then imported inside the proof assistant.

Part IV Conclusion 127

In this concluding part, we discuss the general question of when and how to trust proofs relying on computations. Then, we suggest how the domains of formal methods and cryptocurrencies might benefit each other.

Can we trust computations used in mathematical proofs?

Despite numerical integration being an old problem and despite the relative simplicity of the mathematical ingredients of our methods, we easily found mistakes on elementary examples. One cannot exclude other mistakes in numerical computations in published works. We can speculate that most of the time, as in the case of the slightly off integral we found in the Ternary Goldbach Conjecture, they are only slightly wrong and do not impact the correctness of results. However, not only is there no guarantee that this is always the case, but this level of uncertainty is not acceptable.

In both cases, the mistakes we found happened despite using state-of-the-art integration software. However, computation-based proofs which use customwritten code are not a rare occurrence. In that case, there are even further reasons to be circumspect. We closely examined the proof of correctness of the C++ code in the original proof of the Double Bubble conjecture [START_REF] Hass | Double bubbles minimize[END_REF], and it turns out that the pseudo-code described in the paper has discrepancies with the distributed code. For example in Figure 14.1, the value assigned to the variable R is We re-wrote the pseudo-code in Ocaml using our own interval arithmetic library and barring an error on our side, it is not semantically equivalent. However, when we changed the parameters scrupulously to the actual C++ code we were able to replicate their results. Thus, it is difficult to assess whether these were harmless changes with a justification that did not make it to the paper or if there were typos in the code, rendering the proof invalid. This seems to validate a posteriori the decision of the reviewers of the first proof of the Kepler conjecture which relied on 3 Gigabytes of code [START_REF] Thomas | The Kepler conjecture[END_REF] to assert that they could not certify the correctness of the proof: The Double Bubble code is only 849 lines or 21 kilobytes of code.

We acknowledge that there is no satisfactory tool, as of now, which will meet the computational needs of all modern mathematical proofs at the level of assurance of formal proofs. There are already a range of certified tools which tackle specific domains, such as solving differential equations numerically [START_REF] Immler | Formally verified computation of enclosures of solutions of ordinary differential equations[END_REF], proving nonlinear inequalities [START_REF] Victor Magron | Nlcertify: A tool for formal nonlinear optimization[END_REF], checking satisfiability [START_REF] Ekici | Extending smtcoq, a certified checker for smt[END_REF], or checking primality [START_REF] Grégoire | A Computational Approach to Pocklington Certificates in Type Theory[END_REF]. It is possible that some types of computations will remain out of reach of proof assistants in the near future. For these, one may either develop a new tool to produce formally verified results or fall back on unverified computations.

We contributed to this Swiss-knife-like array of tools destined to the verification of proofs relying on demanding computations, by providing a methodology to validate the recurrences outputted by creative telescoping algorithms. However, in this case, this validation is not as automatic as we had hoped it to be. This work is an example of another problem: the program may be doing exactly what it is supposed to do, but do the semantics of the object it produces match with those of the object we think we are reasoning on? For this reason, using a proof assistant is of precious help.

Some Remarks Concerning Tools for Formal Proofs

In Parts II and III, we commented on the ease or difficulty of using a proof assistant. If it may seem at times that our comments are too negative, this is also because of our enthusiasm and eagerness for improvements to theorem proving. One should also keep in mind that our experience was mainly with the Coq proof assistant and in particular with the Mathematical Components libraries, using the SSReflect proof language; surely other proof assistants provide a better experience on some aspects, and worse on others.

During an internship in the summer of 2016, we used the Lean theorem prover [START_REF] De Moura | The Lean Theorem Prover (System Description)[END_REF] to attempt to port some of Mathematical Components's elementary group theory. Lean is based on the Calculus of Inductive Constructions, like Coq, but it is comparatively very young. For this reason, at the time of this internship, the tactic language was still rudimentary and few libraries were available: this made the task very hard, and we did not succeed in proving Sylow's theorem as we had intended in the beginning. This is of course no fault of Lean's, which has been evolving at a remarkable speed, but it is a testimony to how much combined work has been put into making theorem proving in Coq both effective and agreeable.

Formal Proofs and Cryptocurrencies

We would like to discuss how the domain of computer-aided proofs, with its tools and principles can benefit the emerging domain of blockchains and cryptocurrencies, but also how this domain raises new challenges.

In the domain of formal proofs, as ruthless as the peer review process may sometimes be, authors and reviewers are in a relationship of trust. We mean this in the sense that even though the reviewer carefully checks all of the author's claims, she does not assume that the person who wrote the proof is actively trying to trick them into thinking that they have proved a certain theorem, when in fact they have either proved a different theorem or made use of a bug to build a fake proof.

In 2014, on a website called Proof Market which set bounties in the digital currency Bitcoin for Coq proofs of theorems submitted by users, someone managed to walk away with 1 Bitcoin (about 600 euros at the time) by providing a proof of False along the lines of the following script: This website created a completely new setting for formal proofs, where an automated checker would trigger a transaction upon receiving (what it thought to be) a valid proof. Since no human was in the loop, this created an incentive to find bugs in the checker rather than to write a correct proof.

A Blockchain is a linked list of blocks such that the i + 1-th block contains a cryptographic hash of the i-th block. Blockchains are typically built using a distributed protocol which gives the right to add a new block to a node selected randomly in proportion of their commitment in the form of some rare resource (such as computing power, or money, or storage, etc...). This scheme was invented by Satoshi Nakamoto who created the first decentralized digital currency, Bitcoin, where blocks contain lists of transactions [START_REF] Nakamoto | Bitcoin: A peer-to-peer electronic cash system[END_REF]. Several improvements have been proposed since, including the idea of a general purpose cryptocurrency where transactions contain function calls and blocks additionally maintain the state of a global decentralized computer [START_REF] Wood | Ethereum: A secure decentralised generalised transaction ledger[END_REF]. In particular, the notion of "smart contracts", which are autonomous program-agents imagined in 1995 by Nick Szabo, became a reality in the shape of autonomous pieces of code "running" on the blockchain. In light of several high-profile hacks and thefts of tens of millions of euros in recent years, providing formal guarantees on the behavior of such code, which can typically manipulate considerable amounts of money or handle authorization of access to other contracts, has become an urgent problem.

We have participated in a Microsoft Research hackathon in which we built a proof-of-concept program capable of analyzing both high-level and low-level code of smart contracts in the Ethereum cryptocurrency [START_REF] Bhargavan | Formal Verification of Smart Contracts: Short Paper[END_REF]. We used the F * [START_REF] Swamy | Dependent types and multi-monadic effects in F*[END_REF] functional language.

We also discovered the Easycrypt [START_REF] Barthe | Easycrypt: A tutorial[END_REF] proof assistant, geared towards proofs on cryptographic protocols and more generally on probabilistic programs. We formalized a model called "the Bitcoin Backbone Protocol" [START_REF] Garay | The bitcoin backbone protocol: Analysis and applications[END_REF], an attempt at proving elementary properties of the Bitcoin protocol, together with Pierre-Yves Strub. However, a complete formal proof is likely to necessitate improvements to Easycrypt.

Finally, we did some work in formalizing in Coq the semantics of the smartcontract language of another general-purpose cryptocurrency called Tezos [START_REF] Goodman | Tezos: A self-amending crypto ledger position paper[END_REF]. The main difference with Ethereum is that Tezos proposes formal semantics of this language, which, like Ocaml, is strongly typed and purely functional.

To sum up, we think this is a new type of context for formal proofs, because of the unique adversarial model. When one checks a formal development, one does not usually assume malice from the person who wrote it. By contrast, if someone anonymous manages to take advantage of a weakness in a theorem prover to pretend that a smart contract has a certain behavior, they could walk away with millions with no chance of ever getting caught. Specific tools to tackle such security challenges are being developed.

Figure 1 . 1 :

 11 Figure 1.1: A Coq session stating: ∀x ∈ N, x = x.

  in b * b * (if n mod 2 = 0 then 1 else a);; let fermat n = pow 2 (pow 2 n) + 1;; is_prime (fermat 5);; Listing 2.1: A simple program which finds a factor of F 5 It returns # not prime: 641 is a factor.

Figure 3 . 1 : 8 Conjecture 1 (

 3181 Figure 3.1: Example of a tile with n = 8

Figure 3 . 2 :

 32 Figure 3.2: Tiling solutions for n = 4 and n = 5

Figure 3 . 3 and

 33 Figure 3.4 show solutions for n = 6 and n = 7. We also show symmetric solutions found by Paul Zimmermann and his intern Romain Lemoine in Figure3.5 and Figure3.6. This brings a possible variation of the conjecture: for odd n, there exists a solution which is symmetric with respect to the central vertical line of the grid.In the case of n = 7, the computation takes more than one hour on a standard laptop. Interestingly, when one forces the position of the two lateral bars (bottom left or bottom right) of size n -1, the computation only takes 20 seconds.

Figure 3 . 3 : 6 Figure 3 . 4 : 7 Figure 3 . 5 :

 33634735 Figure 3.3: A tiling solution for n = 6

Figure 3 . 6 :

 36 Figure 3.6: A symmetric solution for n = 7

  K) be matrices over a field K, of respective dimensions m × n and n × p. The matrix C = (c i,j ) k b k,j is called the product A • B of A and B.

4. 4

 4 Scaling the Experiment Up: Implementations from Maple to Ocaml 4.4.1 Maple Implementation: Description, Advantages and Downsides In Section 4.3, we considered a simple question: what do algorithms "produced" with the τ -theorem look like in practice?

  √ 163) ∈ N [104]. The Maple session in Listing 5.1 shows why one might believe this: using 29 digits for intermediary computations in the Maple software, it seems as though exp(π √ 163) = 262537412640768744: > Digits := 29: > evalf(exp(Pi * sqrt(163))); 262537412640768744.00000000000 > Digits := 40: > evalf(exp(Pi * sqrt(163))); 262537412640768743.9999999999992500725972 Listing 5.1: A Maple session examining Ramanujan's "almost integer"On the more trusting side, consider the game of Sudoku which consists in filling out a 9 by 9 grid with integers so that for each row, column, or 3 by 3 square box as delimited in bold in Figure5.1, there be exactly one instance of each integer i ∈ {1,• • • , 9}. If a piece of software fills out the empty squares in a Sudoku grid with numbers, one need only do a quick visual check in order to be convinced that they are correct.

Figure 5 . 1 :

 51 Figure 5.1: A Sudoku grid

  ),ζ(7),ζ(9), and ζ(11) must be irrational, ζ(

  i . (7.1) and the Coq counterpart is Lemma fact_logp_sum_small p j n : prime p -> (n < p ^j.+1) -> logn p n'! = \sum_(i < j) (n %/ p ^i.+1). Listing 7.1: The Coq statement of Equation (7.1).

Proposition 1 . 1 m 3

 113 The sequence z n = n m=1 is a Cauchy real. Lemma creal_z3seq : creal_axiom z3seq. Listing 7.2: Coq counterpart of Proposition 1. The predicate creal_axiom expresses the property of being a Cauchy real.

Chapter 8 Formalizing 1 m 3

 813 Apéry's theorem in Coq 8.1 Recurrence and Creative Telescoping Recall from Equation (6.1) that (z n ) n∈N is the sequence n m=1 , and that

Figure 8 . 1 :

 81 Figure 8.1: Proving that V is C × U

(

  * Coefficients of every recurrence, P, and Q. *) Definition Sn_cf0_0 n k := (n + 1 + k) 2 / (-n -1 + k) 2 . Definition Sk_cf0_0 n k := (-n + k) 2 * (n + 1 + k) 2 / (k + 1)4 . Definition P_cf0 n := (n + 1) 3 . ... (* Conditional recurrences. *) Definition Sn c := ∀ n k, precond.Sn n k -> c (n + 1) k = Sn_cf0_0 n k * c n k Definition Sk c := ∀ n k, precond.Sk n k -> c n (k + 1) = Sk_cf0_0 n k * c n k (* Operators P and Q. *) Definition P c n := P_cf0 n * c n + P_cf1 n * c (n + 1) + P_cf2 n * c (n + 2). ... (* Statement P = ∆ k Q. *) Definition CT_premise c := ∀ n k, precond.CT_premise n k -> P (c ^~k) n = Q c n (k + 1) -Q c n k.

  Module precond. Definition Sn n k:= (k != n + 1) ∧ (n != -1). Definition Sk n k := (k + 1 != 0) ∧ (n != 0). Definition CT_premise n k := (n >= 0) ∧ (k >= 0) ∧ (k < n). End precond.Listing 8.3: Hand-written provisos in annotated rec c Listing 8.2, we also have at our disposal a list representation [:: P_cf0; P_cf1] of the same operator.

============================ 0 Listing 8 . 7 :

 087 P_cf2 n * c (n + 2) n + P_cf1 n * c (n + 1) n + P_cf0 n * c n n + Q_cf0_0 n n * c n n + P_cf1 n * c (n + 1) (n + 1) + P_cf2 n * c (n + 2) (n + 1) + P_cf2 n * c (n + 2) (n + 2) = Terms around the upper bound

0 Listing 8 . 9 :

 089 rewrite (ann_c.Sk (n + 2) (n + 1)); last by intlia. rewrite (ann_c.Sk (n + 2) n); last by intlia. rewrite (ann_c.Sk (n + 1) n); last by intlia. rewrite (ann_c.Sn (n + 1) n); last by intlia. rewrite (ann_c.Sn n n); last by intlia. set cnn := c n n. Fail set no_more_c := c _ _. Listing 8.8: Reduction modulo the Gröbner basis of c ============================ P_cf2 n * Sn_cf0_0 (n + 1) n * Sn_cf0_0 n n * cnn + P_cf1 n * Sn_cf0_0 n n * cnn + P_cf0 n * cnn + Q_cf0_0 n n * cnn + P_cf1 n * Sk_cf0_0 (n + 1) n * Sn_cf0_0 n n * cnn + P_cf2 n * Sk_cf0_0 (n + 2) n * Sn_cf0_0 (n + 1) n * Sn_cf0_0 n n * cnn + P_cf2 n * Sk_cf0_0 (n + 2) (n + 1) * Sk_cf0_0 (n + 2) n * Sn_cf0_0 (n + 1) n * Sn_cf0_0 n n * cnn = Rational function with folded coefficients

Figure 8 .Figure 8

 88 Figure 8.2 describes the global dependencies of the files proving all the steps in Table8.1. In order to complete the formal proof of Lemma 13, we verify

  2) using four evaluations b 0 , b 1 , b 2 , b 3 that we compute in file initial_conds. The reasoning is the following: define b as the unique sequence obtained by fixing b 0 = b 0 , b 1 = b 1 and such that b satisfies the Apéry recurrence (6.2) of order 2. Then prove that b coincides with b for n = 0, 1, 2, 3. Finally, establish that any solution of Apéry's recurrence is also a solution of R 4 . Thus b = b, and this proves that b satisfies Apéry's recurrence.

Corollary 1 .

 1 For any n ∈ N, 0 < ζ(3) -bn an . Proof. Since ζ(3) = b a , it is sufficient to show that for any k < l, we have 0 < b l a l -b k a k . Thus it is sufficient to observe that for any k, we have 0 < b k+1 a k+1 -b k a k , which follows from Lemma 16.

Corollary 2 .

 2 ζ(3) -bn an = O( 1 a 2 n

Remark 3 .

 3 the bound above with k = 4 we get w < w 4 • α For k ≥ 1, we have

[ 3 .

 3 14159; 3.14160].

11. 2 . 1

 21 CoqApprox: Rigorous Polynomial Approximations and Taylor Models in Coq

Figure 12 . 1 :

 121 Figure 12.1: Sollya Session > TL=taylorform(sin(x + exp(x)), 10, 0, [0,8]); > TL[2]; [-4.2589929933332345326623361369068740662104390570317e30; 4.2458685067031149053103989041910036815095504573458e30]

  Axiom completeness : ∀ E : R -> Prop, bound E -> (∃ x : R, E x) -> {m : R | is_lub E m} where is_lub is defined in the following way: Definition is_upper_bound (E:R -> Prop) (m:R) := ∀ x:R, E x -> x <= m. Definition is_lub (E:R -> Prop) (m:R) := is_upper_bound E m ∧ (∀ b:R, is_upper_bound E b -> m <= b).

∞ 0 e -x 2 dx = √ π 2

 02 Rbar_locally x:is_RInt_gen (fun x => exp(-x 2 )) (at_point 0) (Rbar_locally p_infty) but it is not as straightforward to express that ∞ 0 e -x 2 dx ∈ [0; 1]: ∃ l, is_RInt_gen (fun x => exp(-x 2 )) (at_point 0) (Rbar_locally p_infty) l ∧ l ∈ [0; 1].

Definition

  RInt_gen (f : R -> V) (a b : (R -> Prop) -> Prop) := iota (is_RInt_gen f a b).

  Lemma is_RInt_derive (f df : R -> V) (a b : R) : (∀ x : R, Rmin a b <= x <= Rmax a b -> is_derive f x (df x)) -> (∀ x : R, Rmin a b <= x <= Rmax a b -> continuous df x) -> is_RInt df a b (minus (f b) (f a)). which states that if f is continuously differentiable on [a, b], and if df : R → R coincides with f on [a, b], then b a df (x)dx = f (b) -f (a), (12.15) had been proved for f : R → R using the Intermediate Value Theorem and Riemann sums. Such a proof does not generalize to normed modules in general.A more usual way of proving this is to setF (x) := f (a) + x a df (x) dx -f (x) (12.16)and then to observe that ∀x ∈ [a, b], F (x) = 0 (12.17

  )

  Lemma eq_is_derive : ∀ (f : R -> V) (a b : R), (∀ t, a <= t <= b -> is_derive f t zero) -> a < b -> f a = f b.

  Thus, x ∈ A so that [c; c + δ] ⊂ A. We know that u ≤ b. Moreover, u cannot be less than b by the above argument. The only remaining option is that b = u. Hence, for all x ∈ [a; b], |f (x) -f (a)| ≤ ε. Since this is true for arbitrary ε, we conclude that ∀x ∈ [a; b], f (x) = f (a).

  x a f (x) at x ∈ [a; b] is f (x), which in Coquelicot is formulated more generally as: Lemma is_derive_RInt (f If : R -> V) (a b : R) (Hloc : locally b (fun x => is_RInt f a x (If x))) : -> continuous f b -> is_derive If b (f b).

  x a g (x)dx and g . Hence is_derive_RInt tells us that ∀x ∈ [a; b],x a g (x)dx = g (x)(12.19)

  Lemma is_RInt_comp (f : R -> V) (g dg : R -> R) (a b : R) : (∀ x, Rmin a b <= x <= Rmax a b -> continuous f (g x)) -> (∀ x, Rmin a b <= x <= Rmax a b -> is_derive g x (dg x) ∧ continuous dg x) -> is_RInt (fun y => scal (dg y) (f (g y))) a b (RInt f (g a) (g b)).

1 0x 4 +

 14 10x 3 + 19x

I 1 := 1 0 x 4 + 2 := 1 - 1 2048x 12 - 1 exp -x -3 4 2 1 -

 0421112141 10x 3 + 19x 2 -6x -6 e x dx and I 6144x 10 + 6912x 8 -3584x 6 + 840x 4 -72x 2 + x 2 dx.

122For I 2 ,

 2 we reduced the depth from 35 to 10, only computing Simpson with 10 points and still the computation takes 204 seconds, and returns a result with an accuracy of 10 -4 instead of 10 -15 .Int. Time Time (TM) Error Actual Accy Depth Prec. Points I

√ 3 2 -Hi- 1 Y in the C++ code and 1 -Hi- 1 Y

 2111 in the pseudocode.

Figure 14 . 1 :

 141 Figure 14.1: Discrepancy between the pseudo-code in the paper and the distributed C++ code for the Double Bubble Conjecture

  (-(1/p))*bb[0,ebb],ebb=1..p-1))*epsilon3 (1/p)*v[ev2,i],ev2=1..p-1))*epsilon 4 ; p*w[1,i] :=0,0] = -sum(a[ha,0],ha=1..m-1) (* constr_a_1 *) w[0,i:=0] = -sum(w[0,hw0],hw0=1..m-1) (* constr_w_1 *) w[1,i:=0] = -sum(w[1,hw1],hw1=1..m-1) (* constr_w_2 *) y[i:=0,0] = -sum(y[hy,0],hy=1..m-1) (* constr_y_1 *) yy[i:=0,0] = -sum(yy[hyy,0],hyy=1..m-1) (* constr_y_2 *) u[0,k:=0] = -sum(u[0,eu0],eu0=1..p-1) (* constr_u_1 *) u[1,k:=0] = -sum(u[1,eu1],eu1=1..p-1) (* constr_u_2 *) c[i,k:=0] = -sum(c[i,ec],ec=1..p-1) (* constr_c_1 *) cc[i,k:=0] = -sum(cc[i,ecc],ecc=1..p-1) (* constr_c_2 *) x[k:=0,i] = -sum(x[ex,i], ex=1..p-1) (* constr_x_1 *) xx[k:=0,i] = -sum(xx[exx,i],exx= 1..p-1) (* constr_x_2 *) c[i:=0,k] = 0 cc[i:=0,k] = 0 v[k,i:=0] = 0 x[k,i:=0] = 0 xx[k,i:=0] = 0 b[0,k:=0] = 0

r k . For each µ in the direct sum on the right side of Equation (4.14), let us denote by E µ the vector space in which each term lives. Now expand t ⊗k 0 as in Equa- tion (4.14) and choose

  µ ∈ N s with µ 1 +• • • + µ s = k.If we zero out all variables which do not belong to E µ , we get:

						.14)
	Similarly to Equation (4.6), the ∼ = sign hides a bijective rewriting function
	to go from left to right.
	s
	Denote t 0 =	m i , n i , p i , recall that we have by hypothesis:
	i=1
			R(t 0 ) ≤ r
	so that
	R(t ⊗k 0 ) ≤ R s s k i=1 i=1 µ m µi i , n µi i ,	s i=1	p µi i	≤ r k .	(4.15)
	We rename the quantities in Equation (4.15) as :
		R (S k	M k , N k , P k ) ≤ r k .	(4.16)
	At this point, we prove the following
	Lemma 9. If mnp > 1 and s < r,
	R (s	m, n, p ) ≤ r ⇔ s(mnp) ω ≤ r ⇒ ω ≤	log(r/s) log(mnp)	.

  2 -6x -6 e x dx 11.14731055005714On this example, quadrature methods have some troubles: quad gives only 10 correct digits; verifyquad gives a false answer (a tight interval not containing the value of the integral) without warning; 2 quadgk gives only 9 correct digits. VNODE-LP cannot be used because of the absolute value.

	Error Time Accuracy Degree Depth Precision
	10 -3	0.7	14	5	8	30
	10 -6	0.9	24	6	13	40
	10 -9	1.3	34	8	18	50
	10 -12	1.9	44	10	22	60
	10 -15	2.7	54	12	28	70

  25/τ 2 )+4.1396+ln π

				2
		ln τ 1 + 0.25/τ 2	•	ln 2 τ τ 2 dτ 3.17742 • 10 -3 .
	Error Time Accuracy Degree Depth Precision
	10 -3	0.6	2	3	0	30
	10 -4	0.8	5	5	2	30
	10 -5	1.5	8	7	6	30
	10 -6	3.1	11	9	11	30
	10 -7	5.6	15	12	12	30
	10 -8	11.2	18	15	15	30

Symbolic Computation National Days, a yearly meeting of the French community of Computer Algebra and Symbolic Computation.

We have not tested it on a Mac, but there shouldn't be any specific counter-indication.

https://github.com/exPHAT/twitter-sort

This prelude can be seens at https://coq.inria.fr/library/Coq.Init.Prelude.html.

n ∈ R + if and only if x ∈ R + .

A proof of this error bound along with ones for more advanced quadrature methods are detailed in the book Methods of Numerical Integration[START_REF] Davis | Chapter 4 -error analysis[END_REF] byDavis and Rabinowitz.

The only thing that will later distinguish the integration variable t from an actual constant such as π is that its value is placed at the top of the initial evaluation stack.

where [x 0 , x 1 ] represents the interval whose ends are x 0 and x 1 irrespective of their relative order. This interval is included in X 0 .

https://coq.inria.fr/distrib/current/stdlib/

The bug lies in an incorrect implementation of Taylor models for absolute value. It has now been fixed by removing support for absolute values.

although this is changing with a recent pull request https://github.com/coq/coq/pull/

Any

[START_REF] Ahmed | Ahmed's integral: the maiden solution[END_REF]-SAT instance can be reduced to approximating the integral of a multivariate polynomial.

Remerciements

Part II

Certifying Computations

In the last part of the file, see Listing 8.4, a record collects the elements of the Gröbner basis C. Maple indeed often produces a larger set of annihilators for a given sequence, for instance CT_premise in Listing 8.2 is related to a creative telescoping pair but not to the Gröbner basis. Also, the Gröbner basis can be obtained by refining a first set of annihilators, which happens at step 4 of Table 8.1.

Formal proofs of a conditional recurrence

We take as a running example the file ops_for_a, which models step 2 in Table 8.1. This file proves theorems about an arbitrary sequence c satisfying the recurrences in the Gröbner basis displayed on Listing 8.4, and about the sequence a by definite summation over c.

Require Import annotated_recs_c. Variables (c : int -> int -> rat) (ann_c : Annihilators c).

Theorem P_eq_Delta_k_Q : CT_premise c. Proof. ... Qed.

Let a n := \sum_(0 <= k < n + 1) c n k.

The formal proof of lemma P_eq_Delta_k_Q is an instance of Equation (8.12). Using this property, we prove that the sequence a verifies a conditional recurrence associated to the operator P . As suggested in Section 8.1.3, this proof consists in applying the lemma sound_telescoping (see Equation (8.13) and Listing 8.1), which formalizes a sound creative telescoping, and in normalizing to zero the resulting right-hand side of Equation (8.13). Listing 8.5 displays the first lines of the corresponding proof script, which select and name the three components of the right-hand side of Equation (8.13), with self-explanatory names. The resulting proof context is displayed on Listing 8.6.

Theorem recApery_a n (nge2 : n >= 2) : P a n = 0. Proof. rewrite (punk.sound_telescoping P_eq_Delta_k_Q). set boundary_part := (X in X + _ + _). set singular_part := (X in _ + X + _). set overhead_part := (X in _ + _ + X). In this section, we describe how to compute a numerical enclosure of the real number v u f from enclosures of the finite bounds u and v and of the Riemannintegrable function f . In other words, we only deal with proper integrals: the case of improper integrals will be handled later in Section 12.2. We describe two basic methods based respectively on the evaluation of a simple interval extension and on a polynomial approximation of f . They can be combined and improved by a dichotomy process. This is integrated into the interval tactic of CoqInterval.

Listing 8.5: Begining of a proof of sound creative telescoping

A Naive Approach

Our first approach uses an interval extension of the integrand. Let us formalize and extend this notion already mentioned in Section 11.1.2.

Definition 15. For any function f : R n → R, a function F :

In the rest of the section we suppose that F : I → I is an interval extension of the univariate function f , and we want to compute an enclosure of

, hence Formula (12.1). The case v ≤ u is symmetrical.

93

In practice we do not compute with f but only with its interval extension F. Moreover, we want the computations to operate using only enclosures of the bounds. So we adapt Formula (12.1) accordingly.

Lemma 26 (Interval naive integral enclosure). For any intervals u, v such that u ∈ u and v ∈ v, we have

(12.2)

Note that if u and v are point intervals and if F is the optimal interval extension of f , then equation (12.2) reduces to equation (12.1).

Proof. If u ∈ u and v ∈ v, then by (12.1) and reusing notations from the proof, we have

, by definition of the closed convex hull.

The naive_integral Coq function implements (12.2). Given u, v ∈ I and F a function of type I → I, (naive_integral prec F u v) computes an interval i using floating-point arithmetic at precision prec.

Another Approach using RPAs

The enclosure method described in Section 12.1.1 is crude. Better knowledge of the integrated function allows for a more efficient approach. The goal of this section is to use RPAs introduced in Section 11.2.1 to obtain better enclosures of integrals.

Lemma 27 (Polynomial approximation). Suppose f is approximated on [u; v] by p ∈ R[X] and ∆ ∈ I in the sense that ∀x ∈ [u; v], f (x) -p(x) ∈ ∆. Then for any primitive P of p we have

Note that our method and proofs do not depend on the way polynomial approximations are obtained.

Incidentally, we recently realized while reading Rall and Corliss [START_REF] Corliss | Adaptive, self-validating numerical quadrature[END_REF] that by putting c = u+v 2 and choosing p to be of the shape n i=0 a i (x -c) i , P can be chosen so that P (v) -P (u) can be computed in half as many operations. We plan to take benefit from this remark in our code in the near future.

Dichotomy

Both methods presented in Sections 12.1.1 and 12.1.2 can compute an interval enclosing v u f when u and v are proper bounds. As seen in Section 12.1.3, polynomial approximations usually give tighter enclosures of the integral, but not always, so we combine both methods by taking the intersection of their result.

This may still not be sufficient for getting a tight enough enclosure, in which case we recursively split the integration domain in two parts, adding up the results. The function integral_float_absolute performs this dichotomy and the integration on each subdomain. It takes an absolute error parameter ε; it stops splitting as soon as the width of the computed integral enclosure is smaller than ε. The function also takes a depth parameter, which means that the initial domain is split into at most 2 depth+1 subdomains. Note that, because the depth is bounded, there is no guarantee that the target width will be reached.

Let us detail more precisely how the function behaves. It starts by splitting

Otherwise, several cases can occur:

, the first enclosure is sufficient but the second is not. So integral_float_absolute calls itself recursively on [m; v] with depth -1 as the new maximal depth and ε -w(i 1 ) as the new target accuracy, yielding i 2 . The function then returns i 1 + i 2 .

• If w(i 1 ) > ε 2 and w(i 2 ) ≤ ε 2 , we proceed symmetrically. • Otherwise, the function calls itself on [u; m] and [m; v] with depth -1 as the new maximal depth and ε 2 as the new target accuracy, yielding i 1 and i 2 . It then returns i 1 + i 2 .

Accuracy

Both methods described in Sections 12.1.1 and 12.1.2 use a single approximation of the integrand on the integration interval. A decomposition of this interval into smaller pieces may increase the accuracy of the enclosure, if tighter approximations are obtained on each subinterval. In this section we give an intuition of how the naive and polynomial approaches compare, from a time complexity point of view. What we mean here by complexity is the number of operations needed to obtain a result of width less than a parameter 0 < ε given as an input to the algorithm.

The naive (resp. polynomial) approach here consists in using a simple interval approximation (resp. a valid polynomial approximation) to estimate the integral on each subinterval. Let us suppose that we split the initial integration interval, before computing integral enclosures: domain. More precisely, we can prove the following property:

By combining Formulas (11.4) and (12.11), we obtain a numerical method to prove that a function is continuous on a domain. Indeed, we just have to compute an enclosure of the function on that domain, and to check that it is not ⊥ I . A closer look at the way naive integral enclosures are computed provides the following corollary: whenever the enclosure of the integral is not ⊥ I , the function is actually continuous and thus integrable on any compact set of the input domain.

Note that Property (12.11) intrinsically depends on the operations that can appear inside p, i.e. the operations belonging to the class E of Section 11.1.2. Therefore, its proof has to be extended as soon as a new operator is supported in E. This proof relies on the invariant: Lemma 28 (Invariant on operators of E). For all unary operators u included in the definition of E, for all x ∈ R:

For all binary operators b included in the definition of E, for all x, y ∈ R:

Suppose that we want to extend E with the integer part • . Now Lemma 28 is not true any more. An unsolved question for now is: what more general invariant could be used to obtain integrability?

Estimating Improper Integrals

Following usual mathematical terminology, we qualify integrals as improper when their range is unbounded, as for

x 2 dx, or when the integrand is unbounded, such as 1 0 ln x dx. They are defined as limits of proper integrals. In this section, we describe how to compute a numerical enclosure of some improper integrals. Improper integrals are computed by splitting the interval into two parts, a proper part which is treated with the methods described previously in Section 12.1, and the remainder which is handled in a specific way. The splitting is automatically performed by a variant of the adaptive method presented in Section 12.1.4 where the splitting point m for [u; +∞) when u > 0 is chosen to be 2 • u. We now describe how we handle the remainder. On paper, when we want to show that an improper integral

(12.12)

where C > 0 and g is a well-known integrable function such as 1 x 2 , ln x or exp(-x 2 ). On top of establishing that the integral exists, this leads to a natural way of bounding the integral:

Establishing the bound (12.12) is not so easy to automate as it can involve complicated algebraic manipulations of usual functions. For this reason, we restrict ourselves to the case when this work has already been done for us. We consider improper integrals of the shape v u f g where either u = 0 + or v = +∞, and f is bounded. Function g belongs to a catalog of functions with known enclosures of their integral, such as x α ln β x.

In Sections 12.2.1 and 12.2.2 we focus on integrals of the shape +∞ u f g. Section 12.2.1 presents the general theorem, while Section 12.2.2 lists the functions contained in our catalog. Finally, Section 12.2.3 focuses on integrals of the shape v 0 + f g.

Improper integral of a product

Let h : R → R be a function. To determine that +∞ u h exists, we have added a proof of the following Cauchy criterion: this integral exists if and only if for any v ≥ u, v u h exists and for all ε > 0, there exists M > 0 such that for all u, v ≥ M , | v u h| ≤ ε. We use this criterion to show the following lemma. Lemma 29. Let f, g : R → R. Suppose that, on [u; +∞), f is bounded, f and g are continuous, and g has a constant sign. Moreover, suppose +∞ u g exists. Then +∞ u f g exists, and

We provide an effective version of the previous lemma, in the same spirit as Lemma 26, with a similar proof: Lemma 30. Let F, I g : I → I be interval extensions respectively of f and x → +∞ x g. For any interval u such that u ∈ u:

Definition Rbar_locally (a : Rbar) (P : R -> Prop) := match a with | Finite a => locally a P | p_infty => ∃ M : R, ∀ x, M < x -> P x | m_infty => ∃ M : R, ∀ x, x < M -> P x end.

The predicate locally a P expresses that the predicate P is valid in a neighborhood of the real number a.

Given a point x ∈ R, at_point x is a filter for x:

Definition at_point (P : R -> Prop) : Prop := P x.

Coquelicot also features a formalization of the product of two filters:

Lemma 32 (product of two filters). Let F be a filter of the set S and G be a filter of the set

Coquelicot adopts a systematic approach of replacing dependent types used in the standard library (like in Listing 12.2.5) with total functions. It also provides an alternative formalization of the theory of Riemann integration. Let V be a complete normed R-module: When the function f is Riemann-integrable on [u; v], the value

Otherwise it is left unspecified and thus most properties about the actual value of

The library also features a notion of generalized Riemann integrals used to represent integrals such as +∞ 0 + ln x/(1 + x 2 )dx. The bounds of such integrals are not real numbers, but filters. There is a predicate is_RInt_gen similar to is_RInt:

The definition of is_RInt_gen amounts to the following. Let f : R → V be a function, and let u, v be two bounds. If l ∈ R, l is said to be the generalized integral of f between u and v if b a f (t)dt is well-defined for all (a, b) in a neighborhood of (u, v) and if

Contrary to the previous case of proper integrals, the library does not feature a total function computing the improper Riemann integral of f between two filters F and G. The function is_RInt_gen is in fact a relation. Even with a proof that this relation locally describes a function, it is not possible to directly mention the value of this function, i.e. the integral. For example, we can write in Coquelicot Bounding the remainder with a low accuracy is sufficient to prove that the integral on the whole domain is included in [226.849; 226.850] and thus that the upper bound 226.844 used in [START_REF] Helfgott | Major arcs for Goldbach's problem[END_REF] is incorrect, although this does not invalidate the proof. Here we present benchmarks from both the proper and the improper case for estimating integrals.

-1 n (the inverse is considered as belonging to the class of usual functions here, and so it is treated below)

• Let f be any of the other (unary) operators in E. It can be checked that the first derivative of f is also an expression d f ∈ E. Let us define f n recursively:

where :: denotes concatenation.

We need to justify the correctness of this construction. The first two bullet points are straightforward. The third bullet point only replaces a division by a product with the inverse. For the fourth bullet point, suppose that for some g and x,

which is exactly what we wanted.

Remark 4. This technique generalizes to the case when we have a differential equation on f of the shape f (k+1

, where g is defined as an expression of E, and some initial conditions at x: y

We have implemented an interval version of this method in Coq, on top of CoqInterval. It has not been proved correct yet, however it was used to perform benchmarks on Simpson's method for quadrature.

Appendix A

Complexity of Some Matrix Multiplication Algorithms

The following tables were deferred to the Appendix from Chapter 4. Table A.2 concerns the following tensor from Pan, which we present as a table: Abstract : This thesis presents three contributions to the topic of computer-assisted proofs: both in the case of proofs relying on computations and of formal proofs produced and verified using a piece of software called a proof assistant. We first illustrate the theme of experimentation at the service of proofs by programming a symbolic manipulation tool and using it to find and demonstrate a result about matrix multiplications.

In a second part, we use a proof assistant to write a proof of Apéry's theorem, which is based on complex computations. These computations are performed by an efficient external program and are then verified by this proof software.

In the third contribution, we propose a program which computes integrals and simultaneously builds proof that the result is correct. We use this program to find errors in published results.