Christophe Riedinger

Michel Jordan
email: michel.jordan@u-cergy.fr

Hedi Tabia
email: hedi.tabia@ensea.fr

3D MODELS OVER THE CENTURIES: FROM OLD FLOOR PLANS TO 3D REPRESENTATION

Keywords: Cultural heritage, 3D Restitution, Image Processing, Image Segmentation, 3D Application

This paper presents a set of algorithms dedicated to the 3D modeling of historical buildings from a collection of old architecture plans, including floor plans, elevations and cutoffs. Image processing algorithms help to detect and localize main structures of the building from the floor plans (thick and thin walls, openings). The extrusion of the walls allow us to build a first 3D model. We compute height informations and add textures to the model by analyzing the elevation images from the same collection of documents. We applied this pipeline to XVIII th century plans of the Château de Versailles, and show results for two different parts of the Château.

INTRODUCTION

3D models are very common in a wide area of domains, including design and architecture, industries, entertainments, etc. Also in the cultural heritage domain, they are more and more useful for archeologists and historians, as well as curators. Most of 3D models for culture heritage are built from laser scanning and/or photogrammetry [START_REF] Cosmas | 3D MURALE: A Multimedia System for Archaeology[END_REF][START_REF] Yastikli | 3D model generation and visualization of cultural heritage[END_REF][START_REF] Tingdahl | Arc3d: A public web service that turns photos into 3d models[END_REF][START_REF] Remondino | Heritage Recording and 3D Modeling with Photogrammetry and 3D Scanning[END_REF] (in [START_REF] Van | Special issue on 3d acquisition technology for cultural heritage[END_REF], interested readers will find a review of various techniques which were used in the beginning of the XXI st century). Recent 3D laser scanning techniques combined with photography augmented reality algorithms, allow to build very impressive models, involving accurate geometry and realistic texture mappings, as in [START_REF]Three dimensional monuments, programme national de numérisation du patrimoine[END_REF] which main objectives are the preservation and the promotion of the cultural heritage. Laser scanning is also combined with CAD softwares and image processing techniques in order to build 3D models of partially disappeared monuments [START_REF]Cluny numérique[END_REF][START_REF] Landrieu | Reconstitution virtuelle de l'église abbatiale Cluny III: des fouilles archéologiques aux algorithmes de l'imagerie[END_REF]. These techniques involve a great amount of manual operations not only in the 3D digitization steps but also in the virtual/augmented reality steps.

But what could be done when the building or monument has been partially or completely destroyed, or when a lot of transformations and changes have modified its appearance or indoors? Usually, a lot of old floor plans and other graphic documents help historians to know these transformations, and to imagine the disappeared appearance of the most important monuments.

The aim of our work is to propose a pipeline for 3D model generation from old floor plans of historical castles, oldmansions or other interesting monuments. Among the difficulties of the task we are coping with, are some specific features of these old documents (cf. fig 1) due to drawing as well as digitization: straight lines are not so straight, writings look like hand-writings, line thickness is not equal along the same line, symbols are not homogeneous, paper watermarks and foldings appear in the images, etc. The reminder of the paper is organized as follows: in section 2, we give a brief review of automatic analysis of architectural plans; section 3 gives the method overview, the 2D floor plan processing is detailed in section 4, and section 5 presents the 3D model extrusion algorithm. In section 6, we present some 3D model captures of the Château de Versailles, and conclusions and future developments end the paper (section 7).

RELATED WORK

The use of 3D representation is very helpful for architects. It gives them a more intuitive perspective on their work, and moreover allows to conduct large sets of studies related to lighting, heat transfer, radio waves propagation, etc. Automatic generation of 3D models from 2D architectural drawings is a quite active research field since the end of XX th century [START_REF] Xuetao | Generating 3D Building Models from Architectural Drawings: A Survey[END_REF]. Architectural floor plans are generally composed of geometric shapes (straight lines, curves) representing the building structure (external and inner walls), while various symbols represent openings, heating and furniture elements, etc. A lot of textual annotations gives information about the building use, rooms, dimensions. . . The most common automatic approaches consist first in seperating graphics and texts in the images. The detection of straight lines (segments) allow to localize thick and thin walls, either by using mathematical morphology or Hough transforms [START_REF] Dosch | A complete system for the analysis of architectural drawings[END_REF][START_REF] Macé | A System to Detect Rooms in Architectural Floor Plan Images[END_REF]. 3D models can be built therefore by analyzing the connectivities of the segments and searching for their internal relationships [START_REF] Or | Highly Automatic Approach to Architectural Floorplan Image Understanding & Model Generation[END_REF]. Key symbols are recognized [START_REF] Guo | An improved example-driven symbol recognition approach in engineering drawings[END_REF] and interpreted [START_REF] Zhu | A new reconstruction method for 3D buildings from 2D vector floor plan[END_REF]. Looking for approximatively convex regions [START_REF] Macé | A System to Detect Rooms in Architectural Floor Plan Images[END_REF] or SURF key-point extraction [START_REF] Ahmed | Automatic analysis and sketch-based retrieval of architectural floor plans[END_REF] help users to find room boundaries, and applying OCR to the text layer gives the room functions [START_REF] Ahmed | Automatic analysis and sketch-based retrieval of architectural floor plans[END_REF]. Finally, walls are extruded in order to provide the 3D model of the building. In [START_REF] Tong | A new recognition model for electronic architectural drawings[END_REF], Lu et al. present a system for generating 3D models from "architect drawings", i.e. sketches ordinary used by the architects to describe their primary works. These drawings are divided in 3 categories: structural (building structures), functional (doors, windows, elevators, etc.), and decorative (dividing walls, lights, etc.). The presented system uses 3 main steps to process the drawings: shape analysis and recognition allow to detect structural elements; the graphic primitives corresponding to the recognized shapes are removed; and finally a symbol recognition algorithm detect architectural entities.

However, as explained in section 1, we are facing some drawbacks that do not allow us to directly use these algorithms to build our culture heritage 3D models. In particular, we will have to cope with not so straight "straight lines", unnormalized symbols and texts. Nevertheless, we will adapt some of the global architecture of these systems in our context, as we will explain in next section.

METHOD OVERVIEW

We propose an automated method to perform 3D reconstruction of monuments using floor plans designed by architects. The different steps of the process are summarized on figure 2. The first part of the pipeline uses image processing algorithms in order to extract relevant informations from the floor plans (section 4): after some preprocessings (section 4.1) in order to remove noise, and reduce the amount of information (the floor plan digitized images are usually more than 8000 × 8000), we try do detect, localize and classify the building walls in the floor plan (section 4.2). We then extrude the 3D model on the basis of the walls (section 5).

The second part of the pipeline deals with the elevation and cut-off images (section 5.2); these images are used to fix the height of the 3D model, localize openings and add textures to the walls.

2D FLOOR PLAN PROCESSING

In this section, we present the 2D processing pipeline used in this work in order to extract relevant informations from 2D plans. We will particularly detect load-bearing walls as well as dividing walls, that we will next extrude in order to generate the 3D model.

Preprocessings

Image resizing

Most old architectural plans have significant size, up to 2 meters each side. Moreover, the digitization of these plans is performed with a high resolution. Depending on the resolution of the input plan and in order to increase the processing speed, the image dimensions can be reduced by a given factor r. Image reduction is performed in the following way: we keep the darkest pixel in each r ×r vignette. Indeed, the plans that we want to extract are drawn with black pencil on white paper; thus, the dark gray levels carry the useful information for our purpose.

Denoising and binarization

Before the binarization step, we perform a denoising processing, using the Non Local Means algorithm [START_REF] Van De Ville | Sure-based non-local means[END_REF]. After that the binarization step works as follows: we first compute the median value med of the gray level image I ; as most of the plan is white paper, med is close to the background white color of the paper. For each pixel (u, v), we compute, in its 3 × 3 neighborhood, the mean value of the absolute differences between the pixel's gray level and med, and create an image I dif f . Each pixel of the binarized image I b is computed as follows (let s a fixed threshold parameter):

if I(u, v) > med I b (u, v) = 0 if I dif f (u, v) > s × med dif f I b (u, v) = 255 otherwise I b (u, v) = 0
We chose the s value experimentally, and checked that the binarization process is robust against its variations. The binarized image shows well contrasting edges that will be easier to detect.

Detection and localization of building walls

We will detect the main walls and the dividing walls in the binarized image; in this image, the walls seam as dark strokes with different thicknesses, thus we use mathematical morphology operators to detect them. This allow to use the same algorithm to detect thick walls and thin ones, only by tuning a threshold parameter. Because of the "not so straight" character of straight lines in these documents, using Hough transforms would not give relevant results.

Wall detection

We iterate a unit ball morphological erosion process, thus obtaining an image where each pixel belonging to a connex

Wall segmentation

The next step consists in computing the contours of the detected walls, which will be the edges of the final 3D model. We first compute the thresholded local gradient image I grad of the wall image. We obtain the morphological skeleton of I grad by means of the 8 L structuring elements [START_REF] Soille | Morphological Image Analysis: Principles and Applications[END_REF]; finally, we localize in this image the extremities of the segments (red points in figure 4, left).

To extract primary segments we perform the following procedure: we assign a label to each point in the 5 × 5 neighborhood of a junction if the point is at a two pixel distance of the junction (the labels and the junctions are not first neighbors). The labels are dilated iteratively with the unit ball as a structuring element. We do not modify the image points not belonging to the skeleton. After several iterations, when the image does not change any more, each label covers a part of the skeleton that is a segment between two junctions. For each segment S, beginning with the longest segment and finishing with the smaller one, we merge the adjacent segments to S using the following criteria: CRITERIA -Let A be a segment and B the segment that is to be combined with A; the segments are merged if and only if the distance between the barycenter of the points of B and the line passing through A is less than a given threshold.

We iterate the criteria until there is no more candidate segments. The figure 4 (right) shows segments after fusion on a small part of the original image.

Detection and localization of staircases

Staircases are an important architectural and functional element of buildings, and we defined a semi-automatic method to detect them in the floor plans. The semi-automatic method allows to deal with the unnormalized symbols in the old floor plans; we found a semi-automatic method better than a fully automated one, which would result in more erroneous detections. The user first designs a rectangular area around stairs, and the algorithm detect its shape and accurately localize it. In fact, it is easier for the user to design areas where to localize stairs than to post-process a fully automatic detection and remove false detections. Moreover the localization is faster in this way, because it deals only with small vignettes in the image.

In our plans, staircases are represented with plain strokes (for the stairs at the given floor) and hashed strokes (for the ones belonging to the up/bottom floor), see figure 5.

The localization process works as follows: first we detect the different connex components in the set stairs-walls of the floor plan (these are the black pixels figure 6, left); by connex component, we mean a set of contiguous black pixels which are fully disconnected from the other components because white pixels surround them (please note that a dash is a connex component by itself). We then try to merge connex components that belong to the same implicit line segment. In this purpose, we use the well-known Tukey estimator [START_REF] Beaton | The fitting of power series, meaning polynomials, illustrated on bandspectroscopic data[END_REF] to compute the best line approximation for a given connex component; each dash is represented by a connex component and its best line approximation (its angle with the horizontal axis). We store the x axis angles of all dashes in a circular histogram (each bin of the histogram represents a 5 • deviation), then we merge all the dashes in the same bin of the histogram. This merging process is divided in two cases: (i) for the angles of 0 • (resp. 90 • ), we horizontally (resp. vertically) prolongate the dashes; (ii) for the oblique lines, we take all the pixels of the dashes falling in the same bin and compute the angle of the implicit line corresponding with these dashes with the horizontal axis, thus obtaining a better approximation of this line. As a result, the white regions in figure 6 (right) are fully surrounded by the implicit lines joining the dashes; these white regions represent the stairs of the stairwell (figure 6). If two adjacent connex components have two common points, we state that they represent two stairs with a high probability; thus we keep all the pixels belonging to these components, and their bounding box defines the bounding box of the staircase (figure 6, right).

3D RESTITUTION

In section 4, we analyzed the floor plans and obtained a 2D vectorial representation of the main structures of the building. In the following paragraphs, we will build a 3D representation from the footprints detecting in the floor plans and using informations from the elevation and cut-off views. We do not intend to build a realistic 3D model, but only to give it a "pleasant" appearance in order that the user (particularly historians) could better understand the building. Nevertheless, even in very well-documented monuments such as the Château de Versailles, we will rarely find parts of the monument with a complete description for all the walls of the rooms.

3D extrusion

A first 3D model is computed by simply extruding the walls on a given height (figure 7); the height of the walls is arbitrary fixed, and will be tuned thanks to the analysis of the elevation images.

Texturing the 3D model

At this point of the process, each wall of the 3D model consists of two adjacent triangles constituting a rectangle. Additionally to the floor plans there are also several elevation plans of the castle representing vertical views of the rooms and the outside wall. The goal of the work presented in this section is to use the elevation plans (figure 8) to refine the vertical 3D representation of the walls of the monument. The user interactively selects a wall or a set of walls on the primary 3D model and an approximatively corresponding rectangular region on the elevation image (figure 9). A vertical mesh of the selected walls is then automatically constructed and the older walls of the 3D primary model are replaced by the newly formed polygons. The elevation plan is also used to texture the new polygons. The vertices of the polygon are obtained from the points of the selected region in the elevation plan and this set of vertices is triangulated.

In order to compute the polygonalization of the vertical faces of the 3D model, we compute a 2D triangulation of the elevation image. We detect the contours of the interesting part of the elevation map thanks to a flooding algorithm, starting We used the same algorithm to analyze elevation images representing the inner sides of the rooms, in order to detect doors, windows, fireplaces, or other remarkable objects.

RESULTS

We first applied the algorithms presented in this paper to images of the Château de Versailles, thanks to a collaboration with the Research center of the Château de Versailles (CRCV), the french National Archives and the french National Library (BnF). We present in this section results computed on two different parts of the Château: the chapel and the "Appartement du Roy".

Château de Versailles

The algorithm were developed in Java, and tested on a bi-Xeon 2.4 GHz processor. We present in table 1 the computation times for the 3D modeling of the "Appartement du Roy" (original floor plan image resized to 2600 × 3700).

We show on figure 11 different views of the "Appartement du Roy", either general or detailed view of the "Chambre du Roy", the 3D model is textured with the available data for two walls of the "Chambre" (bedroom). In this room, there is a balustrade in front of the king's bed, that we detected in the floor plan as a dividing wall, we fixed its height by computing it in comparison of the height of the room on the elevation image (figure 11, lower-right). The 3D model of the Chapel is in figure 12.

Results on other data

We also checked our algorithms on data from other sources. For instance, we present in figure 13 the result of the floor plan analysis of a house in Aix-en-Provence (document from the Bibliothèque nationale de France). The size of this floor plan is significantly lower than the size of the floor plan of Versailles, thus we did not perform a size reduction. These results show that our algorithms are robust enough to be used on various data sources.

CONCLUSION AND PERSPECTIVES

We presented here a pipeline of an automated process towards 3D models of historical buildings such as royal castles from old floor plans. Image processing algorithms based on morphological mathematics allow to extract 2D informations from the plans, such as exterior and inner walls. Dedicated tools filter text areas.

Documents not only involves plans, but also cross-sections and elevations. A database will register links between the various documents describing each part of the castle, and links between 2D features computed from the digitized plans, thus leading to 3D models of different rooms of the monument. Nevertheless, historians know that these old plans are not technical drawings as architects use nowadays, and that there was often a "gap" between these drawings and the actual building, especially concerning the proportions of the rooms; thus we do not intend to produce "hyper-realistic" 3D views of the buildings, but only to give quite a good idea of their appearance and their internal organization at the time of the plans.

These 3D models will help historians to better understand the life in these castles in the past centuries; in fact, longliving historical buildings have changed over the years, due to the various usage of them (habitation, royal power place, national palace, . . . ), parts of the monuments have been destroyed and rebuilt. Historians have documents about the life in these buildings which can not be explained by looking at their current appearance and organization. We hope that 3D models of their disappeared state(s) will provide a great benefit for those who are working on these old documents. Moreover, 3D models could also be interesting for general purpose, such as virtual visits: people will be able to visit a castle as it was in the past centuries, according to the documents that architects of these times has left for us!

Fig. 1 .

 1 Fig. 1. Floor plans of an historical castle (XVIII th century).

Fig. 2 .

 2 Fig. 2. Method overview.

Fig. 3 .

 3 Fig. 3. Walls detected in the floor plan: (left) main walls; (right) dividing walls.

Fig. 4 .

 4 Fig. 4. Segment junctions (left), and segments after fusion (zoom, right figure).

Fig. 6 .

 6 Fig. 6. Detection of staircase: results of the localization algorithm.

Fig. 7 .

 7 Fig. 7. 3D model of the walls detected in the original plan.

Fig. 5 .

 5 Fig. 5. Detection of staircase: the user first gives an approximate localization of the staircase.

Fig. 8 .

 8 Fig. 8. Original elevation image.

Fig. 9 .

 9 Fig. 9. The delimited boundary of the polygon.

Fig. 10 .

 10 Fig. 10. 3D model with one elevation wall.

Fig. 11 .

 11 Fig. 11. 3D model of the "Appartement du Roy", with general views and details on the "Chambre du Roy".

Fig. 12 .

 12 Fig. 12. 3D model of the Chapel of the Château de Versailles.

Fig. 13 .

 13 Fig. 13. Floor plan of the Hôtel de Caumont, and 3D model.

Table 1 .

 1 Computation times.

	Image size	Preprocessing Segmentation Walls Staircase Extrusion Textures
	2600 × 3700	10 s	25 s	61 s	1,5 s	1 s	20 s

Acknowledgments

VERSPERA is a research project funded by the PATRIMA french Labex and the "Fondation pour les Sciences du Patrimoine" (Sciences for Material Cultural Heritage). We particularly thank the Centre de recherche du Château de Versailles, the Archives nationales and the Bibliothèque nationale de France, for providing us images, data and historical informations.


Sylvie Boldo
email: sylvie.boldo@inria.fr

Jean-Christophe Filliâtre
email: christophe.filliatre@lri.fr

Guillaume Melquiond
email: guillaume.melquiond@inria.fr

Combining Coq and Gappa for Certifying Floating-Point Programs

Formal verification of numerical programs is notoriously difficult. On the one hand, there exist automatic tools specialized in floatingpoint arithmetic, such as Gappa, but they target very restrictive logics. On the other hand, there are interactive theorem provers based on the LCF approach, such as Coq, that handle a general-purpose logic but that lack proof automation for floating-point properties. To alleviate these issues, we have implemented a mechanism for calling Gappa from a Coq interactive proof. This paper presents this combination and shows on several examples how this approach offers a significant speedup in the process of verifying floating-point programs.

Introduction

Numerical programs typically use floating-point arithmetic [START_REF]Microprocessor Standards Subcommittee: IEEE Standard for Floating-Point Arithmetic[END_REF]. Due to their limited precision and range, floating-point numbers are only an approximation of real numbers. Each operation may introduce an inaccuracy and their total contribution is called the rounding error. Moreover, some real operations may not be available as sequences of floating-point operations, e.g., infinite sums or integrals. This introduces another inaccuracy called the method error. Both errors make it somehow complicated to know what floating-point programs actually compute with respect to the initial algorithms on real numbers.

One way to proceed is to give a program a precise specification of its accuracy. Generally speaking, a specification explains what can be expected from the result given facts about the inputs. Typically, it bounds the sum of both rounding and method errors. For example, the specification for a function float cos defined on the double type of floating-point numbers may be the following:

3 ∀x : double, |x| ≤ 2π ⇒ cos(x) -float cos(x) cos(x) ≤ 2 -53 .
Such an inequality is typically proved using pen and paper. It can be deduced from the specification of each floating-point operation and the mathematical properties of the cosine function. Such proofs are notoriously difficult and errorprone.

Ideally, the code of a software would be analyzed by a certification tool, which would answer whether it is correct or not. In order to increase the confidence in the answer (and therefore in the software), the tool may rely on formal methods. Obviously, such a tool is not conceivable, but we aim at making this process as automatized as possible.

In this paper, we will describe how we have combined existing tools to help in the process of formally certifying numerical codes. From a user point of view, the first step is to annotate the source code with the specifications of the software. This annotated code is sent to a first tool that produces proof obligations corresponding to the correctness of the software. The examples of this paper are C programs and we are using the Caduceus tool (Section 2.1). By computing weakest preconditions based on the code and the annotations, it generates theorem statements to be verified by automated theorem provers or proof assistants. Some theorems, hopefully most of them, will be automatically discharged. For instance, numerical properties may be discharged by the Gappa tool (Section 2.3), which is very efficient at proving bounds, especially on rounding errors. But Gappa only tackles the floating-point fragment of a program, so properties that involve more than just floating-point arithmetic may not be handled.

The remaining theorems will have to be manually handled by the user in a proof assistant with a suitable floating-point formalization. The paper focuses on the use of Coq for this task (Section 2.2). When using a proof assistant, the user issues tactics to split the goal into simpler subgoals. Examples of such tactics are logical cut, case analysis, or induction. This will become tedious if the user has to repeat this process until all the subgoals are discharged, which may require a high number of explicit proof steps. Especially frustrating is the fact that, once simplified, the subgoals may fit into specific logic fragments, which some tools, such as Gappa, could handle automatically outside Coq.

In order to benefit from Gappa inside Coq, we have implemented a mechanism for calling the tool from an interactive proof. From a technical point of view, Gappa is called as an external prover. This does not weaken the confidence in Coq formal proofs since Gappa produces a proof trace that is checked by Coq (Section 4).

This combination of Coq and Gappa does not radically change the way to tackle rounding and method errors. It simply eases the use of traditional approaches in a formal setting. C programs illustrating this point are given in Section 3. The combination of all these tools (Caduceus, Coq, Gappa) makes it possible to formally verify a source code while benefiting from automation.

There have been previous work on formally proving numerical components (especially hardware ones) while relying on automated tools. Among them, the certification of the IEEE-compliance [START_REF]Microprocessor Standards Subcommittee: IEEE Standard for Floating-Point Arithmetic[END_REF] of a gate-level design for the Pentium Pro processor used Forte, a combination of two model checkers and a lightweight theorem prover [START_REF] O'leary | Formally verifying IEEE compliance of floating-point hardware[END_REF]. Another work made the ACL2 theorem prover interact with a VHDL verification tool in order to prove the correctness of a hardware multiplier [START_REF] Reeber | Combining ACL2 and an automated verification tool to verify a multiplier[END_REF]. While not at the source-code level, two other proof assistants have been thoroughly used for verification of floating-point properties: both PVS and HOL Light provide some automation for performing error analysis [START_REF] Carreño | Specification of the IEEE-854 floating-point standard in HOL and PVS[END_REF][START_REF] Harrison | A machine-checked theory of floating-point arithmetic[END_REF].

2 Caduceus, Coq, and Gappa

Verification of C Programs

The Why platform4 is a set of tools for deductive verification of Java and C programs [START_REF] Filliâtre | The Why/Krakatoa/Caduceus Platform for Deductive Program Verification[END_REF]. In this paper, we only focus on verification of C programs but the results would apply to Java programs as well. The verification of a given C program proceeds as follows. First, the user specifies requirements as annotations in the source code, in a special style of comments. Then, the annotated program is fed to the tool Caduceus [START_REF] Filliâtre | Multi-prover verification of C programs[END_REF], which is part of the Why platform, and verification conditions (VCs for short) are produced. These are logical formulas whose validity implies the soundness of the code with respect to the given specification. Finally, the VCs are discharged using one or several theorem provers, which range from interactive proof assistants such as Coq to purely automatic theorem provers such as Alt-Ergo [START_REF] Conchon | CC(X): Semantical Combination of Congruence Closure with Solvable Theories[END_REF]. The workflow is illustrated on Figure 1. Annotations are inserted in C source code using comments with a leading @. They are written in first-order logic and re-use the syntax of side-effect free C expressions. For instance, here is a code excerpt where an array t is searched for a zero value.

// @ invariant 0 <= i ❢♦r ( i = 0; i < n ; i ++) { ✐❢ ( t [ i ] == 0) ❜r❡❛❦ ; } // @ assert i < n = > t [ i ] == 0
The for loop is given a loop invariant, as in traditional Hoare logic [START_REF] Hoare | An axiomatic basis for computer programming[END_REF]. (In that case, the invariant could be found automatically.) A loop invariant typically generates two VCs: one to show that it holds right before the loop is entered; and one to show that it is preserved by the loop body. In this example, an assertion is also manually inserted right after the loop, which results in a VC for this program point. Additional VCs are produced to establish the safe execution of the code, i.e., that the program does not perform any division by zero or any array access out of bounds. In this example, a VC requires to show that t[i] is a legal array access, which may or may not be provable depending on hypotheses regarding t and n. Verification with Caduceus is modular: each function is given a contract and proved correct with respect to the contracts of the functions it calls. 5 For instance, a partial contract for a function sorting an array of integer could be /* @ requires @ 0 <= n && \ valid_range (t , 0 , n -1) @ assigns @ t [0.. n -1] @ ensures @ \ forall int i , j ;

0 <= i <= j < n = > t [ i ] <= t [ j ] */ ✈♦✐❞ sort ( ✐♥t *t , ✐♥t n );
The contract contains three clauses. Keyword requires introduces a precondition, that is a property assumed by the function and proved at the caller site. In this example, it states that n is nonnegative and that all indices from 0 to n-1 in t can be safely accessed. Conversely, keyword ensures introduces a postcondition, that is a property provided by the function, right before it returns. Here it states that the array is sorted in increasing order. 6 Finally, keyword assigns introduces the memory locations possibly modified by the function, which means that any other memory location is left unchanged by a call to this function. Here, it states that only the array elements t[i] for 0 ≤ i < n are possibly assigned. Caduceus handles a large fragment of ANSI C, with the notable exception of pointer casts and unions. It handles floating-point arithmetic, using a model where each floating-point number is seen as a triple of real numbers [START_REF] Boldo | Formal Verification of Floating-Point Programs[END_REF]. The first component is the floating-point number itself, as it is computed. The second component is the real number that would have been computed if roundings were not performed. The third component is a ghost variable attached to the floating-point number and which represents the ideal value that the programmer intended to compute. Annotations are written using real numbers only, and the three components of a floating-point variable x can be referred to within annotations: x itself stands for the first component; \exact(x) for the second one; and \model(x) for the third one. Thus the user can refer to the rounding error as the difference between the first two, and to method error as the difference between the last two. Examples are given in Section 3.

Since the general-purpose automatic provers do not support this model of floating-point arithmetic, we have formalized it in the Coq proof assistant.

The Coq Proof Assistant

The Coq proof checker [START_REF] Bertot | Interactive Theorem Proving and Program Development. Coq'Art: The Calculus of Inductive Constructions[END_REF][START_REF]The Coq Development Team: The Coq Proof Assistant Reference Manual -Version V8.2[END_REF] is a proof assistant based on higher-order logic. One may express properties such as "there exists a function which has such and such properties" or "every relation that verifies such hypothesis has a certain property" and check proofs about these. Proofs are built using tactics (such as applying a theorem, rewriting, computing, etc.). A Coq file contains the statement of lemmas and their proofs as a sequence of tactics in the Coq language.

The Coq standard library contains an axiomatization of real numbers [START_REF] Mayero | The Three Gap Theorem (Steinhauss conjecture)[END_REF]. Few automation is provided to reason about real numbers. As a consequence, the proof of a typical lemma such as 0 < 1 -2 -52 is already a few lines long: The proof is done backward, by transforming the conclusion until it trivially derives from the hypotheses. This proof starts by applying the theorem Rlt Rminus (line 3) since 0 < 1 -2 -52 is a consequence of 2 -52 < 1. The definition of powerRZ is then unfolded (line 4) so that 2 -52 is converted to (2 52 ) -1 . We replace 1 by 1 -1 (theorem Rinv 1, line 5). At this point, the goal has become (2 52 ) -1 < 1 -1 . After applying theorem Rinv 1 lt contravar (line 6), the remaining goals are 1 < 2 52 and 1 ≤ 1, which are solved automatically by the tactic auto (line 6).

1 ▲❡♠♠❛ OneMinusUlpPos : (0 < 1 -powerRZ 2 ( -52))% R . 2 Pr♦♦❢ .
A high-level formalization of floating-point arithmetic [START_REF] Daumas | A generic library of floating-point numbers and its application to exact computing[END_REF][START_REF] Boldo | Preuves formelles en arithmétiques à virgule flottante[END_REF] is also available in Coq. A floating-point number is a pair of integers (m, e) which represents the real number m × 2 e . The value of the mantissa m and the exponent e are bounded according to the floating-point format. For example, in IEEE-754 double-precision format [START_REF]Microprocessor Standards Subcommittee: IEEE Standard for Floating-Point Arithmetic[END_REF], the pair verifies |m| < 2 53 and -1074 ≤ e. This library7 contains a large number of floating-point definitions and theorems and has been used to prove many old and new properties [START_REF] Boldo | Pitfalls of a Full Floating-Point Proof: Example on the Formal Proof of the Veltkamp/Dekker Algorithms[END_REF].

Most floating-point proofs rely on computations on real numbers, such as deciding 0 < 1 -2 -52 or bounding method error. Such goals can be addressed using the interval tactic [START_REF] Melquiond | Proving bounds on real-valued functions with computations[END_REF]. This reflexive tactic, based on interval arithmetic, decides inequalities by bounding real expressions thanks to guaranteed floatingpoint arithmetic. Once done with method error, the user is left with VCs related to rounding errors, which Gappa is typically designed for.

The Gappa Tool

Gappa8 is a tool dedicated to proving arithmetic properties on numerical programs [START_REF] De Dinechin | Assisted verification of elementary functions using Gappa[END_REF][START_REF] Melquiond | De l'arithmétique d'intervalles à la certification de programmes[END_REF]. Given a logical proposition expressing bounds on mathematical real-valued expressions, Gappa checks that it holds. The following is such a proposition and below is its transcription in Gappa's input language.

∀x, y ∈ IR, |x| ≤ 2 ∧ y ∈ [1, 9] ⇒ x × x + √ y ∈ [1, 7] ④ | x | <= 2 ✴❭ y ✐♥ [1 ,9] ✲❃ x * x + sqrt ( y ) ✐♥ [1 ,7] ⑥
In order to verify the proposition, Gappa first analyzes which expressions may be of interest. Then it tries to enclose them in intervals by performing a saturation over its library of theorems on interval arithmetic, forward error analysis, and algebraical identities. Gappa stops when it reaches enclosures small enough to be compatible with the right-hand side of the proposition or when the saturation does no longer improve the enclosures.

Once Gappa has verified the proposition, it generates a formal proof. To increase confidence, this proof script can then be mechanically checked by an independent proof system, such as Coq or HOL Light. 9If Gappa fails to prove the proposition, the user can suggest to the tool that it should perform a bisection-splitting input intervals until the proposition holds on each sub-intervals-or augment the library of theorems with new mathematical identities. Gappa will then assume these equalities hold; they will appear as hypotheses of the generated formal proof.

Expressing Floating-Point Programs. In addition to universally-quantified variables on IR, basic arithmetic operators (+, -, ×, ÷, √ •), and numerical constants, Gappa expressions can also contain rounding operators. The integerpart functions, ⌊•⌋ and ⌈•⌉, are instances of such operators. Since the IEEE-754 standard [START_REF]Microprocessor Standards Subcommittee: IEEE Standard for Floating-Point Arithmetic[END_REF] mandates that "a floating-point operator shall behave as if it was first computing the infinitely-precise value and then rounding it so that it fits in the destination floating-point format", having appropriate rounding operators is sufficient to express the computations of a floating-point program.

The following script is similar to the previous one, but all the expressions are now as if they had been computed in single precision with rounding to nearest (tie-breaking to even mantissa).

@rnd = float < ieee_32 , ne >; z = rnd ( rnd ( x * x ) + rnd ( sqrt ( y )));

④ | x | <= 2 ✴❭ y ✐♥ [1 ,9] ✲❃ z ✐♥ [1 ,7] ⑥
Note that Gappa only manipulates expressions on real numbers. As a consequence, infinities and NaNs (Not-a-Numbers) are no part of this formalism: rounding operators return a real value and there is no upper bound on the magnitude of the input numbers. This means that NaNs and infinities will not be generated nor propagated as they would in IEEE-754 arithmetic.

We can, however, use Gappa to prove that a given code does not produce any of these exceptional values. Indeed, if one proves that a Gappa-rounded value is smaller than the biggest floating-point number in the working format, then the actual IEEE-754 computation is guaranteed not to overflow, by definition of overflow. Therefore, in order to check that computations in the previous example are overflow-safe, one can run Gappa on the following script:10 @rnd = float < ieee_32 , ne >; z = rnd ( rnd ( x * x ) + rnd ( sqrt ( y )));

④ | x | <= 2 ✴❭ y ✐♥ [1 ,9] ✲❃ z ✐♥ [1 ,7] ✴❭ | rnd ( x * x )| <= 0 x1 . FFFFFEp127 ✴❭ | rnd ( sqrt ( y ))| <= 0 x1 . FFFFFEp127 ⑥
The absence of infinities and NaNs is not a deficiency, as reasoning about them is usually done by case analysis. This can be easily performed using Coq traditional tactics. For the cases without infinities and NaNs, which are the complicated ones, the Gappa tactic applies.

Verifying Accuracy. While the previous examples show that Gappa can bound ranges of floating-point variables, this is only a small part of its purpose. This tool was designed to prove bounds on computation errors, which also happen to be real-valued expressions. Let us assume that the developer actually needed the infinitely-precise result M z = x 2 + √ y. Is the computed result z sufficiently close to this ideal value M z ? This can be answered by bounding the absolute error z -M z :11 @rnd = float < ieee_32 , ne >; Mz = x * x + sqrt ( y ); z = rnd ( rnd ( x * x ) + rnd ( sqrt ( y )));

④ | x | <= 2 ✴❭ y ✐♥ [1 ,9] ✲❃ | z -Mz | <= 1b -21 ⑥
For the sake of simplicity, M z has the same operations as z, but without rounding. This is not a requirement, as Gappa is also able to bound errors when M z is a completely different expression. Note also that Gappa is not limited to absolute errors; it can handle relative errors in a similar way, which is especially important when proving floating-point properties.

Proving Floating-Point Programs

Before describing the inner working of the Coq-Gappa combination, we illustrate its use on the verification of three typical floating-point programs.

Naive Cosine Computation

The first example is an implementation of the cosine function for single-precision floating-point arithmetic. To present a complete Coq proof, we have simplified the function by removing its argument-reduction step. Thus, input x is required to have already been reduced to a value close to zero; only the polynomial evaluation has to be performed. The specification of the function states that, for |x| smaller than 2 -5 , the computed value \result is equal to cos(x) up to 2 -23 . Given this annotated C code, Caduceus generates a VC stating the accuracy of the result, which can be formally proved with the Coq script below. 121 Pr♦♦❢ . The first part of the proof script (line 2) turns the goal into a user-friendly form: the why2gappa tactic cleans the goal by expanding and rewriting some Caduceus-specific notations. At this point, assuming that •(•) is the rounding operation from a real number to the nearest single-precision floating-point number, the user has to prove the following goal:

∀x : float, |x| ≤ 1 32 ⇒ | •(•(1) -•(•(x × x) × •(5/10))) -cos(x)| ≤ 2 -23 .
As would be done with a pen-and-paper verification, the formal proof of this goal starts by computing and proving a bound on the method error. Since the polynomial is chosen so that the computed result is close to the cosine, the method error is known beforehand. A typical way to obtain a polynomial approximation and its error is to use a computer algebra system.

Here, the method error is smaller than 7 × 2 -27 . So we are asserting this property in Coq (lines 3 and 4) and we prove it (line 5). The assertion is proved by the interval tactic [START_REF] Melquiond | Proving bounds on real-valued functions with computations[END_REF]. Its option i bisect diff tells the tactic to recursively perform a bisection on the interval enclosure [-2 -5 , 2 -5 ] of x, until a first-order interval evaluation of the method error (1 -(x × x) × (5/10)) -cos(x) gives a compatible bound on all the sub-intervals.

Once the assertion is proved and hence available as an hypothesis, the user has to prove the following property:

∀x : float, |x| ≤ 1 32 ⇒ |(1 -(x × x) × (5/10)) -cos(x)| ≤ 7 • 2 -27 ⇒ | •(•(1) -•(•(x × x) × •(5/10))) -cos(x)| ≤ 2 -23 .
This is achieved by the gappa tactic (line 6). It calls Gappa and then uses the Coq script that the tool generates in order to finish the proof. Note that the Gappa tool takes advantage of the inequality proved by the interval tactic as it knows nothing about the cosine.

Discretization of a Partial Differential Equation

The second example is a numerical code about acoustic waves by F. Clément [START_REF] Bécache | Étude de schémas numériques pour la résolution de l'équation des ondes[END_REF]. Given a rope attached at its two ends, a force initiates a wave, which then undulates according to the following mathematical equation:

∂ 2 u(x, t) ∂t 2 -c 2 ∂ 2 u(x, t) ∂x 2 = 0.
The value u(x, t) gives the position of the rope at the abscissa x and the time t. It is discretized both in space and time with steps (∆x, ∆t). The result is a matrix p of size ni × nk where p[i][k] = p k i is the position of the rope at the abscissa i × ∆x and the time k × ∆t. The matrix p is computed by the following piece of code [START_REF] Boldo | Floats & Ropes: a case study for formal numerical program verification[END_REF], where a is an approximation of an exact constant A derived from c, ni, and nk:

/* @ invariant 1 <= k <= nk && analytic_error (p , ni , ni ,k , a ) */ ❢♦r ( k =1; k < nk ; k ++) { p [0][ k +1] = 0.; /* @ invariant 1 <= i <= ni && analytic_error (p , ni ,i -1 , k +1 , a ) */ ❢♦r ( i =1; i < ni ; i ++) { dp = p [ i +1][ k ] -2.* p [ i ][ k ] + p [i -1][ k ]; p [ i ][ k +1] = 2.* p [ i ][ k ] -p [ i ][ k -1] + a * dp ; } p [ ni ][ k +1] = 0.; }
The predicate analytic error states the exact analytical expression of the rounding error. It also states that the rounding error of a single iteration is smaller than a known value. More precisely, it bounds the absolute value of

ε k+1 i := p k+1 i -(2p k i -p k-1 i + A × (p k i+1 -2p k i + p k i-1 )).
Under some hypotheses on A and the ranges of (p k i ), we could prove that |ε k+1 i | ≤ 85 × 2 -52 (and a similar property concerning the initialization of the p 1 i ) [START_REF] Boldo | Floats & Ropes: a case study for formal numerical program verification[END_REF]. The original Coq proof amounts to 735 lines of tactics. Thanks to the gappa tactic, we were able to improve the result: we now have the formal proof that |ε k+1 i | ≤ 80 × 2 -52 ; -drastically cut off the size of the proof script: the 735 lines of tactics reduce to 10. This is a tremendous improvement. Not only is the new proof script dramatically shorter and simpler to write, but it is also more amenable to future changes and maintenance. Indeed, if the program is to be modified in such a way that the error slightly increases, the initial proof would be completely broken and only a small part could be re-used. Using Gappa, the situation is different: while the statement of the theorem would change, the proof would probably be robust enough to remain valid.

Preventing Overflows

Another type of proof that greatly benefits from automation is overflow proofs. Typically, one wants to guarantee that no overflows happen. To do so, it is usually sufficient to bound the program inputs. The resulting VCs are especially tedious to prove. As a consequence, the bounds are often over-estimated in order to simplify the demonstrations. This is the case for the following example. This program computes an accurate discriminant using Kahan's algorithm [START_REF] Kahan | On the Cost of Floating-Point Computation Without Extra-Precise Arithmetic[END_REF]. The accuracy is measured in ulps (unit in the last place), which is the distance between two consecutive floating-point numbers. The discriminant algorithm relies on the exactmult function which computes the rounding error of a multiplication.

/* @ requires xy == round ( x * y ) && @ ( x * y ==0 || 2^^( -969) <= | x * y |) && @ | x | <= 2^^995 && | y | <= 2^^995 && | x * y | <= 2^^1022 @ ensures \ result == x *y -xy @ */ ❞♦✉❜•❡ exactmult ( ❞♦✉❜•❡ x , ❞♦✉❜•❡ y , ❞♦✉❜•❡ xy ); /* @ requires @ ( b ==0 || 2^^( -916) <= | b * b |) && @ ( a * c ==0 || 2^^( -916) <= | a * c |) && @ | b | <= 2^^510 && | a | <= 2^^995 && | c | <= 2^^995 && @ | a * c | <= 2^^1021 @ ensures \ result ==0 || @ |\ result -( b *b -a * c )| <= 2* ulp (\ result ) @ */ ❞♦✉❜•❡ discriminant ( ❞♦✉❜•❡ a , ❞♦✉❜•❡ b , ❞♦✉❜•❡ c ) { ❞♦✉❜•❡ p ,q ,d , dp , dq ; p = b * b ; q = a * c ; ✐❢ ( p + q <= 3* fabs (p -q )) d =p -q ; ❡•s❡ { dp = exactmult (b ,b , p ); dq = exactmult (a ,c , q ); d =( p -q )+( dp -dq ); } r❡t✉r♥ d ; }
The formal proofs for this program (including overflows) have been presented in [START_REF] Boldo | Proof and certification for an accurate discriminant[END_REF][START_REF] Boldo | Kahan's algorithm for a correct discriminant computation at last formally proven[END_REF]. Here we only focus on the overflows of the discriminant; we do not care about the exactmult function.

All the overflow proofs were first done prior to the gappa tactic. For seven proof obligations, it took more than 420 lines of Coq. Using the tactic, the proofs reduce to 35 lines (about 5 lines per theorem). The Coq compilation time, however, is about 5 times greater. 13 Nevertheless, the profit is clear in the verification process as the time for developing the proof overwhelms the time for compiling it.

It is also interesting to note that the specification is also improved. The hypothesis |a × c| ≤ 2 1020 in the original proof [START_REF] Boldo | Kahan's algorithm for a correct discriminant computation at last formally proven[END_REF] was too strong; we proved instead that |a × c| ≤ 2 1021 is sufficient to guarantee that no overflows occur. The proof was not modified at all after changing the annotations. This means that the automation is sufficient to use exactly the same proof when modifying slightly the specification. This is really worthwhile for proof maintenance.

Implementation Details

The gappa tactic is part of the standard V8.2 Coq distribution. 14 It relies on Gappa, which is an external stand-alone tool and comes with its own library of Coq theorems. As usual with a proof assistant, the user issues tactics in order to split the goal into simpler subgoals that can be handled automatically. If one subgoal is in the scope of the Gappa tool, the user can proceed as follows.

First of all, the goal has to be translated to the floating-point model of Gappa. Stored in an auxiliary library, some theorems state that both models are equivalent. In particular, if a floating-point number is the closest to a real number in the Caduceus model, then it is the result of a rounding function in the Gappa model, and reciprocally. The why2gappa tactic automatically applies these theorems to rewrite the goal and its context (Step ②). It also unfolds the Caduceus model of floating-point numbers, with its floating-point, exact, and model parts. At this point, the goal is made of inequalities between real-valued expressions potentially containing rounding operators matching Gappa's ones. Now, the user can launch the Gappa tool to finish the proof, thanks to a single call to the gappa tactic (Steps ③, ④, ⑤, and ⑥). During Step ④, some OCaml code embedded into Coq reads the goal and outputs a text file suitable for Gappa. This code then runs Gappa and asks for a Coq script of the result (Step ⑤). Another OCaml code loads this script into Coq, checks it, and generates the corresponding λ-term, and uses it to finish the proof (Step ⑥).

This process will succeed only if the type of the Gappa λ-term matches exactly the Coq goal of the user. Otherwise, Coq would rightfully complain that the overall proof is not well-typed. For instance, the user goal could mention the inverse x -1 of a variable, while the generated proof would consider 1/x instead. Although equal, these two terms are not convertible, so type-checking would fail. Rather than transforming the generated script afterward, we decided to transform the goal beforehand. The gappa prepare tactic (Step ③), called internally by the gappa tactic, makes sure that the goal will not leave any margin of interpretation to the OCaml code nor to Gappa.

This subtactic is written in Ltac, the tactic language embedded into Coq and available to user scripts. It transforms all the hypotheses and the goal so that they are enclosures of the form m

1 • 2 e1 ≤ expr ≤ m 2 • 2 e2 , with m 1 , e 1 ,
0 ≤ |e + (5 • 2 0 ) × y| ≤ 3 • 2 -3 .
In order to transform the propositions, the tactic could perform some patternmatching to find all the sub-terms that look unadapted and apply rewriting theorems to them. This method is easy to implement but slow, as a huge number of rewriting operation may be needed, especially for constants. (For instance, the real number 11 is implicitly stored by Coq as 1 + (1 + 1) × (1 + (1 + 1) × (1 + 1)).) Instead, the tactic builds an inductive object that represents the syntax tree of the expressions. Some Coq functions (defined in the logic language, not in the tactic language) then implement the previous transformations. We have proved they generate a syntax tree whose evaluation as a real-valued expression gives the same result as the previous expression. Hence applying this single theorem is enough to get a suitable goal. In other words, the tactic simplifies the goal by convertibility and reflexivity [START_REF] Boutin | Using reflection to build efficient and certified decision procedures[END_REF][START_REF] Melquiond | Proving bounds on real-valued functions with computations[END_REF], which is both time-and space-efficient.

Step ④ is then trivial: the OCaml code just has to select the propositions that are enclosures, to visit the nodes of their simplified syntax trees, and to produce the corresponding Gappa script. If Gappa succeeds in verifying the script (Step ⑤), the OCaml code can then load the produced proof and have Coq check it. It takes only a few seconds for the gappa tactic to reach this point after it is called.

We, however, wanted Coq not only to check the generated proof, but also to embed it into the current user λ-term. Therefore, the gappa tactic, while calling an external prover, does produce a complete Coq proof of the goal. Unfortunately, Coq is unable to deal with two scripts at once. So the tactic first launches a separate Coq session that produces a λ-term in the context of Gappa's libraries. Then it runs another separate session to get a λ-term with fully-qualified names and no notations. This last λ-term can finally be loaded in the original user session, without interfering with user-defined names and notations. This incurs a noticeable slowdown for the user. It could be fixed in two ways: enhance Coq so that other scripts can be checked in the same session, or enhance Gappa so that it directly produces a plain λ-term. Embedding the proof script into the user proof hardly increases the confidence though, since the script has already been checked by Coq. So, Step ⑥ could be reduced to type-checking the Gappa proof, creating an axiom with its type, and applying this axiom to the goal. This would ensure that the gappa tactic takes a few seconds only, without having to modify either Coq or Gappa.

Lastly, note that the gappa tactic accesses only a small part of Gappa's features. Indeed, when using the tool directly, the user can pass hints regarding properties of the problem, such as mathematical identities, to guide it. As long as the goals do not need any particular user hint, the tactic is as powerful as the tool.

Conclusion

We have presented an integration of the Gappa automated prover in the Coq proof assistant. This greatly eases the verification process of numerical programs. As shown with realistic examples, the gappa tactic significantly reduces the size of Coq proofs, and improves their maintainability. This tactic is part of the V8.2 Coq standard distribution.

This paper focuses on C programs verified with the Caduceus tool. However, this approach is generic enough to apply to other verification tools, such as Frama-C15 for C programs or Krakatoa for Java programs [START_REF] Filliâtre | The Why/Krakatoa/Caduceus Platform for Deductive Program Verification[END_REF]. Indeed, our work builds upon the Why platform, which provides a common backend for Caduceus, Frama-C/Jessie, and Krakatoa. Said otherwise, any verification technology using Why to produce Coq verification conditions can benefit from the gappa tactic.

The current gappa tactic does not encompass all the features of the Gappa tool. As explained before, there is no way to pass hints to Gappa, such as interval bisection or equalities. Moreover, while the tool can infer enclosures for variables and expressions, the tactic does not offer a way to query them. This feature would relieve the user from the burden of guessing logical cuts.

The Gappa tool is limited to a small logical fragment dedicated to floatingpoint arithmetic, and so is the gappa tactic. A more ambitious perspective is to integrate Gappa to a state-of-the-art SMT solver such as Alt-Ergo [START_REF] Conchon | CC(X): Semantical Combination of Congruence Closure with Solvable Theories[END_REF]. This would result in more VCs discharged automatically but also in more automation when invoked from Coq.

Fig. 1 .

 1 Fig. 1. The Caduceus tool.

6 apply

 6 Rinv_1_lt_contravar ; auto with real .

/

  * @ requires | x | <= 1./32 @ ensures |\ result -cos ( x )| <= 2^^( -23) @ */ ❢•♦❛t toy_cos ( ❢•♦❛t x ) { r❡t✉r♥ 1. f -x * x * .5 f ; } Note that 2 -23 is a tight bound on the error of this function. It ensures that the computed result is one of the floating-point numbers close to the mathematical value cos(x).

2

  intros ; why2gappa ; unfold cos .

3 assert( 4 <= 7 /

 347 Rabs ((1 -( f * f ) * (5/10)) -Rtrigo_def . cos f ) 134217728 )% R 5 by interval with ( i_bisect_diff f ).

Figure 2 Fig. 2 .

 22 Fig. 2. Dataflow in the Coq and Gappa combination

  m 2 , and e 2 explicit integers. Moreover, the expr part should only contain the basic arithmetic operators (+, -, ×, ÷, √ •, and | • |), rounding operators, identifiers, and constants (m • 2 e or m • 10 e ). For instance, if a proposition is | exp(x) + 5 × y| ≤ 3/8, the tactic will generalize exp(x) to a fresh identifier e everywhere. Then it will replace the proposition by the (equivalent yet not convertible) proposition 0 • 2

Available at http://why.lri.fr/.

That means we only establish partial correctness of recursive functions.

For the specification to be complete, the postcondition should also state that the array is a permutation of its initial value. It can be done, but is omitted here for the sake of simplicity.

Available at http://lipforge.ens-lyon.fr/www/pff/.

Available at http://lipforge.ens-lyon.fr/www/gappa/.

The generated proofs depend on a library of facts written for the target system. Currently, this formalization has been proved for Coq only.

The number 0x1.FFFFFEp127 is the biggest finite floating-point numbers for IEEE-754 single-precision format, written with the notation of the standard of the ISO C language (1999).

[START_REF] Bertot | Interactive Theorem Proving and Program Development. Coq'Art: The Calculus of Inductive Constructions[END_REF] The number 1b-21 means 2 -21 , which is almost the optimal upper bound on the specified absolute error.

The Coq script is reproduced verbatim. In particular, some terms are obfuscated due to Coq renaming them to prevent conflicts. So f designates in fact the variable x; and Rtrigo def.cos is the name of the cosine function in Coq's standard library.

Should Coq only check proofs generated by Gappa instead of embedding them, the compilation time would be equivalent. See end of Section 4.

Available at http://coq.inria.fr/.

Available at http://www.frama-c.com/.

⋆ This research was partially supported by projects ANR-05-BLAN-0281-04 "Cer-PAN" and ANR-08-BLAN-0246-01 "F ST". 3 Note that π/2 cannot be represented by a floating-point number, therefore cos(x) cannot be zero.