domain: spi phys

text:
Tomasz Matthia

Baptiste Fix
email: baptiste.fix@onera.fr

Léna Soun

Christophe Dupuis

Nathalie Bardou

Patrick Bouchon

High efficiency second harmonic generation in coupled nano Fabry-Perot thin resonators

In this paper we experimentally demonstrate SHG enhancement in thin 1D periodic plasmonic nanostructures on GaAs in the infrared spectral range. Due to the properly designed coupling of horizontal Fabry-Perot nanoresonators that occur inside these structures, the obtained conversion efficiencies go up to the 10 -7 W -1 range. Moreover, we demonstrate that the engineering of the plasmonic nanoantennas dimensions on the same GaAs layer can lead to SHG enhancement for pump wavelength ranging from 2.8 µm to 3.3 µm.

Thin nanostructured metasurfaces have the property of trapping and locally concentrating the incoming light in a subwavelength volume [1]. This ability to highly enhance local optical fields is particularly interesting for exciting naturally weak nonlinear optical effects. Moreover, these structures tackle some of the most challenging issues of conventional devices based on macroscopic bulk nonlinear media such as the need of phasematching, the limited operational wavelength ranges and the size-related challenges [2][3][4][5][6][7]. Second order generation (SHG) has been demonstrated with all-dielectric devices based on nonlinear materials nanostructured so as to form nanocylinders [8][9][10][11], gratings [12][13][14] or waveguides [15]. These all-dielectric SHG metasurfaces have demonstrated conversion efficiencies going as high as Π SHG = 10 -5 W -1 for a conversion from the telecom to visible wavelength [16]. Such efficiencies have not yet been reached in the mid-infrared. Metallic metamaterials have also been investigated for SHG and among those, two categories can be distinguished. First, those that exploit the centrosymmetric breakdown at the metal surface, where the localized surface plasmon resonance induces strong local field enhancement [17][18][19][20][21][22][23]. However, the efficiency of those devices is limited by the ohmic losses occurring in the metal [3]. Second, the use of plasmonic nanostructures that trap and locally concentrate the optical fields inside the volume of a highly nonlinear layer [24][25][26][27]. This associates the advantages of both worlds, the intrinsic nonlinearities of high χ (2) bulk media with the high field enhancement provided by plasmonic devices. Multiple quantum wells has been integrated to such resonators and demonstrated to reach efficiencies as high as 10 -2 W -1 [28][29][30][31]. However, these devices are inherently spectrally limited to the quantum well intersubband transitions. On the other hand, even if it leads to weaker efficiencies, exploiting the intrinsic nonlinearities of a high χ (2) bulk media included in a plasmonic resonator enables the conversion for any spectral range. For instance, metal-insulator-metal (MIM) nanoresonators, supporting Fabry-Perot resonances and including a gallium arsenide (GaAs) layer, have been demonstrated to reach efficiencies as high as 10 -9 W -1 in the mid-infrared [26]. In this Letter, we investigate coupled GaAs plasmonic nanocavities and experimentally demonstrate a high SHG efficiency in the range of 10 -7 W -1 for a mid-wave infrared pump. We show that coupled nano Fabry-Perot resonances enable a 20-fold increase in SHG efficiency with regards to their classical single cavity counterpart. Finally, we experimentally demonstrate SHG on three arrays fabricated on the same GaAs layer for pump wavelength ranging from 2.8 µm to 3.3 µm which evidences the tunability of the resonators.

The studied periodic nanostructure (see Fig. 1 a)), is composed of three layers: a back gold mirror, a thin (h = 142 nm) layer of (100)-oriented GaAs crystal whose lattice vectors are aligned along (xyz) and a thin (t = 60 nm) top periodic nanostructured gold mirror composed of two equally spaced ribbons. It has been simulated with the B-spline modal method (BMM) [32,33] as infinite along y and infinitely periodic along x. The model of gold used here is a Drude model in agreement with experimental data in the infrared range: ϵ Au (λ) = 1 -1/(λ P /λ(λ P /λ + iγ)) with λ P = 159 nm and γ = 0.0075 [34]. The nanostructure, defined by a period d = 1010 nm, and two ribbons of respective widths w 1 = 260 nm and w 2 = 380 nm, has been optimized to display a resonance at both the pump wavelength λ pump = 2.93 µm and the signal wavelength λ signal = 1.46 µm. Indeed, as it has been previously described [25,35], the two resonances are used to enhance the nonlinear process of light generation. The measured second-harmonic signal power as a function of the incident pump power. The red line is a least-square regression corresponding to the following mathematical expression: log(P SH ) = 2 • log(P inc ) -1.138. c) The second-harmonic signal power as a function of the incident pump polarization angle.

Moreover, the GaAs is a zync-blende crystal whose only nonzero components of the second-order nonlinear susceptibility tensor are [36]:

χ ijk =    150 pm V ⇔ i ̸ = j ̸ = k 0, otherwise (1)
Where i, j, k ∈ {x, y, z}. Thus, the whole SHG process can be mathematically summarized in the form of the following expression:

E i 2ω ∝ ∑ jk χ ijk (2) • E j ω • E k ω ( 2
)
Where E i ω is the i-axis directed electric field at the frequency of ω. With the 1D periodic structure described, the only solution available consists in a TM-polarized (along the x-axis) pump wave and TE-polarized signal wave (along the y-axis).

The structure simulated and measured SHG efficiencies are presented in Fig. 1 c). It has been calculated using the common linearized formula :

Π SHG = P SH P inc 2 W -1 (3)
Where P SH and P inc are respectively the power of the generated second-harmonic beam and the power of the incident pump beam. The peak efficiency of 1.69 × 10 -7 W -1 is achieved at 2.93 µm which represents a 750-fold increase with regard to an equivalent unstructured device.

The optimization of the SHG requires to efficiently trap the pump light inside the GaAs layer and to extract the signal light. This means that, at the two resonances, the cavity modes must be critically coupled to the incident light : the reflectivity of the structure has to be minimized [26]. Among the different available nanostructures, the combination of TE orientated vertical and TM orientated horizontal Fabry-Perot (FP) resonators are well suited for SHG due to their stability with regards to the angle of incidence and the localisation of the concentrated fields. However, the two resonators critical coupling are not compatible due to the high impact of the GaAs layer thickness on both these resonators. To achieve the critical coupling of both resonances, a newly described type of horizontal FP resonance has been used [37]. Indeed, this "coupled Fabry-perot" resonance is based on two spectrally closed under-coupled FP (i.e., with high reflectivity at resonance) that couple with one another. While most properties of the horizontal FP are retained, the resulting resonance is mostly independant from the GaAs layer thickness. The critical coupling is achieved through the optimization of the width difference between the two ribbons.

The impact of using either horizontal FP or coupled FP on the SHG efficiency is illustrated in Fig. 2. On the one hand, the vertical Fabry-Perot resonance is only slightly influenced by the widening of the second ribbons from w 1 = 260 nm to w 2 = 380 nm. On the other hand, the structure reflectivity at the pump wavelength is drastically reduced from a 75% to 5% leading to a 20-fold increase of the simulated SHG efficiency. Those theoretical results are also congruent with the comparison of the experimental results obtained here (peak of 1.69 × 10 -7 W -1 ) and in analogue single Fabry-Perot resonator structures (peak of 0.8 × 10 -8 W -1 ) [26].

The sample was fabricated first by epitaxial growing of a (100)-oriented layer of GaAs, and depositing it on a Ti/Au/Ti/Al stack in order to create the golden mirror. After mechanical polishing and chemical etching, the golden nanoantennas were deposited by an electron beam lithography process. The structure includes respectively: a golden substrate, a 2 nm thick titanium layer, a 142 nm thick GaAs nonlinear layer, a layer of Al 2 O 3 of 10 nm, a grating of 4 nm titanium adhesion layer placed under 60 nm golden nanoantennas. The sample contains several 300 µm width square shaped structures.

The TE and TM linear reflectivity spectrum (see Fig. 3) of the structure have been measured using a FTIR microscope (Bruker Hyperion 2000) and a Cassegrain objective (NA=0.4) with high accordance to the simulated values. The simulated normalized electric field components taking part in the SHG process at both λ pump = 2.92 µm (TM) and λ signal = 1.46 µm (TE) wavelength are shown in Fig. 3. The non linear characterisation set-up is illustrated in Fig. 4 a). The incident pump light is provided by the idler output of an APE nanoLevante optical parametric oscillator (OPO). The laser beam repetition rate is equal to 20 kHz and the peak power density is equal to 0.58 MW cm 2 . The beam direction of polarization is adjusted to be TM, with respect to the sample geometry, by a half-wave plate.

An additional polarizer can be introduced before the halfwave plate to regulate the incident optical power. The beam is focalized on the sample by a 5 cm CaF 2 focal lens, resulting in a focal spot of 37 µm diameter and an aperture angle of 9 degrees. The generated second-harmonic beam is collimated by the same lens and is reflected by dichroic plate (Ge) towards a Xenics Xeva-1.7-320 SWIR camera. The incident pump light is monitored by measuring the 10% pump power reflected by the Ge plate. Due to its thickness, the germanium plate also introduces a walk-off of the pump beam with regards to the center of the focusing lens, resulting in an mean incident angle of about 4.5°. The evolution of the second harmonic beam generation with the incident pump power and the direction of polarization is presented in Fig. 4. The SH generation power follow the expected quadratic dependency with the incident pump. Moreover, the SHG process depends on the TM polarisation of the pump which is confirmed by the relation of the SH power with the direction of polarisation.

Noteworthily, a wide variety of resonance wavelength can be achieved for both cavities while the nonlinear material thickness is maintained constant. Indeed, as illustrated in Fig. 5, the nonlinear generation can be spectrally shifted by varying the period of the structure as well as the width of each of the plasmonic nanorods. The first resonator (green line) is the one that has been presented in Fig. 1. The two additional SHG structures have been obtained for w 1 = 280 nm, w 2 = 365 nm, d = 1015 nm (purple line) and w 1 = 295 nm, w 2 = 350 nm, d = 1215 nm (blue line) to have a peak generation for a pump at 2.8 µm and 3.24 µm respectively. The two structures have been characterised through the same set-up.

In conclusion, we have theoretically and experimentally studied the properties of multiresonant MIM structures optimized for the SHG. Due to an increased pump light absorption, the coupling of two horizontal Fabry-Perot nano-resonators provides a 20-fold increase in the obtained efficiencies. This increase enables us to attain the 10 -7 W -1 efficiency range in bulk GaAs for wavelength ranging from λ = 2.8 µm to λ = 3.3 µm. The measured efficiencies are in good agreement with the theoretical predictions. This advancement paves the way to more efficient nonlinear plasmonic nanodevices. Moreover, those devices can be spectrally adjusted without modifying the substrate and the GaAs layer (see Fig. 5 and the supplement). This gives the possibility, either of creating discrete arrays of different resonators, or combining them into a single wide-spectrum converter. Different approaches may be studied in order to further increase the SHG efficiency. First, the use of a different set of resonances could improve the spatial mode overlap. Finally, the interaction between the electric fields and the nonlinear layer could be enhanced by better fitting the components of the nonlinear susceptibility tensor or by using multiple quantum wells.

Funding. This work is supported by a DGA -Agence Innovation Défense scholarship.

Disclosures.

The authors declare no conflicts of interest.

Fig. 1 .

 1 Fig. 1. a) Schematic representation of the resonator. b) A scanning microscope image of the sample. c) The measured and simulated SHG efficiency spectra of the coupled resonator nanostructure. The resulting efficiency is 750 times higher compared to a non-engineered equivalent GaAs layer reference.

Fig. 2 .

 2 Fig. 2. a) Schematic explanation of the single and coupled FP resonators behavior for the SHG. b) Simulated linear reflectivity of the single resonator (red dashed) and coupled resonators (plain blue) structure at TE-polarized signal frequency. c) Simulated linear reflectivity of the single resonator (red dashed) and coupled resonators (plain blue) structure at TM-polarized pump frequency. d) Simulated nonlinear SHG efficiency of the single resonator (red dashed) and coupled resonators (plain blue) structure at pump frequency.

Fig. 3 .

 3 Fig. 3. Simulated and measured linear reflectivities of the structure. In inset the simulated normalized field component amplitudes in the structure at pump (λ pump = 2.93 µm) and signal (λ signal = 1.46 µm) resonances.

Fig. 4 .

 4 Fig. 4. a) Schematic of the non linear characterisation set-up b)The measured second-harmonic signal power as a function of the incident pump power. The red line is a least-square regression corresponding to the following mathematical expression: log(P SH ) = 2 • log(P inc ) -1.138. c) The second-harmonic signal power as a function of the incident pump polarization angle.

Fig. 5 .

 5 Fig. 5. a) An illustration of the concept of a broad-band SHG nanodevice. b) Measured and simulated SHG conversion efficiency spectra obtained for w 1 = 295 nm, w 2 = 350 nm, d = 1215 nm (blue line), w 1 = 280 nm, w 2 = 365 nm, d = 1015 nm line) and the structure from Fig. 1 c) (green line). All three resonators have been designed using the same thicknesses of the golden nanoantennas t = 60 nm and the GaAs nonlinear layer h = 142 nm.

Data Availability Statement. All numerical and experimental data

are available on request from the authors.

Supplemental document.

See Supplement 1 for supporting content.


Philippe Chevalier

Adrien Moreau

Francis Bessière

Sylvain Richard

Mohamed Chahine

Gilles Millat

Elodie Morel

Franck Paganelli

Nathalie Lesavre

Leslie Placide

Identification of Cx43 variants predisposing to ventricular fibrillation in the acute phase of ST-elevation myocardial infarction

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


Ali Ayad

On formal methods for certifying floating-point C programs

Keywords: static analysis of programs, formal methods, deductive verification, automatic theorem provers, interactive proof assistants, floating-point arithmetic, IEEE-754 standard

. When the Why proof obligations are written in the syntax of the Gappa library, we can use the gappa [3] and interval [15] tactics to achieve the proof. Several examples of floating-point C programs are presented in the paper to prove the efficiently of this implementation.

Sur les méthodes formelles de certification de programmes C flottants

Résumé : Ce papier présente une implémentation d'une extension du langage de spécification ACSL dans l'outil Frama-C dans le but de prouver la certitude de programmes C flottants. Cette implémentation est essentiellement basée sur les deux modèles Why de l'arithmétique flottante de [START_REF] Boldo | Formal Verification of Floating-Point Programs[END_REF]: le premier modèle suppose qu'il n'y a pas de dépassement de mémoire (overflow) durant l'exécution du programme et est appelé le modèle Real, i.e., les nombres flottants sont supposés de réels et les opérations flottantes sont faites avec une précision infinie. Le deuxième modèle vérifie qu'il n'y a pas d'overflows, i.e., des obligations de preuve sont générées par l'outil Why pour prouver que le résultat de l'opération flottante n'est pas plus grand que le flottant maximal du format de destination, ce modèle est appelé le modèle Strict. Dans ce papier, on décrit un troisième modèle, appelé le modèle Full, qui étend le modèle Strict. Le modèle Full permet d'avoir d'overflows et traite de valeurs spéciales: infinis signés, NaNs (Not-a-Number) et zéros signés comme dans la norme IEEE-754 (voir [2,[START_REF] Goldberg | What every computer scientist should know about floating-point arithmetic[END_REF]).

Les conditions de vérification générées par Why sont (partiellement) prouvées par de prouveurs automatiques: Alt-Ergo, Simplify, Yices, Z3, CVC3 et Gappa ou par l'assistant de preuve Coq [START_REF] Coq | The Coq Proof Assistant Reference Manual -Version V8.2[END_REF] en utilisant deux formalisations existantes de l'arithmétique flottante dans Coq: la première vient de l'article [START_REF] Théry | A generic library for floating-point numbers and its application to exact computing[END_REF] et la deuxième est la bibliothèque Gappa [START_REF] Melquiond | Floating-point arithmetic in the Coq system[END_REF]. Quand les obligations de preuve de Why sont écrites dans les syntaxes de la bibliothèque Gappa, on pourrait utiliser les deux tactiques Coq gappa [3] et interval [START_REF] Melquiond | Interval package for coq[END_REF] pour finir la preuve. Plusieurs exemples de programmes C flottants sont présentés dans ce papier pour prouver l'efficacité de cette implémentation.

Mots-clés : analyse statique de programmes, méthodes formelles, vérification déductive, prouveurs automatiques de théorèmes, assistants de preuve interactifs, arithmétique flottante, norme IEEE-754

Introduction

Many critical problems such as numerical analysis, physics, aeronautics (software embedded in satellites, robots), energy (nuclear centers), automotive, etc. deal with floating-point computations. As we rely more and more on softwares, the consequences of a bug are more and more dramatic, causing great financial and even human losses. We cite two historical accidents caused by bugs in softwares, the first one is the overflow bug that caused the failure of the Ariane 5 launcher in 1996: a double precision number which represents the horizontal speed of the launcher was converted into a 16 precision format, this speed was greater than the largest integer representable in 16 bits, the conversion failed without triggering an exception. The second one is the cumulated imprecision errors in a Patriot missile defense that caused it to miss its Scud missile target, resulting in 28 people being killed and around 100 injured in 1992 : the error was an inaccuracy in computing the date of the Patriot anti-missile. Thus for these reasons, one needs good tools to ensure the proper functioning of critical softwares.

There are many techniques to ensure that a program runs well every time we give it new input values and does not create errors at run-time such as division by zero and out-of-bound memory access. One can try to test all possible values of the input variables and verify that they have the correct outputs by the program. This can be done if the number of these tests is small and their running time is reasonable. Dynamic analysis performs tests to study important behaviors of the program. Dynamic analysers try to automatically generate path tests for programs and to extract properties that are true over all executions of the program. However, static analyse performs verification of programs without executing them by applying mathematical techniques such as first order logic and abstract interpretation.

The aim of this paper is to apply mathematical formal methods and theorem provers to static analysis of programs that handle floating-point arithmetic.

We are interested in certifying imperative programs, especially C programs. We will be able to prove the correctness of such a program, i.e., to ensure that the program is free of threats (null pointer dereferencing or out-of-bounds array access, etc ...) and to check that there is no overflow, never NaNs (Not-a-Number) when dealing with floating-point programs. Additionaly, one can prove some functional properties introduced manually as annotations in the source code. The general strategy is based on a classical Hoare logic framework [START_REF] Hoare | An axiomatic basis for computer programming[END_REF] and follows the following steps:

Formal specification We annotate the program according to the properties we would to preserve after execution, i.e., preconditions, postconditions, (global, loop) invariants, termination, etc ... that determine the memory state of the imput and the output variables of the program. These annotations are generaly written in a programming language inspired from ML language. In this paper we use the ACSL language (ANSI C Specification Language [4]) which is partially implemented in the Frama-C verification tool [START_REF] Ayad | Framework for modular analyses of C[END_REF]. Generation of verification conditions (i.e., logical formulas of first order that we should prove their validity). Different verification tools are intended to generate verification conditions from annotated C programs, e.g., the Why platform [START_REF] Filliâtre | The Why verification tool[END_REF]8] and its verification tool Caduceus [START_REF] Filliâtre | Multi-prover verification of C programs[END_REF] for C programs. In this paper, we will apply the Frama-C tool to an annotated C program which transforms it to an annotated program written in the Jessie language (an intermediate language for the verification of Java and C programs [START_REF] Marché | Jessie: an intermediate language for Java and C verification[END_REF]). The Jessie plugin (see e.g., [START_REF] Yannick | Automatic Modular Static Safety Checking for C Programs[END_REF][START_REF] Yannick | Jessie plugin tutorial[END_REF] produces then the associated annotated Why file (with the syntaxes of the Why language). Finally, the Why platform generates the verification conditions that reflect the behavior of the C program specified by its annotations. All these transformations preserve the semantics of the annotated program. Note that Why does not interprete the annotations: they are built from functions and predicates axiomatization or defined in models that can be specific to each prover.

Formal verification These verification conditions are sent to a proof assistant in order to prove them. There are two kind of proof assistants: automated theorem provers such as Simplify, Alt-Ergo, haRVey, Yices, Gappa, etc ... and interactive proof assistants such as Coq, PVS, Isabelle/HOL, etc ... This paper is organized as follows. Section 2 gives an overview of the floatingpoint arithmetic in the IEEE-754 standard [2]. In section 3, we describe the Frama-C compilation chain, i.e., the ACSL, Jessie and Why languages. Section 4 presents an improved version of the Strict model of [START_REF] Boldo | Formal Verification of Floating-Point Programs[END_REF] and a complete description of the Full Why model. A Coq formalization of these two models with an illustration of their efficiency on some examples are described also in Section 4. We terminate the paper by a conclusion with some perspectives in Section 5. Different tools and libraries are necessary to run correctly our examples, they are cited in a technical appendix. where s = 0 or s = 1 (to determine the sign of x), m ∈ N is the significand and e ∈ Z is the exponent of x.

Definition 1 A binary format f is a tuple (p, e min , e max ) where p is a positive integer called the precision of f and e min and e max are two integers which define a range of exponents for f . Definition 2 A binary floating-point number x = (-1) s × m × 2 e is representable in a binary format f = (p, e min , e max ) if m and e verify the following inequalities : 0 ≤ m < 2 p and e min ≤ e ≤ e max .

Of course, there is a finite number of binary floating-points representable in a given binary format f , i.e., f defines a finite subset of R. The representation of a binary floating number x in f is not unique. For example, for f = (3, -5, 6), x = 6 × 2 0 = 3 × 2 has two different binary representations in f . In general, real numbers with denominator not a power of 2 are not representable in any binary formats, they are approximated by a representable binary floatingpoint number w.r.t. a certain rounding direction (see below). For example, for f = (24, -149, +104), 0.1 is approximated by the binary floating-point number 13421773 × 2 -27 which is representable in the format f .

The IEEE-754 standard

The standard IEEE-754-2008 [2] defines methods to perform binary and decimal floating-point arithmetic. It facilitates the portability of programs between computers which support the standard i.e., programs will product the same results of numerical computations on different architectures that implement the standard. The standard specifies binary and decimal formats to represent floating-point numbers with the operations of conversion between these formats. It specifies the elementary operations (addition, subtraction, multiplication, division, square root, Fused-Multiply-Add, remainder, etc ...) and the comparison operators on floating-point numbers. It specifies conversion operators between different formats and between integers and floating-point numbers. It handles floating-point exceptions, including data that are not numbers (NaNs).

Formats

The standard has introduced three basic binary formats: binary32 (Single), binary64 (Double) and binary128 (Quad) which are defined with the above notations by: ❼ binary32 = (24, -149, +104)

❼ binary64 = (53, -1074, +971) RR n➦ 6927 ❼ binary128 = (113, -16494, +16271)
For example, a real number representable in binary32 is encoded by a string of bits of length 32 (1 bit for the sign, 23 bits to represent the significand and 8 bits for the exponent). The standard defines also two decimal formats decimal64 and decimal128, extended and extendable binary and decimal formats.

Normal and denormal binary floating-point numbers

Let f = (p, e min , e max ) be a binary format and x = (-1) s × m × 2 e be a binary floating-point number which is representable in f . In this subsection, we define normal, denormal and canonic binary floating-point numbers as follows:

Definition 3 The binary floating number x is said to be normalized in the format f if 2 × m ≥ 2 p , this is equivalent to that the leading bit in the binary representation of the significand m is nonzero, i.e., m is represented by exactly p bits.

Definition 4

The binary floating number x is said to be denormalized in the format f if the exponent e of x is minimal and 2 × m < 2 p , this is equivalent to that e = e min and the leading bit in the binary representation of the significand m is zero, i.e, m is represented by fewer than p bits.

Definition 5 The binary floating number x is said to be canonic in the format

f if it is normal or denormal in f .
The smallest positive normalized floating-point number in the format f is 2 emin+p-1 and the largest is 2 emax (2 p -1). The smallest positive denormalized floating number in the format f is 2 emin . In the standard, zero is neither normalized nor denormalized and there are distinct representations for +0 and -0 for all formats but +0 = -0.

Rounding direction and basic operations

The standard requires that every floating-point operation shall be performed as if it first produced an intermediate result correct to infinite precision and with unbounded range of exponent, and then rounded that result according to the correspond rounding direction. Rounding a real number x (i.e., a result of an operation with infinite precision) modifies it to fit in the destination's format f of the operation. The standard has defined five rounding directions as follows: one can always find two representable floating-point numbers x 1 and x 2 in the format f such that

x 1 ≤ x ≤ x 2 (if x is representable in f then x 1 = x 2 = x).
1. roundTiesToEven (nearest_to_even): the rounding of x w.r.t. this direction is the nearest floating-point number to x among x 1 and x 2 . If x is in the middle of the interval [x 1 , x 2 ] then one chooses whose significand is even.

2. roundTowardPositive (up): the rounding of x is x 2 .

3. roundTowardNegative (down): it is the opposite of roundToward-Positive, i.e., the rounding of x is x 1 .

INRIA 4. roundTowardZero (to_zero): the rounding of x w.r.t. this direction is

x 1 if x > 0 and x 2 if x < 0.
5. roundTiesToAway (nearest_away): the rounding of x is the closest to x among x 1 and x 2 . If x is exactly the middle of the interval [x 1 , x 2 ] then the rounding is

x 2 if x > 0 and x 1 if x < 0.
If the infinitely precise sum (or minus) of two floating-point numbers is zero, then the result is -0 if the rounding direction is roundTowardNegative and +0 otherwise. The standard classified floating-point operations in two categories:

❼ Homogeneous operation: the operands and the result have the same format in the same base.

❼ Non-homogeneous operation: the operands and the result have different formats in the same base.

Special values and exceptions

The standard defined three special values: -∞, +∞ and N aN (Not-a-Number). These numbers should be treated both in the input and the output of the arithmetic operations as usually. For example, (+∞)

+ (+∞) = (+∞), 1/(+∞) = +0, 1/(-0) = -∞, (+0)/(+0) = N aN , (-0)×(+∞) = N aN and √ -1 = N aN but √ -0 = -0.
Operations propagate NaN operands to their results, then NaNs are not rounded. In addition, the standard handles five exceptions:

Invalid operation exception It holds when the operation is not defined on the operands. For example, the square root of negative numbers and the division of two zeros. The result of an invalid operation is NaN.

Inexact operation exception It holds when the result of the operation is not representable in the destination's format. This result must be rounded.

Division by zero exception It happens in the case of division of a nonzero and non NaN number by zero. The result of this division is infinity with sign the product of the signs of the operands.

Overflow Overflow is said to occur when the exact result (i.e., when the exponent range unbounded) of a floating-point operation is finite but with an absolute value that is larger than the maximal floating-point number represented in the destination's format. As with division by zero, in the days before IEEE arithmetic was available the usual treatment of overflow was to set the result to (±) the largest floating-point number in the destination's format or to interrupt or terminate the program. The IEEE-754 standard response to overflow is to deliver the correctly rounded result, either (±) the maximal floating-point number or ±∞. The result of an overflow operation depends on the rounding direction and the sign of the intermediate result as follows: in the case of roundTiesTo-Even and roundTiesToAway, this is the same as saying that overflow occurs when an exact finite result is rounded to ±∞ according to the sign of the result, but it is not the same for the other rounding modes. If the rounding direction is roundTowardPositive, the result is +∞ (if it is positive) or the opposite of the largest normalized floating-point number of the destination's format (if it is negative). For roundTowardNegative, the result is the largest normalized floating-point number of the destination's format (if it is positive) or -∞ (if it is negative). roundTowardZero is analogue to roundTowardNegative if the result is positive and to roundTowardPositive if the result is negative.

Underflow Underflow is said to occur when the exact result of an operation is nonzero but with an absolute value that is smaller than the smallest normalized floating-point number represented in the destination's format. In the days before IEEE arithmetic, the response to underflow was typically, though not always, flush to zero, i.e., return the result 0. The IEEE-754 standard response to underflow is to return the correctly rounded value, which may be a denormal number, ±0 or (±) the minimal denormal number. This is known as gradual underflow.

Rounding errors

One of the major problems encountered in floating-point computations is the propagation of rounding errors. Because of the fact that not all floating-point operations are exacts, the final result of a chain of floating-point instructions can be very far from the exact result (i.e., the expected result if the computation is done on real numbers). There are three concepts of measurement to bound the errors in operations: the absolute error, the relative error and the ulp error (unit in the last place). If x = (-1) s ×m×2 e is the result of a floatingpoint operation rounded in the destination's binary format f = (p, e min , e max ) of the operation w.r.t. a certain rounding direction.

❼ The absolute error of x is the distance between x and the exact result of the operation (i.e., without rounding).

❼ The relative error ǫ of x is the division of the absolute error by the exact result y of the operation, i.e., x = (1 + ǫ)y. One can bound this error by |ǫ| ≤ 2 -p if the rounding direction is roundTiesToEven and by |ǫ| ≤ 2 1-p otherwise.

❼ ulp(x) is the weight of the last bit in the significand m of x, i.e., ulp(x) = 2 p . For a real r, ulp(r) = 2 ⌊log2|r|⌋+1-p . ACSL is a specification language of C programs. It is inspired from the specification language of the Caduceus tool and the Java Modeling Language (JML).

A part of ACSL is implemented in the Frama-C framework, another part is under experimentation. In this paper, we extend the syntaxes of ACSL to handle floating-point computations. Annotations in C programs are written in comments that start by /*@ or //@ and end as in C comments. In the sequel, we show only a brief overview of the ACSL language, people interested in details can return to the tutorial [START_REF] Prevosto | ACSL mini-tutorial[END_REF] or the reference manual [4].

Types in ACSL

❼ The types integer, real and boolean represent respectively mathematical integers, real numbers and booleans (with just two values \true and \false).

❼ Users can introduce new logic types in ACSL.

❼ C integral types char, short, int and long, signed or unsigned, are all subtypes of the type integer which is itself a subtype of the real type.

Annotations of C functions

There are built-in constructs that evaluate the pre-state and the post-state of a C function. The keywords \old(e) and \result are reserved in ACSL to denote the value of the expression e before the execution of the function and the value returned by the function respectively. Simple contracts are written immediately before the function field.

❼ For a set P of predicates, the contract requires P; means that the predicates of P are satisfied by the input variables of the function.

❼ For a set L of memory locations, the contract assigns L; enforces that no memory locations outside L can be modified. It does not enforce modifications of locations in L.

❼ For a set Q of predicates, the contract ensures Q; asserts that the predicates of Q hold on the output variables of the function.

Annotations of C loops

They are written immediately before the C loops for, while, etc ... ❼ For a set I of inductive predicates, the contract loop invariant I; means that the predicates of I are satisfied at every loop iteration.

❼ For a set L of memory locations, the contract loop assigns L; has the same meaning as assigns L; for every loop iteration.

❼ For an integer expression e, the contract loop variant e; ensures that the loop terminates after a finite number of steps if e decreases during the loop iterations

Global properties

❼ Global invariants: for a set P of predicates, the contract global invariant P; means that the predicates of P are always true on global variables of the program.

❼ Type invariants: for a type τ and a set P of predicates defined over τ , the contract type invariant P; means that the predicates of P are satisfied by any variable of type τ .

❼ Logic specifications: users of ACSL are able to introduce new logic types, logic definitions, functions, predicates and axioms.

The Jessie intermediate language

Jessie language is an intermediate language presented in [START_REF] Marché | Jessie: an intermediate language for Java and C verification[END_REF]. It provides primitive types such as integers, booleans, reals, etc ... and abstract datatypes.

Jessie has not a notion of arrays, structures, pointers and memory heap. Jessie programs can be annotated using pre-and postconditions, loop invariants and intermediate assertions. Logical formulas are written in a typed first-order logic with built-in equality, booleans, infinite precision integers and real numbers with usual operations. We have extended this language by intoducing two new types, float and double as in C programs to deal with floating-point computations.

Ordinary operations and comparison over float and double variables are defined. We have introduced rounding modes and casting operators to Jessie. Special predicates are devoted to test if a float or double variable is finite, NaN or infinity.

The Why language

The platform Why aims at being a verification conditions generator (VCG) back-end for other verification tools. It provides a powerful input language including higher-order functions, polymorphism, references, arrays and exceptions.

The Why language is inspired from the Ocaml [START_REF] Leroy | The Objective Caml system release 3.01, Documentation and user's manual[END_REF] language for functionnal programs. So there is no distinction between instructions and expressions. Special keywords are reserved in Why, for example, result denotes the value returned by a function and is used only in the postcondition. Programs in Why language can be annotated by preconditions, assertions, postconditions, loop invariant and loop variant to ensure the termination. The language contains basic built-in types: bool, unit, int and real. The expressions in Why are formed by: ❼ The integer, boolean constants. The Why tool implements this programming language. It takes annotated programs as input and generates proof obligations for a wide set of interactive proof assistants (Coq, PVS, Isabelle/HOL, HOL 4, HOL-Light, Mizar) and decision procedures (Simplify, Ergo, Yices, CVC Lite, CVC3, haRVey, Zenon). It can verify algorithms rather than programs, since it implements a rather abstract and idealistic programming language. Several non-trivial algorithms have already been verified using the Why tool, such as the Knuth-Morris-Pratt string searching algorithm for instance. It computes Dijkstra's weakest preconditions in the Floyd-Hoare logic. The Why tool can be used to write axiomatizations and goals and to dispatch them to several existing provers. It is independent of the back-end prover that will be used, it makes no assumption regarding the logic used. It uses a syntax of polymorphic first-order predicate logic for annotations with no particular interpretation. Function symbols and predicates can be declared in order to be used in annotations, together with axioms, and they may be given definitions on the prover side later, if needed. For more details on the Why language and tool, you can refer to the reference manual [START_REF] Filliâtre | The Why verification tool[END_REF]8]. The extension to handle floating-point programs is given in the next section.

Why models for floating-point C programs

In this section, we describe an improved version of the Strict Why model of [START_REF] Boldo | Formal Verification of Floating-Point Programs[END_REF] where any operation on floating-point numbers produces a proof obligation to ensure that there is no overflow. We introduce also the Full Why model, i.e., the complete model which implements the IEEE-754 standard, e.g., infinities, signed zeros and NaNs are supposed to be special values in the program. We give an implementation of these two models in the Coq proof assistant with the syntax of the Gappa librabry. Some examples will be treated at the end of this section.

The Strict model

This model is based on the Why theory and written in the syntaxe of the Why language. In [START_REF] Boldo | Formal Verification of Floating-Point Programs[END_REF], a floating-point value is viewed as three parts:

1. the floating-point number of the given type, as computed by the program; 2. a real number, which would be the value if all previous computations were performed on real numbers and thus exact;

3. another real number, which is the value which should be ideally computed.

We define new types, constants, functions and predicates to be able to write annotated floating-point Why programs. The complete description of this model can be found in the next Why distribution in the file lib/why/floats_strict.why.

In the logic:

❼ An abstract type gen_float for generic floats.

❼ A type float_format with three elements: Single, Double and Quad.

❼ A type mode for the five rounding modes of the revision of the IEEE-754 standard: nearest_even, to_zero, up, down, nearest_away as described in Section 2.

❼ Functions float_value, exact_value, model_value to access respectively to the floating-point, the exact and the model part of variables of type gen_float.

❼ A rounding operator gen_float_of_real_logic which rounds a real to a gen_float according to a certain float format and rounding direction. In the program:

❼ A rounding operator gen_float_of_real which rounds a real to a gen_float according to a certain float format and rounding direction.

❼ Basic binary and unary operations: add_gen_float, sub_gen_float, mul_gen_float, div_gen_float, sqrt_gen_float, neg_gen_float, abs_gen_float on the gen_float type.

❼ Basic operators for comparison: lt_gen_float,gt_gen_float, etc ... ❼ A casting operator cast_gen_float between float formats according to rounding directions.

In this model, pre-and postconditions are written for each basic operation of the program. These conditions product proof obligations for testing that there is no overflow at run-time. For example, multiplication is specialized by: All the above definitions are passed to the ACSL language to be able to write annotated C programs. They are defined as follows:

❼ Constants: floating-point constants are interpreted as mathematical real numbers in the specifications, they can be written in decimal or in hexadecimal format. Integers are also seen as mathematical real numbers if necessary. Variables of type float or double in the C program are implicitly converted to the real type when used in the annotations.

❼ Casts: a built-in logic type is defined by the five rounding_mode of the IEEE-754 standard: \NearestEven, \ToZero, \Up, \Down, \NearestAway.

Casting from a C integer type or a float type to a float or a double is as in C. Conversion of real numbers to float or double values is given by the following two logic functions: logic float \round float(rounding mode m, real x); logic double \round double(rounding mode m, real x);

Cast operators (float) and (double) applied to a mathematical integer or real number x are equivalent to apply the above rounding functions with the \NearestEven rounding_mode (which is the default rounding mode in C programs). Suffixes f and l are meaningful, because they implicitly add a cast operator as above.

❼ Floating-point operations: in the Strict model, the sum (resp. minus, product and division) of two C floating-point variables x and y is equivalent to the sum (resp. minus, product and division) of their associated real numbers with the mathematical sense.

❼ A special predicate to verify if there is an overflow in floating-point operations is declared by: predicate \no overflow(f:float format,rounding mode m,real x);

where float_format is a built-in logic type with three elements: \Single, \Double and \Quad.

❼ Rounding errors: according to the Strict model, a floating-point expression e is represented by three parts: the floating-point part, i.e, the value computed by the program, the exact part which is a real value accessed by the construct \exact(e) and the model part, a real value accessed by the constructs \model(e). The rounding errors are computed by the construct \round_error(e) which is the difference between the floating-point part and the exact part of e, the construct \total_error(e) which is the difference between the floating-point part and the model part of e and the construct \relative_error(e) which is the quotient of \round_error(e) by \exact(e).

A coq formalization of the Strict model

After the specification of the Strict model in the Why syntaxes, the Why tool can now generate verification conditions for programs according to specific annotations. The automatic and interactive provers supported by the Why platform are not until now able to prove these properties. One has to translate our model to them with their syntaxes. We have choose to formalize our models in the Coq proof assistant by using an existing implementation of floating-point arithmetic, the Gappa library [START_REF] Melquiond | Floating-point arithmetic in the Coq system[END_REF]. Why translates automaticaly definitions and predicates from its language to that of Coq.

The gen_float type is defined to be the type float2 of the Gappa library, i.e., the set of couples of integers. The floating-point value is then given by Definition float value x := float2R x.

The round_float operator for rounding reals to floating-point numbers is defined using the functions gappa_rounding and rounding_float as follows: where R is the type of reals in Coq and round_mode is a kind of translate function to the rounding modes of Gappa.

The Full model

The Full model allows overflows, NaNs, infinities and signed zeros. In this model, a floating-point number is still viewed with three parts: the real, exact and model parts. The complete description of this model can be found in the next Why distribution in the file lib/why/floats_full.why.

In the logic: We keep the same types (gen_float, float_format, mode) and functions (float_value, exact_value, mode_value) as in the Strict model. (is finite(x) and is finite(y) and float value(x) <= float value(y)) or (is minus infinity(x) and not is NaN(y)) or (not is NaN(x) and is plus infinity(y)) predicate float lt float(x:gen float,y:gen float) = (is finite(x) and is finite(y) and float value(x) < float value(y)) or (is minus infinity(x) and not is NaN(y)

and not is minus infinity(y)) or (not is NaN(x) and not is plus infinity(x) and is plus infinity(y))

predicate float eq float(x:gen float,y:gen float) = not is NaN(x) and not is NaN(y) and ((is finite(x) and is finite(y) and float value(x) = float value(y)) or (is infinite(x) and is infinite(y) and same sign(x,y))) ❼ a sign_invariant to satisfy the finite_sign axiom as above;

❼ the floating-point value defined by: float_value := float2R genf;

❼ an exact value of type R;

❼ a model value of type R.

Examples

In this subsection, we illustrate the efficiently of our models on several short floating-point C programs. Pragmas allow the user to switch from one model to another. There are three pragmas: the first one is the default model real, the pragma JessieFloatModel(strict) for the Strict model and JessieFloatModel(full) for the Full model. There are also five other pragmas for the five rounding directions of the IEEE standard: the default pragma is #pragma JessieFloatRoundingMode(nearest) and the pragmas:

#pragma JessieFloatRoundingMode(downward) #pragma JessieFloatRoundingMode(upward) #pragma JessieFloatRoundingMode(towardzero) #pragma JessieFloatRoundingMode(towardawayzero)

Polynomial approximation of transcendental functions

The cosine function:

The following annotated C program implements a Taylor polynomial approximation of order 2 of the cosine function in a neighborhood of zero.

#pragma JessieFloatModel(strict)

/*@ requires \abs(x) <= 0x1p-5; @ ensures \abs( \result -\cos(x)) <= 0x1p-23; @*/ float my cosine(float x) { //@ assert \abs(1.0x*x*0.5 -\cos(x)) <= 0x1p-24; return 1.0fx * x * 0.5f; }

Technical appendix

In order to prove the above floating-point C programs, one needs to compile and install the following tools and libraries:

❼ The SVN version of the Frama-C tool [START_REF] Ayad | Framework for modular analyses of C[END_REF].

❼ The CVS version of the Why platform [START_REF] Filliâtre | The Why verification tool[END_REF].

❼ The last version of Coq (Coq 8.2) available in the new official website of the Coq proof assistant: http://logical.saclay.inria.fr/coq/?q= node/60.

❼ The Gappa tool [START_REF] Melquiond | Automatic proof generation of arithmetic properties[END_REF] with its libraries: gappalib-coq-0.9 (for the gappa tactic) and interval-0.8 [START_REF] Melquiond | Interval package for coq[END_REF] (for the interval and interval_intro tactics).

❼ The Coq library on Floating-Point Arithmetic (PFF) which is avalaible on the website: http://lipforge.ens-lyon.fr/www/pff/ ❼ The Coq module JessieGappa.v to translate Jessie goals into Gappa goals.

INRIA 2

 2 Binary floating-point arithmetic 2.1 Binary floating-point representation A binary floating-point number is a real number x represented by x = (-1) s × m × 2 e

INRIA 3

 3 The intermediate languages in the Frama-C chain 3.1 The ACSL language

❼❼❼❼❼❼

  The variables and the function call f e_1 ... e_n. INRIA ❼ The basic unary and binary operations. ❼ The local definition let x = e_1 in e_2 and the sequence of instructions let _ = e_1 in e_2. The local reference let x = ref e_1 in e_2. The dereferenciation !x, the assigment x := e. The conditionnal if e_1 then e_2 else e_3. The loop while e_1 do e_2 done. The function declaration fun (x_1 : t_1) ... (x_n : t_n) -> e. ❼ The recursive function let rec f (x_1 : t_1,...,x_n : t_n) in e.

❼

  A function round_float which is the rounding operator of reals with unbounded range of exponents. ❼ A function max_gen_float which is the maximal real representable number for each float_format. Its values are axiomatized by (constants are written in hexadecimal format): INRIA axiom max single: max gen float(Single) = 0x1.FFFFFEp127 axiom max double: max gen float(Double) = 0x1.FFFFFFFFFFFFFp1023 axiom max quad: max gen float(Quad) = 0x1.FFFFFFFFFFFFFFFFFFFFFFFFFFFFp16383 ❼ A predicate no_overflow which indicates that there is no overflow when the result of an operation is computed with infinite precision and unbounded range of exponents. It is defined by predicate no overflow(f:float format, m:mode, x:real) = abs real(round float(f,m,x)) <= max gen float(f) ❼ Functions to handle rounding errors: gen_round_error for the absolute error, gen_relative_error for the relative error and gen_total_error for the total error. They are defined by function gen round error(x:gen float) : real = abs real(float value(x)exact value(x)) function gen relative error(x:gen float) : real = abs real(float value(x)exact value(x))/exact value(x) function gen total error(x:gen float) : real = abs real(float value(x)model value(x))

  parameter mul gen float : f:float format -> m:mode -> x:gen float -> y:gen float -> { no overflow(f,m,float value(x) * float value(y)) } gen float { float value(result) = round float(f,m,float value(x) * float value(y)) and exact value(result) = exact value(x) * exact value(y) and model value(result) = model value(x) * model value(y) } In some particular functions, like division and square root, this gives conditions on their parameters. For example, div_gen_float is parametrized by: parameter div gen float : f:float format -> m:mode -> x:gen float -> y:gen float -> { float value(y) <> 0.0 and no overflow(f,m,float value(x)/float value(y)) } gen float { float value(result) = round float(f,m,float value(x)/float value(y)) and exact value(result) = exact value(x)/exact value(y) and model value(result) = model value(x)/model value(y) } Comparisons between two gen_float variables x and y is equivalent to the comparison of their associated float values with the mathematical sense. For example, lt_gen_float is specialized by: parameter lt gen float: x:gen float -> y:gen float -> {} bool { if result then float value(x) < float value(y) else float value(x) >= float value(y) } The rounding and the cast operators are specified by: parameter gen float of real : f:float format -> m:mode -> x:real -> { no overflow(f,m,x) } gen float { float value(result) = round float(f,m,x) and exact value(result) = x and model value(result) = x } parameter cast gen float : f:float format -> m:mode -> x:gen float -> { no overflow(f,m,float value(x)) } gen float INRIA { float value(result) = round float(f,m,float value(x)) and exact value(result) = exact value(x) and model value(result) = model value(x) } 4.1.1 Floating-point arithmetic in ACSL with the Strict model

  Definition round float (f : float format) (m : mode) (x:R) := match f with | Single => gappa rounding (rounding float (round mode m) 24 149) x | Double => gappa rounding (rounding float (round mode m) 53 1074) x end.

  In addition, we introduce two new types: the type Float_class with three elements Finite, Infinite and NaN and the type sign with two elments Negative and Positive. Two additional functions float_class and float_sign are necessary to indicate respectively the class and the sign of a gen_float. Special predicates are defined to test if a gen_float variable is finite, NaN, infinities INRIA or signed zero: is_finite, is_NaN, is_plus_infinity, is_minus_infinity, is_gen_zero, is_gen_zero_plus and is_gen_zero_minus. They are naturally defined by: predicate is finite(x:gen float) = float class(x) = Finite predicate is infinite(x:gen float) = float class(x) = Infinite predicate is NaN(x:gen float) = float class(x) = NaN predicate is minus infinity(x:gen float) = is infinite(x) and float sign(x) = Negative predicate is plus infinity(x:gen float) = is infinite(x) and float sign(x) = Positive predicate is gen zero(x:gen float)= is finite(x) and float value(x) = 0.0 predicate is gen zero plus(x:gen float) = is gen zero(x) and float sign(x) = Positive predicate is gen zero minus(x:gen float) = is gen zero(x) and float sign(x) = Negative Comparison between two gen_float numbers is given by the predicates float_le_float, float_lt_float, float_eq_float and float_ne_float: predicate float le float(x:gen float,y:gen float) =

  predicate float ne float(x:gen float,y:gen float) = not float eq float(x,y)where same_sign is defined by: predicate same sign(x:gen float,y:gen float) = float sign(x) = float sign(y){ no overflow(f,m,x) -> is finite(result) and float value(result) = round float(f,m,x) and not no overflow(f,m,x) -> same sign real(result,x) and overflow value(f,m,result) and exact value(result) = x and model value(result) = x } parameter add gen float :f:float format -> m:mode -> x:gen float -> y:gen float -> { } gen float { (is NaN(x) or is NaN(y) -> is NaN(result))and ((is finite(x) and is infinite(y)) -> (is infinite(result) and same sign(result,y))) and ((is infinite(x) and is finite(y)) -> (is infinite(result) and same sign(result,x))) and ((is infinite(x) and is infinite(y) and same sign(x,y)) -> (is infinite(result) and same sign(result,x))) and ((is infinite(x) and is infinite(y) and diff sign(x,y)) -> is NaN(result)) and ((is finite(x) and is finite(y) and no overflow(f,m,float value(x)+float value(y))) -> (is finite(result) and float value(result) = round float(f,m,float value(x)+float value(y)) and sign zero result(m,result))) and ((is finite(x) and is finite(y) and not no overflow(f,m,float value(x)+float value(y))) -> (same sign real(result,float value(x)+float value(y)) and overflow value(f,m,result))) and exact value(result) = exact value(x) + exact value(y) and model value(result) = model value(x) + model value(y) } where diff_sign is defined by: predicate diff sign(x:gen float,y:gen float) = float sign(x) <> float sign(y) parameter mul gen float : RR n➦ 6927 f:float format -> m:mode -> x:gen float -> y:gen float -> { } gen float { ((is NaN(x) or is NaN(y)) -> is NaN(result)) and ((is gen zero(x) and is infinite(y)) -> is NaN(result)) and ((is finite(x) and is infinite(y) and float value(x) <> 0.0) -> is infinite(result)) and ((is infinite(x) and is gen zero(y)) -> is NaN(result)) and ((is infinite(x) and is finite(y) and float value(y) <> 0.0) -> is infinite(result)) and ((is infinite(x) and is infinite(y)) -> is infinite(result)) and ((is finite(x) and is finite(y) and no overflow(f,m,float value(x)*float value(y))) -> (is finite(result) and float value(result) = round float(f,m,float value(x) * float value(y)))) and ((is finite(x) and is finite(y) and not no overflow(f,m,float value(x)*float value(y))) -> (overflow value(f,m,result))) and product sign(result,x,y) and exact value(result) = exact value(x) * exact value(y) and model value(result) = model value(x) * model value(y) } where product_sign is defined by: predicate product sign(z:gen float,x:gen float,y:gen float) = (same sign(x,y) -> float sign(z)= Positive) and (diff sign(x,y) -> float sign(z)= Negative) parameter div gen float : f:float format -> m:mode -> x:gen float -> y:gen float -> { } gen float { ((is NaN(x) or is NaN(y)) -> is NaN(result)) and ((is finite(x) and is infinite(y)) -> is gen zero(result)) and ((is infinite(x) and is finite(y)) -> is infinite(result)) and ((is infinite(x) and is infinite(y)) -> is NaN(result)) INRIA and ((is finite(x) and is finite(y) and float value(y) <> 0.0 and no overflow(f,m,float value(x)/float value(y))) -> (is finite(result) and float value(result) = round float(f,m,float value(x)/float value(y)))) and ((is finite(x) and is finite(y) and float value(y) <> 0.0 and not no overflow(f,m,float value(x)/float value(y))) -> overflow value(f,m,result)) and ((is finite(x) and is gen zero(y) and float value(x) <> 0.0) -> is infinite(result)) and ((is gen zero(x) and is gen zero(y)) -> is NaN(result)) and product sign(result,x,y) and exact value(result) = exact value(x)/exact value(y) and model value(result) = model value(x)/model value(y) } parameter cast gen float : f:float format -> m:mode -> x:gen float -> { } gen float { (is NaN(x) -> is NaN(result)) and (is infinite(x) -> (is infinite(result) and same sign(result,x))) and ((is finite(x) and no overflow(f,m,float value(x))) -> (float class(result) = Finite and float value(result) = round float(f,m,float value(x)))) and ((is finite(x) and not no overflow(f,m,float value(x))) -> same sign(result,x) and overflow value(f,m,result)) and exact value(result) = exact value(x) and model value(result) = model value(x) } parameter lt gen float : x:gen float -> y:gen float -> { } bool { ((is NaN(x) or is NaN(y)) -> result = false) and ((is finite(x) and is infinite(y)) -> if result then float sign(y) = Positive 4.2.2 A coq formalization of the Full model With the Full model, the Coq implementation of the gen_float type is a seven fields record composed of: ❼ the field genf of binary floating-points float2 of the Gappa library; ❼ the Float_class field; ❼ the sign field;

RR n➦ 6927

supported by the French national projects: CerPan (Certification of numerical programs, ANR-05-BLAN-0281-04), Hisseo (Static and dynamic analysis of floatingpoint programs, Digiteo 09/2008-08/2011) and U3CAT (Unification of Critical C Code Analysis Techniques, ANR-09-ARPEGE)

Ali Ayad

In the program: In the full model there are no preconditions on the unary and binary operations. If the result of a floating-point operation is finite and nonzero, the axiom finite_sign ensures that the sign of this result is the same as that of its floating value (using a predicate same_sign_real): axiom finite sign : forall x:gen float.

(is finite(x) and float value(x) <> 0.0) -> same sign real(x,float value(x))

and if the result is zero and the operation is an addition or a subtraction, the predicate sign_zero_resut indicates the sign of the result according to the rounding direction as it is specified in the IEEE-745 standard:

predicate sign zero result(m:mode,x:gen float) = float value(x) = 0.0 -> ((m = down -> float sign(x) = Negative) and (m <> down -> float sign(x) = Positive))

If the result is greater than the maximal floating-point number representable in the destination's format, its class depends on its sign and the current rounding direction as described in Section 2 (paragraph on overflows). This is given by the predicate overflow_value which is defined by: The rounding operator gen_float_of_real, the basic binary and unary operations, the basic operators for comparison and the casting operator cast_gen_float between float formats are specialized according to the Full model. Here we show the complete specification of some of them: ((is finite(x) and is finite(y)) -> if result then float value(x) < float value(y) else float value(x) >= float value(y)) }

Floating-point arithmetic in ACSL with the Full model

In addition to the constructs defined in the Strict model, special predicates are provided on float and double numbers, which check that their arguments are finite numbers, NaNs or infinities respectively: predicate \is NaN(float x); predicate \is NaN(double x); predicate \is finite(float x); predicate \is finite(double x); predicate \is infinite(float x); predicate \is infinite(double x); predicate \is minus infinity(float x); predicate \is minus infinity(double x); predicate \is plus infinity(float x); predicate \is plus infinity(double x);

Useful predicates are also available to express the comparison operators of float (resp. double) numbers as in C: \le_float, \lt_float, \ge_float, \gt_float, \eq_float, \ne_float (resp. for float and double numbers). They are translated to the operators float_le_float, etc ... defined in the Why side.

The sign of a floating-point number is given by the built-in logic function: logic sign \sign(float x); logic sign \sign(double);

where sign is an ACSL built-in logic type with two elements \Positive and \Negative.

Casting a too large real may result into one of the special values ±∞, ±max_float and ±max_double accroding to the current rounding_mode where max_float and max_double are respectively the maximal floating-point numbers represented in the float and double formats.

INRIA

Ali Ayad

The exponential function:

The following annotated C program implements a polynomial approximation of order 2 of the exponential function in a neighborhood of zero by Remez algorithm (the minmax approximation).

#pragma JessieFloatModel(strict)

/*@ requires \abs(x) <= 1.0; @ ensures \abs( \result -\exp(x)) <= 0x1p-4; @*/ double my exp(double x) { /*@ assert \abs(0.9890365552 + 1.130258690*x + 0.5540440796*x*x @ -\exp(x)) <= 0x0.FFFFp-4; @*/ return 0.9890365552 + 1.130258690 * x + 0.5540440796 * x * x; } When compiling these two annotated programs by Frama-C, the proof obligations generated by Why are partially proven by automatic theorem provers thanks to our axiomatization: axiom bounded real no overflow : forall f:float format. forall m:mode. forall x:real. abs real(x) <= max gen float(f) -> no overflow(f,m,x)

Gappa tool proves all these proof obligations except the two assertions which are proven by Coq using the tactic interval [START_REF] Melquiond | Interval package for coq[END_REF] from the Gappa library.

Interval arithmetic

The following annotated C program implements arithmetic floating-point operations on intervals with floating-point bounds. The current rounding mode in this program is Down.

#pragma JessieFloatModel(full) #pragma JessieFloatRoundingMode(downward)

/*@ predicate dif sign(double x, double y) = @ \sign(x) != \sign(y); @ predicate sam sign(double x, double y) = @ \sign(x) == \sign(y); @ predicate double le real(double x,real y) = @ ( \is finite(x) && x <= y) || \is minus infinity(x); @ predicate real le double(real x,double y) = @ ( \is finite(y) && x <= y) || \is plus infinity(y); @ predicate in interval(real a,double l,double u) = @ double le real(l,a) && real le double(a,u); @ predicate is interval ( @ && ( \is infinite(x) ==> y != 0.0 && dif sign(x,y)) @ && ( \is infinite(y) ==> x != 0.0 && dif sign(x,y)); @ assigns \nothing; @ ensures double le real( \result,x*y); @*/ double mul dn(double x, double y) { return x * y; } /*@ requires ! \is NaN(x) && ! \is NaN(y) && sam sign(x,y) && @ ( \is infinite(x) ==> y != 0.0 && \abs(y) >= 0x2.0p-1074 ) @ && @ ( \is infinite(y) ==> x != 0.0 ) @ && @ ( \is finite(y) && ! \no overflow( \Double, \Down,-y) && @ \sign(y) == \Positive ==> x != 0.0 ); @ ensures real le double(x * y, \result); @*/ double mul up(double x, double y) { return -(x * -y); } /*@ requires ! \is NaN(x) && ! \is NaN(y); @ ensures \le float( \result,x) && \le float( \result,y); @ ensures \eq float( \result,x) || \eq float( \result,y); @*/ double min(double x, double y) { return x < y ? x : y;

Ali Ayad } /*@ requires ! \is NaN(x) && ! \is NaN(y); @ ensures \le float(x, \result) && \le float(y, \result); @ ensures \eq float( \result,x) || \eq float( \result,y); @*/ double max(double x, double y) { return x > y ? x : y; } When compiling this annotated program by Frama-C, the proof obligations generated by Why are completely proven by automatic theorem provers thanks to our axiomatization of the two Why models: In this paper, we have presented two formal models to deal with floating-point C programs: the Strict and the Full models. Also, an extention of ACSL and a Coq formalization are given according to these two models. We have seen in Section 4.3 that automatic theorem provers are able to prove the correctness of some floating-point interval arithmetic by using a fruitful axiomatization of the models. Proofs become very long and very complicated when we try to use interactive proof assistants. For example, the proof of the add function of intervals in