Chapter 1

Introduction

The problem of program correctness is always an issue, especially as software applications are more and more complex. Moreover, software correctness plays a crucial role in the overall software quality especially in the automatic control of engines like airplanes or subways. A bug in such systems can not only lost a lot of money but also human life. Thus, the validation of program implementation, which checks whether the code conforms to developer's intentions, is necessary.

A common method for reducing bugs in programs is testing. However, this technique cannot cover all cases of the input space, and thus, although the tests are successful, the implementation may still contain errors. Another possibility is to use formal verification methods. It mathematically proves the correctness of the specification. Although formal verification is used in critical software, there are some aspects which are not studied fully such as computations on floating-point numbers.

Before describing the contributions and the organization of this thesis, we introduce some reasons which cause the discrepancies of results in floating-point programs and give a state of the art of formal verification in general and a state of the art of formal verification for floating-point arithmetic.

Floating-point arithmetic

Nowadays, floating-point computations are widely used in critical systems from domains such as physics, aerospace system, nuclear simulation, etc. In systems which use floating-point computations, the calculations are not exact. This means that there is an approximation for each operation and thus rounding error occurs. There were some famous disasters caused by calculation bugs. The first one is Patriot Missile Failure [START_REF]Patriot Missile Defense: Software Problem Led to System Failure at Dhahran, Saudi Arabia[END_REF] on February 25, 1991. The cause was an inaccurate calculation of the time since boot due to computer arithmetic errors. The second one is the explosion of the Ariane 5 on June 4, 1996, just forty seconds after lift-off. "The internal SRI1 software exception was caused during execution of a data conversion from 64-bit floating point to 16-bit signed integer value. The floating point number which was converted had a value greater than what could be represented by a 16-bit signed integer.", explained the report of the Inquiry Board [START_REF] Lions | Ariane 5 flight 501 failure[END_REF]. The third one is a famous FDIV bug in 1994 in the Intel Pentium processor that costs $500 millions [START_REF] Sharangpani | Statistical analysis of floating point flaw in the Pentium processor[END_REF]. The effort to reduce errors is thus necessary for the safety and reliability of programs.

As we have presented, inaccurate calculation is an issue in critical systems. The reason is that there exist inconsistencies between program executions on different architectures. The IEEE-754 is the IEEE standard for floating-point arithmetic. This standard specifies interchange, arithmetic formats and methods for binary and decimal floating-point arithmetic in computer programming environments. An implementation of a floating-point system conforming to this standard may be realized entirely in software, hardware or in any combination of software and hardware. This standard ensures that we get exactly the same results from the floating-point calculations of a program when it executes in any processor that follows strictly this standard. In Java, the floating-point calculations in a class or a method containing the modifier strictfp2 will follow strictly the standard. In C, the compiler gcc for x86-64 architectures uses -msse and -msse2 by default to enable SSE extensions. For example, a double number is represented in 64 bits; a calculation on double is always rounded in 64 bits. However, the computations do not always follow strictly the IEEE-754 standard because of some architecture-dependent issues. This means that for a double type in C for example, by changing options of the compiler, the calculation is not always done directly in 64 bits. The first architecture-dependent issue is the x87 floating-point unit (FPU) featured in processors of IA32 architecture. It uses the 80-bit internal floating-point registers on the Intel platform. The second one is the fused multiply-add (FMA) instruction, supported by the PowerPC and the Intel Itanium architectures, computes xy ± z with a single rounding. Besides, compiler issues like optimization may also cause discrepancies in results. Because of the issues above, the floating-point computations of a program running on different architectures may be different [START_REF] Monniaux | The pitfalls of verifying floating-point computations[END_REF].

Formal verification

Formal verification is more and more studied and used in both academy and industry research in order to guarantee the safety of programs. There are many approaches for verifying programs formally. Among of them are approaches which verify the source code such as static analysis, others formally verify the machine code implementation.

Static analysis is an approach for checking a program without running it. Deductive verification techniques which perform static analysis of code, rely on the ability of theorem provers to check validity of formulas in first-order logic or even more expressive logics. They usually come with expressive specification languages such as JML [START_REF] Burdy | An overview of JML tools and applications[END_REF][START_REF]JML -Java Modeling Language[END_REF] for Java, ACSL [9] for C, Spec# [6] for C#, etc. to specify the requirements. For automatic analysis of floating-point codes, a successful approach is abstract interpretation based static analysis, that includes Astrée [START_REF] Cousot | The ASTRÉE analyzer[END_REF][START_REF] Monniaux | Analyse statique : de la théorie à la pratique[END_REF] and Fluctuat [START_REF] Delmas | Towards an industrial use of FLUCTUAT on safety-critical avionics software[END_REF].

We introduce here a short history about the formalization of floating-point arithmetic for the specification and the verification. Floating-point arithmetic has been firstly formalized in 1989 by G. Barrett [8]. This is only the formalization of the IEEE standard for binary floating-point arithmetic in the specification language Z. However, there are not any proofs using this specification. Then, in 1995, Carreño and Miner [START_REF] Carreño | Specification of the IEEE-854 floating-point standard in HOL and PVS[END_REF] presented the specification of the IEEE-875 Floating-point standard in HOL and PVS [START_REF] Owre | PVS: A prototype verification system[END_REF]. Two year later, Harrison [START_REF] Harrison | Formal verification of floating point trigonometric functions[END_REF] presented his works on the specification of floating-point in HOL Light -a tool written by himself -about the specification of floating-point arithmetic and demonstrated how to use this specification for proving. The works of Russinoff [START_REF] Russinoff | A mechanically checked proof of IEEE compliance of the floating point multiplication, division and square root algorithms of the AMD-K7 processor[END_REF] are the formal verification of floatingpoint multiplication, division and square root instructions of the AMD-K7 microprocessor in ACL2 [START_REF] Kaufmann | Computer-Aided Reasoning: An Approach[END_REF]. The formalization of IEEE-754 standard in Coq [START_REF]The Coq Proof Assistant[END_REF] was first done in 2001 by Théry, Rideau and Daumas [START_REF] Daumas | A generic library for floating-point numbers and its application to exact computing[END_REF]. The common point of these works is the formalization of the floating-point standard for a specification language.

There are few attempts in specifying and proving behavioral properties of floating-point programs in deductive verification system. In 2004, Miné used relational abstract domains to dectect floating-point run-time errors [START_REF] Miné | Relational abstract domains for the detection of floating-point run-time errors[END_REF]. Leavens presented how to handle NaN problems in floating-point operations for JML in Java in 2006 [START_REF] Leavens | Not a number of floating point problems[END_REF]. Another proposal has been made in 2007 by Boldo and Filliâtre [11]. Authors introduced a methodology and proofs in Coq to perform formal verification of floating-point C programs which was implemented in a tool called Caduceus [START_REF] Filliâtre | Multi-prover verification of C programs[END_REF]. In 2010, Ayad and Marché [START_REF] Ayad | Multi-prover verification of floating-point programs[END_REF] extended this to the support of special values and to the use of automatic theorem provers. This work was implemented in the Jessie plug-in of Frama-C platform. However, these works only follow strictly the IEEE-754 standard, with neither FMA, nor extended registers, nor considering optimization aspects.

Formal verification of machine-code programs is an approach for verifying a program at machine-code level. The specification and correctness of machine-code programs has been first studied in 1961 by Goldstein and von Neumann [START_REF] Goldstine | Planning and coding of problems for an electronic computing instrument[END_REF]. After that, program verification was applied to machine code by Maurer [START_REF] Maurer | Proving the correctness of a flight-director program for an airborne minicomputer[END_REF] in 1976. This paper presented the effort to put techniques for proving the correctness of assembly language and machine language programs into practice. More precisely, it is used for programs written for the Litton C400 airborne computer including overflow analysis, round-off and truncation analysis, fixed-point scaling considerations, etc.. A year later, Floyd-Hoare-style verification condition generator was applied to machine code by Clutterbuck and Carré [START_REF] Clutterbuck | The verification of low-level code[END_REF]. In 1996, Boyer and Yu [START_REF] Boyer | Automated proofs of object code for a widely used microprocessor[END_REF] formalized the MC68020 Instruction Set in the automated reasoning system Nqthm (also known as "the Boyer-Moore theorem prover" and predecessor of ACL2). In this paper, bit vector operations are defined with non-negative integer arithmetic and not for floating-point arithmetic. Another point is that they were able to verify few programs only specified in Nqthm. In 2006, Matthews et al. [START_REF] Matthews | Verification condition generation via theorem proving[END_REF] published a paper on mechanized Hoare logic reasoning for machine code in the form of verification condition generator (VCG) inside the ACL2 theorem prover. Meanwhile, Leroy [START_REF] Leroy | Formal certification of a compiler back-end or: programming a compiler with a proof assistant[END_REF] used Coq for the correctness of an optimizing compiler which takes a significant subset of C as input and produces PowerPC assembly code as output from a simple imperative intermediate language called Cminor. In a separate work on compilation, Li and Slind et al. [START_REF] Li | Structure of a proof-producing compiler for a subset of higher order logic[END_REF] showed that one can compile programs with proof, directly from the logic of the HOL4 theorem prover. In 2007, Pavlova [START_REF] Pavlova | Vérification de bytecode et ses application[END_REF] proposed a framework which allows the verification of Java bytecode programs. In that thesis, she defined both a Bytecode Modeling Language (BML) and a verification condition generator for Java bytecode which is independent from the source code. A year later, Myreen [START_REF] Myreen | Formal verification of machine-code programs[END_REF] made contributions both to approaches for verification of programs and methods for automatically constructing correct code. This work was implemented only for the HOL4 theorem prover.

Although there are a lot of former works on formal verification at low-level programs, none of them considers any aspect of floating-point computation behavioral verification at assembly level.

Contributions

As our goal is to verify floating-point programs in considering architecture and compiler aspect, this thesis consists of two approaches: first, proving floating-point properties by using static analysis for multiple architectures; second, analyzing assembly code and translating it to verification conditions proved by automatic or interactive provers.

The first approach proves numerical programs whatever the environment, or in other words, it proves numerical programs with few restrictions on the compiler and the processor. More precisely, considering rounding-to-nearest mode, double precision numbers and computations, we calculate the rounding error for each operation independent to hardware and the choice of compiler. This approach is implemented in Frama-C platform3 associated with Why [START_REF] Filliâtre | The Why/Krakatoa/Caduceus platform for deductive program verification[END_REF] for static analysis of C code 4 .

The second approach is to prove a program on its assembly code. This means that we prove a program depending on the architecture and compiler. Indeed, with the first approach, we assure that the result is correct in an interval and thus we cannot prove an exact value. Once we compile the program into assembly code, all the necessary information is known such as the precision of each operation, the order of operations, etc. By analyzing assembly code and translating its instructions to verification conditions in the intermediate verification language Why, we can prove a numerical program totally dependent to the architecture and the compiler. A modified version of GAS (GNU Assembler) is used to generate Why verification conditions (VCs) and then these VCs are proved by automatic or interactive provers. Our contributions consist of a translator which allows us to translate all annotations in C program under the forme of inline assmemly and another translator which interpret assembly code into Why. This implementation is integrated in a modified version of GAS that we've mentioned above 5 .

The structure of this dissertation is as follows:

Chapter 2 presents basic knowledge about the IEEE-754 floating-point standard and some architecture-dependent issues. We also present in this chapter a few words about Why intermediate verification language; Frama-C framework and ANSI/ISO C Specification Language (ACSL).

Chapters 3, 4 and 5 present the approach for proving numerical programs for multiple architectures. Chapter 3 studies how to bound a floating-point operation and how to prove a program independent to hardware. Next, Chapter 4 consider the bound when the additions/subtractions re-organize in case of optimization. Chapter 5 then illustrates this approach by examples.

Chapters 6 to 10 will present step by step how to translate assembly code of a numerical program into Why. Chapter 6 introduces an overview of how to prove a C program by using our approach and presents essential elements of assembly language. Chapter 7 presents how to translate assembly instructions to Why with simple programs. Then this translation will be extended for floating-point programs in Chapter 8. Chapter 9 considers programs containing conditional and loop statements. Chapter 10 presents a memory model for programs which contain arrays and pointers. Chapter 11 is a preliminary attempt to support programs with operations at the bit level.

The last chapter gives the summary of contributions and discusses directions of future works.

Chapter 2

Background

This chapter includes the background on the floating-point arithmetic, architecture-dependent and compiler optimization issues. We also presents the Why intermediate language, the ACSL specification language and the Frama-C framework as they will be used in the next chapters.

Floating-point arithmetic

This section consists of two subsections: firstly, we will talk about the IEEE-754 standard; secondly, we will present architecture-independent issues.

The IEEE-754 floating-point standard

The IEEE-754 standard [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Std[END_REF] for floating-point arithmetic was developed to define formats and behaviors for floating-point numbers and computations. There are five basic formats defined in this standard:

• Three binary formats, with encodings in lengths of 32, 64 and 128 bits.

• Two decimal formats, with encodings in lengths of 64 and 128 bits.

In this thesis, we take into account only 32-bit and 64-bit binary formats.

Floating-point data

The set of finite floating-point numbers representable within a particular format is determined by the following integer parameters:

• b = the radix, 2 or 10,

• p = the number of digits in the significand (precision),

• emax = the maximum exponent,

• emin = the minimum exponent where emin = 1 -emax for all formats.

Within each format, the following floating-point data is represented:

• Signed zero or non-zero floating-point numbers of the form • Two infinities: +∞ and -∞,

(-1) s × b e × m ( 2
• Two NaNs: qNaN (quiet) and sNaN (signaling).

There are two kinds of floating-point number: normal and subnormal. The smallest positive normal floating-point number is b min and the largest one is b max × (b -b 1-p ). The non-zero floating-point numbers for a format with magnitude less than b min are called subnormal numbers.

Binary format encodings

In this thesis, we use only binary formats. The representation of floating-point data in the binary format is coded in k bits in the three fields ordered as follow:

• 1-bit sign S,

• w-bit biased exponent E = e + bias, The values of k, p, emax, t, w and bias are listed as follows: The range of the encoding biased exponent E include:

• (t = p-
Parameter binary32 binary64 k,
• every integer in [1, 2 w -2] to encode normal numbers.

• the reserved value 0 to encode ±0 and subnormal numbers.

• the reserved value 2 w -1 to encode ±∞ and NaNs.

The value v of the floating-point datum represented are:

E T v 2 w -1 = 0 NaN 2 w -1 0 (-1) S × (+∞)
[1, 2 w -2] (-1) S × 2 E-bias × (1 + 2 1-p × T ) (normal numbers have an implicit leading significand bit of 1) 0 = 0 (-1) S × 2 emin × (0 + 2 1-p × T ) (subnormal numbers have an implicit leading significand bit of 0) 0 0 (-1) S × (+0) (signed zero)

Rounding error

Let x is a real number, x 1 and x 2 are the two floating-point numbers closest to x such that x 1 ≤ x ≤ x 2 . Let •(x) the rounding value of x to a floating-point number. IEEE-754 standard includes five standard rounding modes:

• Round-to-+∞:

•(x) = x 2 .
• Round-to--∞: •(x) = x 1 .

• Round-to-zero:

•(x) = x 1 if |x 1 | ≤ |x 2 |, x 2 if |x 1 | > |x 2 |.
• Round-to-nearest (ties to even):

•(x) =    x 1 if |x 1 -x| < |x 2 -x|, x 2 if |x 1 -x| > |x 2 -x|,
the one with the bit p -1 of the significand being 0 if |x 1 -x| = |x 2 -x|.

• Round-to-nearest (ties away from zero):

•(x) =        x 1 if |x 1 -x| < |x 2 -x| or (|x 1 -x| = |x 2 -x| and |x 1 | > |x 2 |), x 2 if |x 1 -x| > |x 2 -x| or (|x 1 -x| = |x 2 -x| and |x 1 | ≤ |x 2 |).
When approximating a real number x by its rounding •(x), a rounding error usually happens. We here consider only round-to-nearest mode, that includes both the default rounding mode (ties to even) and the new round-to-nearest, ties away from zero. In radix 2 and round-to-nearest mode, a well-known bound on the error is given by D. Goldberg [START_REF] Goldberg | What every computer scientist should know about floating point arithmetic[END_REF] as follows:

• If a floating-point f = •(x) is such that |f | ≥ 2 e min then we bound the rounding error by the relative error:

x -f x ≤ 2 -p
• For smaller f , the value of the relative error becomes large (up to 0.5). In that case, f is a subnormal number and we prefer a bound based on the absolute error:

|x -f | ≤ 2 e min -p

Architecture-dependent issues

With the same program containing floating-point computations, the result may be different depending on the compiler and the processor. We present in this subsection some architecturedependent issues.

x87 floating-point unit

A well-known cause of discrepancy happens in the IA32 architecture (Intel 386, 486, Pentium etc.) [START_REF] Monniaux | The pitfalls of verifying floating-point computations[END_REF]. The IA32 processors feature a floating-point unit called "x87". This unit has 80-bit registers in "double extended" format (64-bit significand and 15-bit exponent), often associated to the long double C type. When using the x87 mode, the intermediate calculations are computed and stored in the x87 registers (80 bits). The final result is rounded to the destination format. Extended registers may lead to double rounding, where floating-point results are rounded twice. For instance, the operations are computed in the long double type of x87 floating-point registers, then rounded to IEEE double precision type for storage in memory. Double rounding may yield different result from direct rounding to the destination type.

An example is given in Figure 2.1: the long vertical lines are floating-point numbers representable in 64 bits, the short ones are floating-point numbers representable in 80 bits. We assume x is near the midpoint c of two consecutive floating-point numbers a and b in the destination format. Using round-to-nearest, with single rounding, x is rounded to b. However, with double rounding, it may firstly be rounded towards the middle c and then be rounded to a (if a is even). The two obtained results are different.

A C program illustrated for Figure 2.1 is presented in Figure 2.2: the values x = 1.0 and y = 2 -53 + 2 -64 are exactly representable in double precision. With strict IEEE-754 double precision computations for double type, the result obtained is z = 1 + 2 -52 . Otherwise, on IA32, if the computations on double is performed in the long double type inside x87 unit, then converted to double precision, z = 1.0.

Another example which gives inconsistencies between x87 and SSE is presented in Figure 2.3. It is is part of the KB3D [START_REF] Dowek | Conflict detection and resolution for 1[END_REF], an aircraft conflict detection and resolution program. In this example, we have a function int sign(double x) which returns a value which is either -1 if i n t s i g n ( double x ) { i f ( x >= 0 ) r e t u r n 1 ; e l s e r e t u r n -1; } i n t e p s _ l i n e ( double sx , double sy , double vx , double vy ) { r e t u r n s i g n ( s x * vx + s y * vy ) * s i g n ( s x * vys y * vx ) ; } i n t main ( ) { double s x = -0x1 . 0 0 0 0 0 0 0 0 0 0 0 0 1 p0 ; // s x = -1 -2 -52 double vx = -1.0; double s y = 1 . 0 ; double vy = 0 x1 . f f f f f f f f f f f f f p -1; // vy = 1 -2 -53 i n t r = e p s _ l i n e ( sx , sy , vx , vy ) ; p r i n t f ( " R e s u l t = %d\n " , r ) ; } x < 0, or 1 if x ≥ 0. The function int eps_line(double sx, double sy, double vx, double vy) then makes a direction decision depending on a sign after few floating-point computations. We execute this program on SSE unit and obtain that Result = 1. When it is performed on IA32 inside x87 unit, the result is Result = -1.

FMA

Another cause of discrepancy in the result of a floating-point program is the fact that some processors (IBM PowerPC or Intel/HP Itanium) have a fused multiply-add (FMA) instruction which computes x × y ± z as if with unbounded range and precision, and rounds only once to the destination format. This operation can speed up and improve the accuracy of dot product, matrix multiplication and polynomial evaluation. But how should a × b + c × d be computed? When a FMA is available, the compiler may choose either •(a×b+•(c×d)), or •(•(a×b)+c×d), or •(•(a × b) + •(c × d)) which may give different results.

Optimization issues

The last cause for discrepancies occurs when compilers optimize floating-point computations. This includes re-organizing additions or multiplications, use of distributivity, etc. Those mathematically correct identities are usually refuted by floating-point operations.

When we compile a C program with optimization options1 , the assembly code may be changed. Here are some options of gcc:

• -O0: no optimization. This is the default.

• -O, -O1: the compiler tries to reduce code size and execution time, without performing any optimizations.

• -O2: Optimize even more. As compared to -O, this option increases both compilation time and the performance of the generated code. -O2 turns on -finline-small-functions that double f ( double x ) { r e t u r n x /1 e308 ; } double s q u a r e ( double x ) { double y = x * x ; r e t u r n y ; } i n t main ( ) { p r i n t f ("%g \n " , f ( s q u a r e ( 1 E308 ) ) ) ; r e t u r n 0 ; } allows to integrate functions into their callers when their body is smaller than expected function call code.

• -O3: Optimize yet more. It turns on all optimizations specified by -O2 and also other options including -finline-functions.

• -funsafe-math-optimizations: Allow optimizations for floating-point arithmetic, e.g. use associative math, use reciprocal instead of division, disregard floating-point exceptions (division by 0, overflow, underflow, etc).

The options above may change the assembly code and thus affect the result of floating-point calculations except the first and the second one.

Inline function means that the function is integrated into their caller. On x87, inlining may change the result. The program in Figure 2.4 which is presented by D. Monniaux [START_REF] Monniaux | The pitfalls of verifying floating-point computations[END_REF] illustrates this case. When compiling it with -mfpmath=387 -O2, function inlining occurs, the two functions f and square are integrated in the function main and intermediate values are stored in 80-bit registers, not in 64-bit registers. This means that the result is 1E308 and not overflow although it is overflow (the result is Inf) with default option of gcc.

Why-An Intermediate Verification Language

Why2 is a software verification platform [START_REF] Filliâtre | The Why/Krakatoa/Caduceus platform for deductive program verification[END_REF]. It is also known as an intermediate language for program verification. In this dissertation, we use Why as a tool to express our results in the logic world.

Figure 2.5 shows us the process of how Why receives Java and C programs as inputs and generates verification conditions which are proved by automatic or interactive provers. Why has front-ends for C and Java in which the front-end for C is in fact integrated in the Frama-C environment for static analysis of C programs.

One advantage of Why is that it allows us to declare logical models (types, functions, predicates, axioms, lemmas) that can be used in programs and annotations. Another advantage is that the Why tool supports a set of existing provers: • Automatic theorem provers: Alt-Ergo [START_REF] Conchon | CC(X): Efficiently combining equality and solvable theories without canonizers[END_REF], CVC3 [7], Gappa [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], etc.

In the input language of Why, we can define a pure model in the logic world by declaring abstract sort names, declaring logic symbols operating on these sorts and posing first-order axioms to axiomatize the behavior of these symbols. Equality and both integer and real arithmetic are built-in in the logic. We can then declare a set of references which are mutable variables denoting logic values. Finally, we can define procedures which can modify these references. The body of such a procedure is made of statements in a while-style language.

Why includes also logical annotations (pre-and post-conditions, loop invariants and variants, intermediate assertions). Procedures are equipped with pre-and post-conditions. The Why VC generator then produces the necessary VCs to ensures that the body respects the post-condition. One can alternatively just declare procedures by only giving pre-and post-conditions, but also declaring the set of modified references. This feature allows to declare how the atomic operations on a given data type behave.

The detail of the syntax of Why is found at: http://why.lri.fr/. In this section, we will present in a few words axiomatizations, logic functions and parameters which we will use in the thesis.

Logical declarations in Why

Axiomatization

The axiomatization includes abstract types, declarations of functions and predicates, axioms. For example: predicate is_int32(x:int) = -2147483648 <= x and x <= 2147483647 type int32 logic integer_of_int32: int32 -> int axiom int32_coerce: forall x:int32. is_int32(integer_of_int32(x))

Predicate gives a definition based on propositions. For example: predicate is_int32(x:int) = -2147483648 <= x and x <= 2147483647

This predicate says that its argument ranges between -2147483648 and 2147483647.

Type is a declaration of an abstract type. For example:

type int32

The type int32 is an abstract type which represents a 32-bit integer.

Logic function is a function symbol which declares type of inputs and one of output. The body of the function is not defined. Normally, we need to declare axioms to specify logic functions. For example:

logic integer_of_int32: int32 -> int
The function logic integer_of_int32 returns an integer value from an abstract type int32.

Axiom is a logical statement that is assumed to be true.

axiom int32_coerce: forall x:int32, is_int32(integer_of_int32(x))
This axiom says that all value x having type int32 always denotes an integer value in the range of 32-bit word.

Programs in Why

A program in Why can be defined as a function with a body, or declared as a parameter.

Reference declaration with parameter

Parameter can be used to declare a reference (mutable value) in Why. For example: parameter value: int ref

Function with let

A function in Why may contain annotations (preconditions, post-conditions and assertions) and statements.

let f (n:int) = { } r := min !r n { r <= r@ }
The function f takes an integer n as input and assigns the value of min !r n to r. It has no precondition and a post-condition which specifies that the final value of r is not greater than its initial value (denotes by r@).

Function declaration with parameter

Parameter can be used for the declaration of an "abstract program" with precondition and postcondition. For example:

parameter add_int: a:int -> b:int ref -> { } unit writes b { b = b@ + a }
The parameter add_int takes an integer a and a variable b having type int with no precondition, returns unit (nothing). The post-condition is that the new value of b is the sum of a and the previous value of b. Notice that in parameter, there is a declaration of side-effects: the reads and writes keywords declare the set of references possibly accessed and modified, respectively.

Example

In order to know how to write a Why program, a small program is presented as follows:

parameter x: int32 ref parameter inc: n:int32 ref -> { is_int32(integer_of_int32(n) + 1) } unit writes n { integer_of_int32(n) = integer_of_int32(n@) + 1 } let main() = { is_int32(integer_of_int32(x) + 2)} inc x; inc x; void { integer_of_int32(x) = integer_of_int32(x@) + 2 }
All the logical declarations of this program are explained in subsection 2.2.1. We declare a parameter x having type int32, another parameter inc which returns the value of input after adding 1.

In the function main, we assign to x the value of inc !x. This function is used twice. The post-condition of the function main assures that the new value of x is equal to the sum of old There are three obligations generated. The two first ones check the precondition of the parameter inc. As the parameter inc has a precondition, each time it is invoked, the precondition must be proved. In the function main, inc is called twice, this is the reason why we have two obligations that check the precondition of inc. The last one checks the post-condition of the function main.

We present below the VCs for these obligations and show how to prove it manually.

VCs for the first obligation

The VCs for proving the precondition of the first inc x in the function main are:

x: int32 H1: is_int32(integer_of_int32(x) + 2) --------------------------------------------------- is_int32(integer_of_int32(x) + 1)
The hypothesis of this obligation is the precondition of main. As in H1 we have int(x) + 2 which is in range of 32-bit integer and as int(x) is also in the range thanks to axiom int32_coerce, easily we prove that int(x) + 1 is also in this range.

VCs for the second obligation

The VCs for proving the precondition of the second inc x in the function main are:

x: int32 H1: is_int32(integer_of_int32(x) + 2) H2: is_int32(integer_of_int32(x) + 1) x0: int32 H3: integer_of_int32(x0) = integer_of_int32(x) + 1 --------------------------------------------------- is_int32(integer_of_int32(x0) + 1)
They are explained by: H1 : int(x) + 2 is in range of 32 -bit integer.

H2 : int(x) + 1 is in range of 32 -bit integer.

H3 : int(x0) = int(x) + 1

Goal : int(x0) + 1 is in range of 32 -bit integer.

We need to prove that: int(x0) + 1 is in range of 32-bit integer. This means that (int(x) + 1) + 1 or int(x) + 2 (from H2, replace int(x0) with int(x) + 1) must be in range of 32-bit integer. As int(x) is in the range and from H1, this obligation is proved.

VCs for the post-condition obligation

The VCs generated for post-condition of the function main are shown below:

x: int32 H1: is_int32(integer_of_int32(x) + 2) H2: is_int32(integer_of_int32(x) + 1) x0: int32 H3: integer_of_int32(x0) = integer_of_int32(x) + 1 H4: is_int32(integer_of_int32(x0) + 1) x1: int32 H5: integer_of_int32(x1) = integer_of_int32(x0) + 1 --------------------------------------------------- integer_of_int32(x1) = integer_of_int32(x) + 2
This means that H1 : int(x) + 2 is in range of 32 -bit integer.

H2 : int(x) + 1 is in range of 32 -bit integer.

H3 : int(x0) = int(x) + 1 H4 : int(x0) + 1 is in range of 32 -bit integer. H5 : int(x1) = int(x0) + 1 Goal : int(x1) = int(x) + 2
We need to prove that: int(x1) = int(x) + 2. From H5, we replace int(x1) with int(x0) + 1, then from H3, we replace int(x0) with int(x)+1. Now what we have is:

((int(x)+1)+1) = int(x)+2
and the obligation is proved.

Use of Why in this thesis

We all know about the syntax of Why. The question now is how it works in this thesis. In the whole thesis, we use the platform Why to generate verification conditions which are then proved by automatic or interactive provers. / * @ r e q u i r e s \ abs ( x ) <= 0 x1p -5; @ e n s u r e s \ abs ( \ r e s u l t -\ c o s ( x ) ) <= 0 x1p -23; @ * / f l o a t my_cos ( f l o a t x ) { //@ a s s e r t \ abs ( 1 . 0 -x * x * 0 . 5 -\ c o s ( x ) ) <= 0 x1p -24; r e t u r n 1 . 0 fx * x * 0 . 5 f ; }

ACSL

The ANSI/ISO C Specification Language (ACSL)3 is a behavioral specification language for C programs. It is inspired by Caduceus [START_REF] Filliâtre | Multi-prover verification of C programs[END_REF][START_REF] Filliâtre | The Caduceus tool for the verification of C programs[END_REF], itself inspired by JML [START_REF] Burdy | An overview of JML tools and applications[END_REF][START_REF]JML -Java Modeling Language[END_REF]. Its goal is to allow us to formally verify that the implementation of a C function respects its specifications.

The most important ACSL concept is the function contract. A function contract for a C function f is a set of requirements over the arguments of f and/or a set of properties that are ensured at the end of the function.

ACSL annotations are written in special C comments:

• /*@ */ for one or multiple lines,

• //@ for a single-line annotation.

There are three important annotations in ACSL:

Precondition Precondition begins by the keyword requires. It contains the conditions which must hold before calling a function.

Post-condition Post-condition begins with the keyword ensures. It is the conditions which must hold after calling a function.

Assertion Properties which must hold at a program point.

The example in Figure 2.8 contains a C function float my_cos(float x) with ACSL annotations. It provides an approximation of the cosine function on a small interval around 0. The precondition assumes that the absolute value of x is not greater than 2 -5 (=

). The postcondition says that the absolute value of total error (which includes method error and rounding error) is not greater than 2 -23 . The method error is stated in an assertion which says that the bound of method error is 2 -24 . Notice that all the operations in the annotations are done in real, they are not floating-point operations.

Frama-C

Frama-C is a framework for static analysis of C code. This framework is flexible as it is easy to add a new plug-in. There are some ready-to-use plug-ins: Jessie This is deductive verification plug-in. It is based on weakest precondition computation techniques. It aims at deductive verification of behavioral properties of the code, specified using the ACSL language.

In Jessie, there are some ready-to-use models for floating-point arithmetic called strict model and full model. The strict model includes rounding modes; single and double types; basic operations, etc. However, this model forbids exceptional values. The full model is an extension of the strict one. It allows to verify the programs in considering special values such as infinities, NaN, etc.. Figure 2.9 shows us how to prove a numerical program with Frama-C/Jessie. The input of this plug-in is an annotated C program. Its annotations are specified in ACSL. The output of Frama-C/Jessie plug-in is then the input of Why verification condition generator. This generator creates VCs which can be proved by automatic provers (Alt-Ergo, Gappa, CVC3, etc.) or interactive provers (Coq, PVS). The Gappa tool is an automatic prover, which handles formulas made of equalities and inequalities over expressions involving real constants and arithmetic operations. While Gappa does not handle quantifiers, SMT solvers which are a decision procedure that can handle quantifiers, various types of arithmetic and other decidable theories. This is the reason why in our examples, we often combine several automatic provers to prove a program. Now let's back to example in Figure 2.8. All the needed properties about floating-point arithmetic are specified in ACSL. The analysis of this program is done by using Frama-C/Jessie with the strict model which follows strictly the IEEE-754 standard. The verification conditions generated are proved by Gappa except the assertion which is proved using Coq and interval tactic [START_REF] Melquiond | Proving bounds on real-valued functions with computations[END_REF].

Part I

Hardware-independent proofs of numerical programs

Chapter 3

Hardware-independent bounds for floating-point operations

This chapter presents how to bound the rounding error of floating-point operations whatever the architecture and the compiler.

Bounds for one floating-point operation

As we want both correct and interesting properties on a floating-point computation without knowing which rounding will be in fact executed, the chosen approach is to consider only the rounding error. This will be insufficient in some cases, but we believe this can give useful and sufficient results in most cases.

As explained, the choice between 64-bit, 80-bit and double rounding is the main reason that causes the discrepancies of result. We prove a rounding error bound that is valid whatever the hardware, and the chosen rounding (corresponds to interval arithmetic [START_REF] Hickey | Interval arithmetic: From principles to implementation[END_REF][START_REF] Moore | Interval Analysis[END_REF] and standard model [START_REF] Higham | Accuracy and stability of numerical algorithms[END_REF]).

We denote by • 64 the round-to-nearest in the double 64-bit type, by • 80 the round-to-nearest to the extended 80-bit registers.

Theorem 3.1 For a real number x, let (x) be either • 64 (x), or • 80 (x), or the double rounding

• 64 (• 80 (x)). We have either |x| ≥ 2 -1022 and x -(x) x ≤ 2050 × 2 -64 or |x| ≤ 2 -1022 and |x -(x)| ≤ 2049 × 2 -1086 .
The proof is done by case analysis in the following subsections. This theorem is the basis of our approach to correctly prove numerical programs whatever the hardware. These bounds are tight as they are reached in all cases where is the double rounding. They are a little bigger than the ones for 64-bit rounding (2050 and 2049 instead of 2048) for both cases. These bounds are therefore both correct, very tight, and just above the 64-bit's.

As 2 -1022 is a floating-point number, we have (2 -1022 ) = 2 -1022 . As all rounding functions are monotone1 , (x) is also monotone. Then |x| ≥ 2 -1022 implies | (x)| ≥ 2 -1022 and vice versa. Now let us prove the bounds of Theorem 3.1 on the rounding error for all possible values of .

Case 1: Rounding error in 64-bit rounding

From the formulas given by D. Goldberg [START_REF] Goldberg | What every computer scientist should know about floating point arithmetic[END_REF] (explained in subsection 2.1.1), with 64-bit rounding, we have p = 53 and e min = -1022. Therefore:

• With | • 64 (x)| ≥ 2 -1022 : x -• 64 (x) x ≤ 2 -53 • With | • 64 (x)| ≤ 2 -1022 : |x -• 64 (x)| ≤ 2 -1075
Let us see Figure 3.1, the rounding error bounds of 64-bit rounding are above the line and the ones of Theorem 3.1 (T3.1) are below the line. The rounding error bounds of D. Goldberg are smaller than the desired ones, therefore Theorem 3.1 holds in 64-bit rounding.

64

T3.1 0 In Figure 3.2, the rounding error bounds in 80-bit rounding are above the line and the ones of Theorem 3.1 are below the line. Intuitively, the rounding error bounds in case of 80-bit rounding are much smaller in this case than Theorem 3.1's except in the case where |x| is between 2 -16382 and 2 -1022 . There, we have

2 -1022 +∞ x |x -(x)| ≤ 2049 × 2 -1086 |x -• 64 (x)| ≤ 2 -1075 x -• 64 (x) x ≤ 2 -53 x -(x) x ≤ 2050 × 2 -64
• With | • 80 (x)| ≥ 2 -16382 : x -• 80 (x) x ≤ 2 -64 • With | • 80 (x)| ≤ 2 -16382 : |x -• 80 (x)| ≤ 2 -16446 80 T3.1 0 2 -1022 2 -16382 +∞ x |x -(x)| ≤ 2049 × 2 -1086 |x -• 80 (x)| ≤ 2 -16446 x -• 80 (x) x ≤ 2 -64 x -(x) x ≤ 2050 × 2 -64
|x -• 80 (x)| ≤ 2 -64 × |x| ≤ 2 -64 × 2 -1022 = 2 -1086 .
All bounds are smaller than that of Theorem 3.1, so Theorem 3.1 holds in 80-bit rounding.

Case 3: Rounding error in double rounding

The bounds here will be the one of Theorem 3.1. We split in two cases depending on the value of |x|.

Normal range

We first assume that |x| ≥ 2 -1022 . We bound the relative error by some computations and the previous formulas:

x -• 64 (• 80 (x)) x ≤ x -• 80 (x) x + • 80 (x) -• 64 (• 80 (x)) x ≤ 2 -64 + • 80 (x) -• 64 (• 80 (x)) • 80 (x) × • 80 (x) x ≤ 2 -64 + • 80 (x) -• 64 (• 80 (x)) • 80 (x) × • 80 (x) -x x + 1 ≤ 2 -64 +2 -53 × • 80 (x) -x x + 1 ≤ 2 -64 +2 -53 × (2 -64 + 1) = 2050 × 2 -64
(• 80 (x))| ≤ |x -• 80 (x)| + |• 80 (x) -• 64 (• 80 (x))| ≤ 2 -16446 + 2 -1075 < 2 -1086 + 2 -1075 = 2049 × 2 -1086 . If x is in the 80-bit normal range, then 2 -16382 ≤ |x| < 2 -1022 and |x -• 64 (• 80 (x))| ≤ |x -• 80 (x)| + |• 80 (x) -• 64 (• 80 (x))| ≤ 2 -64 × |x| + 2 -1075 ≤ 2 -64 × 2 -1022 + 2 -1075 = 2 -1086 + 2 -1075 = 2049 × 2 -1086 .
Then Theorem 3.1 holds in double rounding.

In conclusion, Theorem 3.1 is proved for all 3 roundings. In practice, instead of dividing into two cases: normal range and subnormal range, we can use the following formula:

|x -(x)| ≤ ε × |x| + η (3.1) with ε = 2050 × 2 -64 and η = 2049 × 2 -1086 .
In strict IEEE-754, where inputs and outputs are on 64 bits, we can set η = 0 for addition and subtraction. Unfortunately here, inputs may be 80-bit numbers so η cannot be set to 0. Note also that absolute value and negation may produce a rounding if we put a 80-bit number into a 64-bit number.

Proof in Coq

We used the Gappa Coq library2 (version 0.13) for Coq (version 8.2) [12] to prove the correctness of Theorem 3.1. The corresponding theorem (see Figure 3.3) and proof (228 lines) are in Coq3 . As we use the Gappa tactic in Coq to prove the theorem, a bound on x is necessary. This is the reason why we add the requirement that |x| ≤ 2 35000 . This value is large enough to satisfy all operations (addition, subtraction, multiplication, division, square root, negation and absolute value) in all types (64 or 80-bit rounding).

Theorem p o s t _ c o n d i t i o

The Coq proof exactly corresponds to the one described. It is not very difficult, but needs many computations and a very large number of subcases. The formal proof gives a very strong guarantee on this result, allowing its use without doubt in the Frama-C platform.

Proofs of numerical programs

In the previous section, we gave a bound for one operation. If we have a program containing a sequence of operations, how can we prove it? In this section, we will use Theorem 3.1 to prove such a program.

When FMA occurs

Theorem 3.1 gives rounding error formulas for various roundings denoted by (64-bit, 80-bit and double rounding). Now, let us consider the FMA that computes x × y + z with one single rounding. The question is whether a FMA was used in a computation. We therefore need an error bound that covers all the possible cases.

The idea is very simple: we consider a FMA as a rounded multiplication followed by a rounded addition. And we only have to consider another possible "rounding" that is the identity:

(x) = x. This specific "rounding" magically solves the FMA problem: the result of a FMA is (x × y + z), that may be considered as 1 ( 2 (x × y) + z) with 2 being the identity. So we handle in the same way all operations even in presence of FMA or not, by considering one rounding for each basic operation (addition, multiplication, etc.). Of course, this "rounding" easily verifies the formulas of Theorem 3.1.

What is the use of this odd rounding? The idea is that each basic operation will be considered as rounded with a that may be one of the four possible roundings. Let us go back to the computation of a*b+c*d: it becomes ( (a × b) + (c × d)) with each being one of the 4 roundings. It gives us 64 possibilities. In fact, only 45 possibilities are allowed (for example, the addition cannot be exact). But all the real possibilities are included in all the considered possibilities. And all considered possibilities have a rounding error bounded by Theorem 3.1.

So, by considering the identity as a rounding like the others, we handle all the possible uses of the FMA in the same way as we handle multiple roundings.

Bounds of a sequence of operations

Theorem 3.2 Let be an operation among addition, subtraction, multiplication, division, square root, negation and absolute value. Let x = (y, z) be the exact result of this operation (without rounding). Then, whatever the architecture and the compiler, the computed result x is such that

If |x| ≥ 2 -1022 , then x ∈ x -2050 × 2 -64 × |x| , x + 2050 × 2 -64 × |x| \ -2 -1022 , 2 -1022 . If |x| ≤ 2 -1022 , then x ∈ x -2049 × 2 -1086 , x + 2049 × 2 -1086 ∩ -2 -1022 , 2 -1022 .
This is straightforward as the formulas of Theorem 3.1 subsume all possible roundings (64 or 80-bit) and operations (using FMA or not), whatever the architecture and the compiler choices.

Theorem 3.3 If we define each result of an operation by the formulas of Theorem 3.2, and if we are able to deduce from these intervals an interval I for the final result, then the really computed final result is in I whatever the architecture and the compiler that preserves the order of operations.

This is proved by using Theorem 3.2.

Implementation

What we have until now is the bound of the rounding error for a sequence of operations. The question is how to use Theorem 3.1 to prove a program in Frama-C/Jessie?

An interesting point of the Jessie plug-in is that we can change the interpretation of floatingpoint operations easily. We define a new "pragma" called multirounding in the same way as the strict model of Jessie plug-in and change the interpretation of each basic floating-point operation (addition, subtraction, multiplication, division, absolute value, square root and negation). In order to do this, we put as post-conditions the formulas of Theorem 3.1 for operations in the Frama-C platform to look into the rounding error of the whole program. Ordinarily, the pragma directive is the method specified by the C standard for providing additional information to the compiler, beyond what is conveyed in the language itself. In our pragma as in Jessie, each floating-point number is represented by two values, an exact one (a real value, as if no rounding occurred) and a rounded one (the true floating-point value). At each computation, we are only able to bound the difference between these two values, without knowing the true rounded value.

For example, the post-condition of the addition is defined by the following predicate: double is an abstract type, double_value and double_exact are logic functions which return a rounded and an exact value from a double. This abstract type and the two logic functions have already defined in the strict and full model and we reuse them in our approach. In this predicate, we also define the exact result which is the addition of two exact values. Notice that the value 0x1.004p-53 = 2050 × 2 -64 and 0x1.002p-1075 = 2049 × 2 -1086 .

The post-conditions for other operations (subtraction, multiplication division, square root, negation, absolute value) are defined in the same way as the addition.

Each operation is declared by a parameter and the post-condition will use the predicate above. The parameter add_double for addition in double, for example, is declared as follows: This parameter returns a double value. The precondition assures that there is no overflow in the addition. The post-condition is the predicate add_double_post.

Notice that the rounding modes that we handle in this approach are rounding-to-nearest modes (ties to even and ties away from zero). Thus, in this Why parameter, either m = nearest_even or m = nearest_away are allowed.

Chapter 4

When the compiler re-organizes a computation

In the previous chapter, we bound the rounding error of each operation provided that the order of the operations is preserved. We also showed that in case of FMA, our method is still valid. Our problem here is that in case of optimizations, the re-organization of operations in a computation may happen. How do we deal with it? We will explain in this chapter how to bound the error of a computation when a reorganization of additions/subtractions occurs.

Associativity for the addition

For the sake of simplicity we will denote floating-point addition by ⊕ and floating-point subtraction by , when the precision is unknown (it may be 64-or 80-bit or double rounding).

Of course, floating-point addition is not associative even if compilers may re-associate additions. For example, if |e| |x|, then (e ⊕ x) x gives zero while e ⊕ (x x) gives e. This catastrophic cancellation is the main problem for the reorganization of additions. The idea here is that we will change the rounding error formula for the addition in order to guarantee that, even if (a + b) + c is transformed into a + (b + c) by the compiler, the bound on the rounding error will still hold. For that, instead of using the formula 64 and η = 2049 × 2 -1086 (apply Formula (3.1) on page 31 for the addition), we will use a formula (with a given ε and η ):

|a ⊕ b -(a + b)| ≤ ε × |a + b| + η with ε = 2050 × 2 -
|a ⊕ b -(a + b)| ≤ ε • (|a| + |b|) + η (4.1)
The new error proportional to |a| + |b| may increase the bound on the rounding error but it handles the cancellations.

To prove a program in multiple environments, we will change the definition of the result of an addition. As we have said in Section 3.3, the post-condition describes how to calculate an operation result. Here, we modify the post-conditions to cover all cases, including the fact that there are several possible results: Note that the strict IEEE-754 definition implies a rounding error formula of the same type but is moreover deterministic. The advantage of modifying the operation post-condition is that it also handles reordering when intermediate values are handled. For example, x=a+b; y=x+c; can be reordered into y=a+(b+c) if x is unused and b+c already computed.

To reason about any ordering of the additions, let us consider a generic algorithm for adding a sequence of numbers [START_REF] Higham | Accuracy and stability of numerical algorithms[END_REF].

Algorithm 1 Let S = {a 0 , . . . , a n }. Repeat while S contains more than one element Remove two numbers x and y from S and add their sum x ⊕ y to S. Return the remaining element.

This generic algorithm is instantiated by the choice of the two numbers that are removed from S. We will denote by o an ordering and by S o n the result of Algorithm 1 for the ordering o. For example, if you choose the preceding computed value and the a i of smaller index, you get the left-associated summation (((a 0 + a 1 ) + a

2 ) + • • • ) + a n .
To ensure the correctness of the approach of taking Formula 4.1 as post-condition, we proved the following theorem for a positive ε. We pose ε n = (1 + ε) n -1.

Theorem 4.1 Let n be an integer such that n ≤ 1 ε , (a i ) 0≤i≤n be a sequence of real numbers and I be a real. Assume that we set the addition post-condition as: x ⊕ y is any real number r such that |r - This means that, if we are able to prove a bound on the rounding error for a sum in a program using our loose post-conditions (Formula (4.1)), then this bound is still correct whatever the compiler reorganization. What is proved using Frama-C and the loose post-conditions (Formula (4.1)) still holds with another ordering (in that case, we use the correct tight post-condition (Formula (3.1) on page 31) proved in the preceding Chapter).

(x + y)| ≤ ε n • (|x| + |y|) + n • η.
Proof.First, we prove an overestimation of |S o 2 n -n 0 a i | with ⊕ having the ε • |x + y| property. We prove by induction on n, that whatever the ordering,

|S o 2 n - n 0 a i | ≤ ε n n 0 |a i | + n 2 η.
If n = 0, then

|S o 2 0 - 0 0 a i | = |a 0 -a 0 | = 0 = ε 0 |a 0 | + 0 2 η
so the property holds.

If n = 1, then

|S o 2 1 - 1 0 a i | = |a 0 ⊕ a 1 -(a 0 + a 1 )| ≤ ε|a 0 + a 1 | + η ≤ ε • (|a 0 | + |a 1 |) + η ≤ ε 1 1 0 |a i | + 1 2 η. If n = 2, then |S o 2 2 - 2 0 a i | = |(a 0 ⊕ a 1 ) ⊕ a 2 -(a 0 + a 1 + a 2 )| ≤ ε|a 0 ⊕ a 1 | + ε|a 2 | + η ≤ ε • (|a 0 | + |a 1 | + |a 2 |) + ε 2 (|a 0 | + |a 1 |) + η + εη ≤ ε 2 • (|a 0 | + |a 1 | + |a 2 |) + 2η ≤ ε 2 2 0 |a i | + 2 2 η.
The other orderings of course give the same property so the overestimation of

|S o 2 n -n 0 a i | holds for n = 2.
Let us assume that the property holds for any value k < n and that n ≥ 2 and let us consider a sequence (a i ) 0≤i≤n+1 and an ordering o 2 . As n + 1 > 1, the value S o 2 n+1 is computed as the sum of two preceding computed values x and y. And x is a computed sum with a known ordering deduced from o 2 of a part of the {a 0 , . . . , a n+1 }. Let I 1 be such that

x ≈ i∈I 1 a i . If k = |I 1 |, then 1 ≤ k ≤ n + 1. Let us denote I 2 = {0, . . . , n + 1} \ I 1 , then |I 2 | = n + 1 -k. |S o 2 n+1 - n+1 0 a i | = |x ⊕ y - n+1 0 a i | ≤ |x ⊕ y -(x + y)| + |x - i∈I 1 a i | + |y - i∈I 2 a i | ≤ ε|x + y| + η + |x - i∈I 1 a i | + |y - i∈I 2 a i | ≤ ε|x - i∈I 1 a i + y - i∈I 2 a i + ( i∈I 1 a i + i∈I 2 a i )| + η + |x - i∈I 1 a i | + |y - i∈I 2 a i | ≤ ε(|x - i∈I 1 a i | + |y - i∈I 2 a i | + | n+1 0 a i |) + η + |x - i∈I 1 a i | + |y - i∈I 2 a i | ≤ ε| n+1 0 a i | + η + (1 + ε)   |x - i∈I 1 a i | + |y - i∈I 2 a i |   ≤ ε n+1 0 |a i | + η + (1 + ε)   |x - i∈I 1 a i | + |y - i∈I 2 a i |  
And x is the sum of (a i ) i∈I 1 with k -1 numbers that is less or equal to n so the induction hypothesis can be used. In a similar way, y is the sum of (a i ) i∈I 2 with n + 1 -k numbers that is less or equal to n. Both are using an ordering that can be deduced from o 2 .

|S o 2 n+1 - n+1 0 a i | ≤ ε n+1 0 |a i | + η + (1 + ε) •   ε k i∈I 1 |a i | + k 2 η + ε n+1-k i∈I 2 |a i | + (n + 1 -k) 2 η   ≤ (ε + (1 + ε) • max(ε k , ε n+1-k ) • n+1 0 |a i | +η • 1 + (1 + ε) • k 2 + (n + 1 -k) 2
As (ε i ) is an increasing sequence, and as k 2 + (n + 1 -k) 2 has its maximum value for k = 1 or k = n, we have:

|S o 2 n+1 - n+1 0 a i | ≤ (ε + (1 + ε)ε n ) • n+1 0 |a i | + η • (1 + (1 + ε)(n 2 + 1)) = ε n+1 • n+1 0 |a i | + η • (1 + (1 + ε)(n 2 + 1))
The last equality is due to this fact:

ε + (1 + ε)ε n = ε + (1 + ε)((1 + ε) n -1) = ε + (1 + ε) n+1 -(1 + ε) = (1 + ε) n+1 -1 = ε n+1 .
Now we bound the η term:

1 + (1 + ε)(n 2 + 1) = n 2 + 1 + (1 + ε) + εn 2
As n ≥ 2, we have

1 + (1 + ε)(n 2 + 1) ≤ n 2 + 1 + n + n • (nε).
And as n ≤ 1 ε , we deduce

1 + (1 + ε)(n 2 + 1) ≤ n 2 + 1 + n + n = (n + 1) 2 .
Therefore,

|S o 2 n+1 - n+1 0 a i | ≤ ε n+1 • n+1 0 |a i | + (n + 1) 2 η,
so this overestimation property holds.

Next, we prove that

ε n n 0 |a i | + n 2 η ≤ I
For that, we use the first hypothesis. The idea is that, if we were able to prove I with the given post-condition, then we may choose each result of an operation (fulfilling this post-condition) and see which error it creates.

For that, we will pose each operation result. More precisely,

• if neither x, nor y is an a i , then we choose for x ⊕ y the value x + y + nη;

• if x = a i and y is not an a i , then we choose for x ⊕ y the value a i + y + ε n |a i | + nη;

• if y = a i and x is not an a i , then we choose for x ⊕ y the value x + a i + ε n |a i | + nη;

• if x = a i and y = a j , then we choose for x ⊕ y the value

a i + a j + ε n |a i | + ε n |a j | + nη.
All those results fulfill the post-condition requirements. Note also that there will be exactly n additions (whatever the ordering). Therefore,

I ≥ |S o 1 n - n 0 a i | = ε n n 0 |a i | + n • nη = ε n n 0 |a i | + n 2 η ≥ |S o 2 n - n 0 a i |
that ends the proof.

We will use the preceding value ε = 2050 × 2 -64 and η = 2049 × 2 -1086 to handle any rounding of one operation. What is proved is that, if we put Formula (4.1) as post-condition of the addition and subtraction with ε = ε n and η = n • η for a given n, then the produced properties will be correct, even if the compiler re-associate the additions. Note that Formula (4.1) subsumes Formula (3.1) for n ≥ 1. Note also that ε n = nε + O(ε 2 ) and that a similar value for bounding the rounding error of a sum can be found by Higham [START_REF] Higham | Accuracy and stability of numerical algorithms[END_REF].

The result is reasonable as the rounding error ε ≈ nε. This is an intuitive demonstration of the optimality where we discard the ε 2 terms (which is reasonable as ε ≈ 2 -53 ) and we discard underflows. We consider we are only allowed to modify the addition post-condition. In that case, if we consider the post-condition of Formula (3.1) and if we study (((a 0 ⊕a 1 )⊕a

2 )⊕a 3 ) • • • , then the final error is at least n • ε • |a 0 | + (n -1) • ε • |a 1 | + • • • + ε • |a n |.
To justify this, we consider an example using 64-bit computations: let a 0 = 1 and a i = 2 -53 , then (((a 0 ⊕ a 1 ) ⊕ a 2 ) ⊕ a 3 ) • • • = a 0 and the error is n • 2 -53 . We are only interested in the first term (n • ε • |a 0 |). Now let us assume we use Formula (4.1) as post-condition and that the program was written a 0 ⊕ (a 1 ⊕ (a 2 ⊕ (a 3 • • • ))) (but the compiler rewrote it in the inverse order). Then the error will be about ε

• |a 0 | + 2 • ε • |a 1 | + • • • + n • ε • |a n |.
As we want this last error to subsume the previous one, we need ε n • ε to make this approach work. We need that the ε of Formula (4.1) be ε n and η will be nη but we do not know n beforehand. The question left is the choice of n. A solution is to look into the program before to have an overestimation of n. We did not put this idea in practice and decided that 16 will be enough. Of course, for linear algebra, it will be insufficient, but for our examples, it will be correct. Moreover, this value can be changed if a bigger value is needed. We will therefore put ε = 2051 • 2 -60 so that ε ≥ ε 16 = 16ε + 256ε 2 (1 + ε) 16 and we will put η = 16η = 2049 × 2 -1082 .

If we constructed a post-condition for a n-term floating-point addition, the results would be better in some cases, but it would not handle the reordering that goes into intermediate values. This is why we chose to only modify the basic block of the numerical program and did so in the best possible way.

|x ⊕ y -(x + y)| ≤ ε • (|x| + |y|) + η |x y -(x -y)| ≤ ε • (|x| + |y|) + η |x ⊗ y -(x * y)| ≤ ε • |x * y| + η |x y -(x/y)| ≤ ε • |x/y| + η • √ x - √ x ≤ ε • √ x + η
and if we are able to deduce a property (such as a rounding error), then this property holds whatever the architecture and the compiler optimizations among commutativity, addition/subtraction associativity (for less than 16 additions/subtractions), use of FMA, use of extended registers, expression factorization and unfactorization.

Without any further work, our method handles other optimizations:

• commutativity: As we have only symmetric formulas for defining the result of an operation, an optimization such as a+b -→ b+a does not endanger our analysis. And commutativity inside associativity is handled by Theorem 4.1.

• expression factorization: As we only consider rounding errors for one operation, the fact that the compiler factorizes or un-factorizes expressions is not a problem. This includes reordering inside intermediate results.

Implementation

Similarly to Chapter 3, we create a new pragma multiroundingR and define the post-conditions by using the formulas of Theorem 4.2 for basic operations in the Frama-C platform to look into the rounding error of the whole program. With this pragma, we let Frama-C know that floating-point computations may be done with extended registers and/or FMA and/or compiler optimizations.

The parameter declared for each operation is the same as the one in the previous chapter. What is changed in this pragma is the predicate in the post-condition of each operation. Here is the predicate for the post-condition of addition. With subtraction, we do the same way as the one of addition.

predicate add_double_post(m:mode,x:double,y:double,result:double)= abs_real(double_value(result)-(double_value(x)+double_value(y))) <=0x1p3000 and abs_real(double_value(result)-(double_value(x)+double_value(y))) <= 0x1.006p-49*(abs_real(double_value(x)) +abs_real(double_value(y)))+0x1.002p-1071 and double_exact(result) = double_exact(x) + double_exact(y)

We use Gappa tool to prove floating-point properties, and Gappa need to bound correctly |x ⊕ y -(x + y)| as a hint. For this reason, we assume that |x ⊕ y -(x + y)| ≤ 2 3000 and this bound is large enough for all cases.

The post-condition of the multiplication is defined as follows:

predicate mul_double_post(m:mode,x:double,y:double,result:double)= abs_real(double_value(result)-double_value(x)*double_value(y)) <=0x1p3000 and abs_real(double_value(result)-double_value(x)*double_value(y)) <=0x1.004p-53*abs_real(double_value(x)*double_value(y))+0x1.002p-1075 and double_exact(result) = double_exact(x) * double_exact(y)

The predicates and parameters of division, negation, square root and absolute value are defined similarly.

Chapter 5

Experimentations

Double rounding example

This first example is the same as the one in Figure 2.2 and is presented in subsection 2.1.2. We add to the original program an assertion written in ACSL (See Figure 5.1). In the modified program, we define two bounds A and B. The assertion we need to prove is that z must be between A and B. Thanks to Gappa, we can find the value of A and B depending on the pragma we use. With the multirounding model, we automatically prove that in every architecture or compiler, the value of z in this program is always in [1, 1 + 2 -52 ]. If we use the multirounding model with addition reorganization, z will be in [1.0 -(2 -49 -2 -53 ), 1.0 + 2 -49 ]. This bound is higher than the previous one but it ensures that with double rounding or not, with reordering of addition or not, the bound is still hold. With the strict model, the value z is proved to be 1 + 2 -52 where only 64-bit rounding is used.

#pragma J e s s i e F l o a t M o d e l ( m u l t i r o u n d i n g ) #d e f i n e E 0 x1 . aap -42 / * @ l o g i c i n t e g e r l _ s i g n ( r e a l x ) = ( x >= 0 . 0 ) ? 1 : -1; * / / * @ r e q u i r e s e1<= x-\exact ( x ) <= e2 ; @ e n s u r e s \ abs ( \ r e s u l t ) <= 1 && @ ( \ r e s u l t != 0 == > \ r e s u l t == l _ s i g n ( \ exact ( x ) ) ) ; @ * / i n t s i g n ( double x , double e1 , double e2 ) { i f ( x > e2 ) r e t u r n 1 ; i f ( x < e1 ) r e t u r n -1; r e t u r n 0 ; } / * @ r e q u i r e s @ s x == \ exact ( s x ) && s y == \ exact ( s y ) && @ vx == \ exact ( vx ) && vy == \ exact ( vy ) && @ \ abs ( s x ) <= 1 0 0 . 0 && \ abs ( s y ) <= 1 0 0 . 0 && @ \ abs ( vx ) <= 1 . 0 && \ abs ( vy ) <= 1 . 0 ; @ e n s u r e s \ r e s u l t != 0 @ == >\ r e s u l t==l

_ s i g n ( \ exact ( s x ) * \ exact ( vx ) + \ exact ( s y ) * \ exact ( vy ) ) @ * l _ s i g n ( \ exact ( s x ) * \ exact ( vy ) -\ exact ( s y ) * \ exact ( vx ) )
; @ * / i n t e p s _ l i n e ( double sx , double sy , double vx , double vy ) { i n t s 1=s i g n ( s x * vx+s y * vy ,-E , E ) ; i n t s 2=s i g n ( s x * vy-s y * vx ,-E , E ) ; r e t u r n s 1 * s 2 ; }

KB3D example

This example includes possible FMA and/or extended registers use but not addition reordering (only one addition). It is part of KB3D [START_REF] Dowek | Provably safe coordinated strategy for distributed conflict resolution[END_REF] 1 , an aircraft conflict detection and resolution program. The aim is to make a decision corresponding to value -1 and 1 to decide if the plane will go to its left or its right. The inputs are the relative position and speed of the other aircraft. Note that KB3D has been formally proved correct using PVS and under the assumption that the calculations are exact [START_REF] Dowek | Provably safe coordinated strategy for distributed conflict resolution[END_REF]. However, in practice, when the computed value is small, the result may be inconsistent or incorrect. The original code is in Figure 2.3 and may give various answers depending on the architecture/compilation. To prove the correctness of this program which is independent to the architecture/compiler, we need to modify this program to detect when the answer may be incorrect.

The modified program (See Figure 5.2) provides an answer that may be 1, -1 or 0. The idea is that, if the result is nonzero, then it is correct. If the result is 0, it means that the result may be under the influence of the rounding errors and the program is unable to give a certified answer. The correctness of the modified program is proved with respect to the following specification: if the result is nonzero, it is the same as if the computations were done on real numbers.

In the original program, the discrepancy of the result is derived from the function int sign(double x)

To use this function only at the specification level, we define a logic function logic integer l_sign (real x)

with the same meaning. Then we define another function int sign (double x, double e1, double e2)

that gives the sign of x provided we know its rounding error is between e 1 and e 2 . In the other cases, the result is zero. The function int eps_line (double sx, double sy, double vx, double vy)

of Figure 5.2 then does the same computations as the one of Figure 2.3, but the result may be different. More precisely, if the modified function gives a nonzero answer, it is the correct one (it gives the correct sign). But it may answer zero (contrary to the original program) when it is unable to give a certified answer. As in interval arithmetic, the program does not lie, but it may not answer. About the other assertions, the given values of sx, vx, sy, vy are reasonable for the position and the speed of the plane. The assertions about s1 and s2 are here to help the automatic provers.

The most interesting part is the value chosen for e 1 and e 2 : they need to bound the rounding error of the computation sx * vx + sy * vy (and its counterpart). Thanks to Gappa, we can find the bound of the rounding error e 1 and e 2 and we can also prove that no overflow occurs.

The following table shows us the value of e 1 and e 2 depending on each model. In practice, e 1 = e 2 = E. Thus, we will show only E in this table:

E Strict model 0x1p-45
Multi-rounding 0x1.90641p-45

Multi-rounding + re-organization 0x1.aap-42

In the usual formalization where all computations directly round to 64 bits, the value E = 0x1p-45 is correct. With multi-rounding approach, we have proved that the value E = 0x1.90641p-45 is correct whatever the architecture. We can also have the value E = 0x1.aap-42 when applying Theorem 4.2 (although it is useless because there is only one addition). This means that the rounding error of sx * vx + sy * vy will always be smaller than this value whatever the architecture or the compiler choices. This means that, even if a FMA is used or if extended registers are used somewhere, this function does not lie.

The analysis of this program (obtained from the verification condition viewer gWhy [START_REF] Filliâtre | The Why/Krakatoa/Caduceus platform for deductive program verification[END_REF]) is given in Figure 5.3. By combining different automatic theorem prover: Alt-Ergo, CVC3, Gappa, we successfully prove all proof obligations in this program.

Summation example

To demonstrate our choices about summation reordering, we use an example by Ogita, Rump and Oishi [START_REF] Ogita | Accurate sum and dot product[END_REF]. Take δ = 2 -54 . Then we add 1, δ, -1, δ 2 and -δ. We denote by ⊕ the 64-bit addition:

• exact computation: 1 + δ + (-1) + δ 2 + (-δ) = δ 2
• left-associated floating-point additions:

(((1 ⊕ δ) ⊕ (-1)) ⊕ δ 2 ) ⊕ (-δ) = -δ • right-associated floating-point additions: 1 ⊕ (δ ⊕ ((-1) ⊕ (δ 2 ⊕ (-δ)))) = 0
From this example, we make a small program (see Figure 5.4). Here, the reordering is critical. With the strict model, we are able to prove that a = -δ, that b = 0 and that the exact values of a and b are equal to δ 2 , and that no overflow occur. But if the compiler reorders these additions (if we had not put parentheses for example), then these proved properties are fallacious. In the multiroundingR pragma, we are only able to prove that the rounding error of a and b is smaller than 0x1.0041p -47 and that no overflow occur. The rounding error is the same for a and b as this rounding error is big enough to cover all possible orderings, including the left-and Of course, the obtained error is bigger than what may really happen as there are cancellations, but this is correct whatever the order of operations. We noticed that Formula (4.1) is especially loose when cancellations happen as the error is proportional to |x| + |y| instead of to |x + y|.

This program is fully proved by Gappa using Why/Frama-C with the multiroundingR pragma.

Clock drift example

Clock drift is a phenomena where two clocks do not run at the same speed. This causes a drift of normal clock compared to the actual time. It seems not very important if our clock drifts only some seconds. However, this becomes dramatic in some cases such as the error in a Patriot missile launcher, and caused several deaths2 . The interesting point of this example is that we try to verify a property on the bound between the computed clock and the exact one.

The example of clock drift [START_REF] Boldo | Formal verification of numerical programs: from C annotated programs to mechanical proofs[END_REF] 3 is presented in Figure 5.5. We assume that the initial value of t is 0.0. At each step, t increases by 0.1. This is a classical example to illustrate rounding errors, that is 0.1 is not a representable floating-point number. Here, the rounding error of 0.1 in double is A = 5.55114e-18. The bound B is the rounding error of t -(t'+(double)0.1) where t' is the previous value of t. At each step i, we assure that |t-i×0.1| is bound by i*(A+B). The post-condition says that after n steps:

|result -n × 0.1| ≤ n × (A + B)
There are several lemmas which are added to help the automatic provers to solve the VCs. This example is proved completely and automatically by the combination of Gappa, Alt-Ergo and CVC3. The value of B changed in function of NMAX, is determined using Gappa, shown in the following table :  #pragma J e s s i e F l o a t M o d e l ( m u l t i r o u n d i n g ) #d e f i n e NMAX 100 #d e f i n e NMAXR 1 0 0 . 0 //@ lemma r e a l _of_int_inf_NMAX : \ f o r a l l i n t e g e r i ; i <= NMAX== > i <= NMAXR; //@ lemma r e a l _of_int_succ : \ f o r a l l i n t e g e r n ; n+1 == n + 1 . 0 ; //@ lemma i n f _ m u l t : \ f o r a l l r e a l x , y , z ; x<=y && 0<=z == > x * z <= y * z ; / * @ r e q u i r e s 0 <= n <= NMAX; @ e n s u r e s \ abs ( \ r e s u l tn * 0 . 1 ) <= n * C ; @ * / double f _ s i n g l e ( i n t n ) { double t = 0 . 0 ; i n t i ;

// A i s a bound o f ( d o u b l e ) 0 . 1 -0 . 1 #d e f i n e A 5 .
/ * @ loop i n v a r i a n t 0 <= i <= n ; @ loop i n v a r i a n t \ abs ( ti * 0 . 1 ) <= i * C ; @ loop v a r i a n t n-i ; @ * / f o r ( i =0; i < n ; i ++) { L :

//@ a s s e r t 0 . 0 <= t <= NMAXR * (0.1+C) ; t = t + 0 . 1 ; //@ a s s e r t \ abs ( t -( \ a t ( t , L ) + ( d o u b l e ) 0 . 1 ) ) <= B ; } r e t u r n t ; }

Scalar product example

The annotated C program in Figure 5.6 computes the scalar product (also known as dot product) of two vectors represented as arrays of doubles:

0≤i<n x i y i
Each vector has n elements. We assume that the values of all the elements of two vectors are not greater than 1.0. Because of the rounding error, the value of p -exact(p) increases after each step of p = p +x i ×y i . This value may change depending on the options of compiler/architecture: it might either follow strictly the standard IEEE-754, or use x87 with 80-bit internal registers, or x87 with optimization, or FMA. In this example, we use a bound B for the rounding error of p-(p +x i ×y i ) where p is the previous value of p. The post-condition expresses a bound B on the accumulated rounding error in function of a bound NMAX on the size of the vectors. Several extra assertions are added in the body of the loop: these are needed to help the automatic provers to solve the generated VCs. In particular, to make Gappa solve the VCs on the accumulated rounding error, it is necessary to guarantee that p remains bounded: it is bounded by NMAX(1+B).

The value of B in function of NMAX when proving the program with the strict model and our multirounding and multirounding with the addition re-organization are shown in the table below: / * @ lemma bound_int_to_ r e a l : @ \ f o r a l l i n t e g e r i ; i <= NMAX == > i <= NMAXR; * / / * @ lemma t r i a n g u l a r _ i n e q u a l i t y : @ \ f o r a l l r e a l x , y , z ; \ abs ( x-z ) <= \ abs ( x-y ) + \ abs ( y-z ) ; * / / * @ r e q u i r e s 0 <= n <= NMAX; @ r e q u i r e s \ v a l i d _ r a n g e ( x , 0 , n-1) && \ v a l i d _ r a n g e ( y , 0 , n-1) ; @ r e q u i r e s \ f o r a l l i n t e g e r i ;

0 <= i < n == > @ \ abs ( x [ i ] ) <= 1 . 0 && \ abs ( y [ i ] ) <= 1 . 0 && @ x [ i ] == \ exact ( x [ i ] ) && y [ i ] ==\exact ( y [ i ] ) ; @ e n s u r e s \ abs ( \ r e s u l t -\ exact ( \ r e s u l t ) ) <= n * B ; * / double s c a l a r _ p r o d u c t ( double x [ ] , double y [ ] , i n t n ) { double p = 0 . 0 ; / * @ loop i n v a r i a n t 0 <= i <= n ; @ loop i n v a r i a n t \ abs ( \ exact ( p ) ) <= i ; @ loop i n v a r i a n t \ abs ( p -\ exact ( p ) ) <= i * B ; @ loop v a r i a n t n-i ; * / f o r ( i n t i =0; i < n ; i ++) {
// bounds , n e e d e d by Gappa //@ a s s e r t \ abs ( x [ i ] ) <= 1 . 0 ; //@ a s s e r t \ abs ( y [ i ] ) <= 1 . 0 ; //@ a s s e r t \ abs ( p ) <= NMAXR * (1+B) ;

L : p = p + x [ i ] * y [ i ] ;
// bound on t h e r o u n d i n g e r r o r s i n t h e s t a t e m e n t above , p r o v e d by gappa / * @ a s s e r t \ abs ( p -( \ a t ( p , L ) + x [ i ] * y [ i ] ) ) <= B ; * / // t h e p r o p e r i n s t a n c e o f t r i a n g u l a r i n e q u a l i t y t o show t h e main i n v a r i a n t / * @ a s s e r t \ abs ( p

-\ exact ( p ) ) <= \ abs ( p -( \ a t ( p , L ) + x [ i ] * y [ i ] ) ) + \ abs ( ( \ a t ( p , L ) + x [ i ] * y [ i ] ) -( \ exact ( \ a t ( p , L ) ) + x [ i ] * y [ i ] ) ) ; * / // a lemma t o show t h e i n v a r i a n t \ abs ( \ exact ( p ) ) <= i / * @ a s s e r t \ abs ( \ exact ( p ) ) <= \ abs ( \ exact ( \ a t ( p , L ) ) ) + \ abs ( x [ i ] ) * \ abs ( y [ i ] ) ; * / //

Hardware-dependent proofs of numerical programs

Chapter 6

Principle of proofs on assembly code with Why

This chapter begins a new part of the thesis. Remind that the goal of this part is to show us how to prove a C program on its assembly code. In order to do it, we present firstly the principle of the proofs on assembly code with Why, then we present in each chapter how to translate an assembly program into Why from the simplest cases to the complicated ones. Before talking about the translation of assembly code into Why, this chapter will give a general view of our approach and some basic knowledge about assembly code.

Steps of proofs

Our approach for proving a C source via analyzing its assembly follows several steps illustrated on In a C program with annotations written by ACSL, all annotations are put in comments. When gcc generates assembly code, these annotations will be ignored. As annotations are important to prove the program, a preparation step is needed. This step puts all annotations into inline assembly in order to keep them in assembly code. This will be detailed in Section 6.3.

Once the preparation step is done, another C file containing inline assembly is generated. The regular GNU compiler gcc is called with option -S to generate assembly code from this C file.

The translation from assembly code to Why is implemented in our own modified version of the GNU assembler GAS. This step generates a file containing proof obligations in Why. These obligations are then attempted to be proved by automatic or interactive provers.

Essential elements of assembly language

An assembly language is a low-level programming language. It is directly influenced by the instruction set and architecture of the processor. A program written in assembly language consists of a series of statements. They are translated by an assembler to a stream of executable instructions that can be loaded into memory and executed.

There exist many different assemblers which translate assembly source code into binary programs such as NASM1 , MASM2 , GAS3 , etc. The GNU Assembler, commonly known as GAS, is the default back-end of gcc and it is a part of the GNU Binutils package. By default, on the x86 and x86-64 architecture, it uses the AT&T assembler syntax. We use gcc to generate assembly code. This chapter will give a background on assembly language using the AT&T syntax.

Assembly statements are entered one per line in the source file. They use the same format:

[label] mnemonic [operands] [comment]
The fields in the square brackets are optional.

Operands and Instruction Naming

Operands

An operand in assembly language may be a register, a memory reference or a constant. Note that with AT&T syntax, the first operands are the sources and the last one is the destination (if exists).

Registers are preceded by '%'. For example: the EAX register is specified as %eax Memory references Memory references in AT&T syntax has the following form: section:disp(base, index, scale) where base and index are the optional 32-bit base and index registers, disp is the optional displacement, and scale, taking the values 1, 2, 4, 8 and multiplies index to calculate the address of the operand. If there is no scale specified, it takes 1. section specifies the optional section register for memory operand.

The address of a memory reference is calculated by

addr(m) = base + scale × index + disp
For example:

• -4(%rbp): base is '%rbp'; disp is '-4'. index, scale are both missing.

• foo(,%eax,4): index is '%eax'; scale is '4'; disp is 'foo'. All others fields are missing.

In our model, we suppose that there is no section, this means that there is only one section in the memory.

The x86-64 architectures add an RIP (instruction pointer relative) addressing. This addressing mode is specified by using '%rip' as a base register. For example:

• 1234(%rip) points to the address 1234 bytes past the end of the current instruction.

• .LC0(%rip) points to the symbol .LC0 in RIP relative way [START_REF] Elsner | Using as[END_REF].

Immediate operands are preceded by '$' and written in decimal or in hexadecimal.

For example: $1, $0x3f800000

Instruction Naming

In AT&T syntax, instruction mnemonics are suffixed with one character modifiers which specify the size of operands. The letter 'b', 'w', 'l', and 'q' specify byte, word, long and quadruple words operands. If no suffix is specified, GAS will try to fill in the missing suffix based on the destination register operand. For example, in the instruction movl -4(%rbp), %eax movl is an instruction mnemonic with the suffix 'l'. This means that the instruction copies data from 32-bit source (-4%rbp)) to 32-bit destination (%eax).

movq %rdi, %rsi movq is an instruction mnemonic with the suffix 'q'. This instruction moves data from 64-bit source (%rdi) to 64-bit destination (%rsi).

EFLAGS register

The 32-bit EFLAGS register holds the state of the processor. It is modified by many instructions and is used for comparing some parameters, conditional loops and conditional jumps. This register contains a group of status flags, a control flag, and a group of system flags. The status flags are presented as follows:

Carry Flag(CF) Set if an arithmetic operation generates a carry or a borrow out of the mostsignificant bit of the result; cleared otherwise.

Parity Flag(PF) Set if the least-significant byte of the result contains an even number of 1 bits; cleared otherwise.

Adjust Flag(AF) Set if an arithmetic operation generates a carry or a borrow out of bit 3 of the result; cleared otherwise. This flag is used in binary-coded decimal (BCD) arithmetic

Zero Flag(ZF) Set if the result is zero; cleared otherwise.

Sign Flag(SF) 0 indicates a positive value and 1 indicates a negative value.

Overflow Flag(OF) Set if the integer result is too large a positive number or too small a negative number (excluding the sign-bit) to fit in the destination operand; cleared otherwise.

General-purpose instructions

In this section, we only talk about the general-purpose instructions, floating-point instructions will be presented in Chapter 8. General-purpose instructions are divided into several groups as follows:

Data transfer instructions

The mov instruction transfers data from source operand to destination operand. It requires two operands and has the syntax:

mov src, dest
The data is copied from src to dest and the src operand remains unchanged. Both operands should be of the same size. The mov instruction can take one of the following five forms:

• mov register, register

• mov immediate, register

• mov immediate, memory

• mov register, memory

• mov memory, register

Binary Arithmetic Instructions

The following instructions make a binary calculation. They can be used to add/sub/mul/div two 8-, 16-, 32-or 64-bit operands.

add src, dest dest = dest + src sub src, dest dest = dest -src mul src, dest dest = dest * src div src, dest dest = dest / src Notice that with multiplication and division operations, there exist the signed one and the unsigned one. The unsigned instructions are mul and div. The signed ones use the prefix 'i': imul, idiv.

Calling procedures using call and ret

Stack frame

For function handling in assembly code, the following elements are needed:

• Stack has one stack frame per active function invocation

• Stack pointer register (ESP) points to top (low memory) of current stack frame

• Base pointer register (EBP) points to bottom (high memory) of current stack frame Each stack frame contains

• Return address (Old EIP4 (instruction pointer register, points to next instruction to be executed))

• Old EBP

• Saved register values

• Local variables

• Parameters to be passed to callee function call and ret instructions

The call instruction allows control transfers to procedures within the current code segment (near call) and in a different code segment (far call). Near calls usually provide access to local procedures within the currently running program or task. Far calls are usually used to access operating system procedures or procedures in a different task.

The ret instruction also allows near and far returns to match the near and far versions of the call instruction.

In this thesis, we only talk about near call and ret operation.

When executing a near call, the processor does the following steps:

1. Pushes the current value of the EIP register on the stack.

2. Loads the offset of the called procedure in the EIP register and begins execution of the called procedure.

When executing a near return, the processor performs these actions:

1. Pops the top-of-stack value (the return instruction pointer) into the EIP register.

2. Resumes execution of the calling procedure.

Some assembler directives

All directives are specified by D. Elsner et al. [START_REF] Elsner | Using as[END_REF]. In this section, we present some of them which often appear in our examples.

Comm directive

The syntax of Comm directive is:

.comm symbol, length

It declares a common symbol named symbol. When linking, a common symbol in one object file may be merged with a defined or common symbol of the same name in another object file. If the GNU linker ld does not see a definition for the symbol -just one or more common symbols -then it will allocate length bytes of uninitialized memory. The argument length must be an absolute expression. If ld sees multiple common symbols with the same name, and they do not all have the same size, it will allocate space using the largest size.

Global directive

Its syntax is:

.globl symbol .global symbol

The directive .global makes the symbol visible to ld. If we define symbol in our partial program, its value is made available to other partial programs that are linked with it. Otherwise, symbol takes its attributes from a symbol of the same name from another file linked into the same program. Both spellings .globl and .global are accepted, for compatibility with other assemblers.

CFI directives

The directive .cfi_startproc is used at the beginning of each function that should have an entry in .eh_frame. It initializes some internal data structures. The directive .cfi_endproc is used at the end of a function where it closes its unwind entry previously opened by .cfi_startproc, and emits it to .eh_frame.

Inline assembly

Inline assembly is a way to put directly assembly in high-level source code. We present here the basic syntax of inline assembly which will be used in the next chapters.

Simple Inline Statement

The form of a basic inline statement is: asm("assembly code"); For example: asm("move %eax, %ebx");

Extended Inline Statements

In basic inline assembly, we have only instructions. In extended assembly, we can also specify the operands. The format of the asm statement consists of four components below: asm( assembly template :outputs /* optional */ :inputs /* optional */ :clobber list /* optional */ ); where each component is separated by a colon (:). The last three components are optional.

Assembly template consists of the assembly language statements to be inserted into the C code. This may be a single instruction or a sequence of instructions.

Outputs specify the output operands for the assembly code. The format specifying each operand is "=option-constraint" where option-constraint may be:

r : register operand constraint m : memory operand constraint rm : register or memory ri : register or immediate g : general X : any operand whatsoever is allowed 5 .

Inputs are specified in the same way, except for the '=' sign.

Clobber list is the list of registers modified by the assembly instructions.

The operands specified in the output and input parts are assigned sequence numbers 0, 1, 2, etc. For example: asm( "movl %0, %1" :"=r"(sum)/* output */ :"r"(number)/* input */ );

The C variables sum and number are both mapped to registers. In assembly statement, sum is identified by %0 and number by %1.

We can put the keyword volatile after asm if our assembly statement must execute where it is put. Its form is: asm volatile (...: ...: ...: ...)

Preparation of source code

As said in Section 6.1, the goal of this step is to keep annotations when compiling to assembly code. The idea of this step is that we create a new C file in which all annotations are put in inline assembly statements. The assembly code will be generated from this new C file, not the original one.

By using inline assembly in C program, all the variables in annotations will be replaced by a memory reference or a register in assembly code.

For example, an annotation in ACSL :

/*@ requires n >= 0 && n < 100;*/ is put in inline assembly under the following format:

asm volatile("/*requires #int#%0# >= 0 && #int#%1# < 100;*/"::"X"(n),"X"(n));

where the first occurrence and the second one of n in this annotation are represented by "%0", "%1". Notice that in assembly code, we do not know the type of variables. This is the reason why in this step, we add type of variables into annotations. In the example above, each variable is replaced by #type_of_variable#index#. We indicate that this variable is input option of inline assembly statement with the constraint "X". The inline assembly above will then be translated in assembly code as follows:

/*requires #int#-20(%rbp)# >= 0 && #int#-20(%rbp)# < 100;*/
In the programs containing labels like in the following piece of code: for(i=0; i < n; i++) { L: t = t + 0.1f; //@ assert \abs(t-(t@L+\round_single(\nearest_even,0.

1)))<=B; }

There is a label L and in the assertion, we use t@L which means that we get the value of t at label L. When compiling this code into assembly, the label L will disappear. In order to keep this label in assembly code, we also put it in inline assembly statement. In this case, the inline assembly is asm volatile("L:");

An another point we process in this step is to have only one returned value. This is not a new work, it have been described in CIL (C Intermediate Language) 6 . This ensures that each function has one return statement. For example, the function sign below has three returns statements:

int sign(double x, double e1, double e2) { if (x > e2) return 1; if (x < e1) return -1;
return 0; } The output after the preparation step is int sign(double x, double e1, double e2) { int res;

if (x > e2) {res = 1; goto resL;} if (x < e1) {res = -1; goto resL;} {res = 0; goto resL;} resL: return res; }

In the output of this step, we define a local variable res (if this variable does not exist in the program) with its type is the returned type of the function. In this case, the returned type is int. Then, we replace each return statement return A; by {res = A; goto resL;}. Finally, at the end of the function, we add resL: return res;.

In brief, what we do in this step are

• We replace all the macros in annotations

• We translate annotations under the form of inline assembly statements. We also keep the type of variables in annotations. This is very important because at assembly level, we do not know the type of the registers/memory references.

• We put the labels in C program into inline assembly so that it will appear in assembly code with the same name.

• We have only one return statement in each function.

• Precondition and post-condition are put at the right place. This is important when the function is integrated into the caller. The precondition and post-condition are also integrated with inline function in the right position in assembly code.

/ * @ r e q u i r e s n >= 0 && n < 1 0 0 ; * / i n t f ( i n t n ) { i n t tmp = 100 -n ; //@ a s s e r t tmp > 0 ; //@ a s s e r t tmp <= 1 0 0 ;

r e t u r n tmp ; } Figure 6.2: A simple program i n t f ( i n t n ) {
asm v o l a t i l e ( " / * r e q u i r e s #i n t#%0#>=0 && #i n t #%1#<100; * / " : : " X" ( n ) , "X" ( n ) ) ; i n t r e s ; i n t tmp = 100 -n ; asm v o l a t i l e ( " / * a s s e r t #i n t#%0# > 0 ; * / " : : " X" ( tmp ) ) ; asm v o l a t i l e ( " / * a s s e r t #i n t#%0# <= 1 0 0 ; * / " : : " X" ( tmp ) ) ; { r e s = tmp ; goto r e s L ; } r e s L : r e t u r n r e s ; } The advantage of using inline assembly statements is that when the program is compiled with gcc -S, the compiler will replace all the variables in the annotations automatically by registers/memory refenrences, even with optimization options. However, the assembly code generated may be modified if we use the complex statement of inline assembly. Adding more move instructions is one of the modifications.

Examples

Now let us see a first simple example (Figure 6.2) to know about a program in assembly code and what we obtain in assembly code after passing the preparation step.

The example in Figure 6.2 has a function int f(int n) that returns the value of 100 -n. The precondition of this function is 0 ≤ n < 100. We have two assertions in the body of the function. These are tmp > 0 and tmp ≤ 100.

After the preparation step, a new C file containing the program as in Figure 6.3 is created. Precondition and two assertions are put into inline assembly statement. A new variable int res is added for returning.

The assembly code of this program (See Figure 6.4) is generated from program in Figure 6.3 by the default option of gcc. There are only three basic instructions used: transfer data with mov instruction, sub instruction and instructions for returning a function.

The function f in assembly code is defined as a global symbol with type @function(line 1-2). This means that this function is visible in other files. A label f begins this function. The body is between two directives .cfi_startproc and .cfi_endproc.

As one can see, ACSL annotations appear between #APP and #NO_APP in assembly code. / * @ r e q u i r e s 0 <= x <= 1000 ; @ e n s u r e s \ r e s u l t == x * x ; * / i n t s q u a r e ( i n t x ) { i n t tmp = x * x ; r e t u r n tmp ; } i n t main ( ) { i n t a = 5 ; i n t b = s q u a r e ( a ) ; //@ a s s e r t b== 2 5 ; r e t u r n 0 ; } The second example we illustrate here is a program containing a function which is called by another (See Figure 6.5). In this example, there are two functions:

• The function int square(int x): calculates the square value of an integer; i n t s q u a r e ( i n t x ) { asm v o l a t i l e ( " / * r e q u i r e s 0 <= #i n t#%0# <= 1000 ; * / " : : " X" ( x ) ) ; i n t r e s ; i n t tmp = x * x ; { r e s = tmp ; goto r e s L ; } r e s L : asm v o l a t i l e ( " / * e n s u r e s #i n t#%0# == #i n t#%1# * #i n t #%2#; * / " : : " X" ( r e s ) , "X" ( x ) , "X" ( x ) ) ; r e t u r n r e s ; } i n t main ( ) { i n t r e s ;

i n t a = 5 ; i n t b = s q u a r e ( a ) ; asm v o l a t i l e ( " / * a s s e r t #i n t#%0#== 2 5 ; * / " : : " X" ( b ) ) ; { r e s = 0 ; goto r e s L ; } r e s L : r e t u r n r e s ; } The new C program generated after the preparation step is illustrated in Figure 6.6. Like previous example, all annotations are put in inline assembly statement. One point we want to emphasize in this example is that the keyword \result is replaced by the variable res in post-condition. We insist that this translation is done by our translator called "preparation of source code". The assembly code of this program is in Figure 6.7. The function square is called in assembly code by the instruction call square (line 41). Thank to the new variable res in the new program in Figure 6.6, the returned value is replaced by the reference memory -4(%rbp) without any further step (line 24). If we do not replace \result in annotation, it is difficult to determine which register/memory reference is used for returning value.

. f i l e " s q u a r e _ i n l i n e . c "

. t e x t . g l o b l s q u a r e . t y p e s q u a r e , @ f u n c t i o n s q u a r e : .LFB0 :

. c f i _ s t a r t p r o c

Case of Simple programs

In the previous chapter, we have presented basic knowledge about assembly code and a general view about how to prove a C program from their assembly code with our approach. We also talked about the preparation step for obtaining the registers/memory references corresponding to the variables in annotations.

In this chapter, we describe how to translate a very simple program C into Why including the translation of instructions and annotations.

Definition of the class of "simple" C programs

Simple C programs considered in this chapter are made of a set of function definitions, specified with ACSL-style annotations, which satisfy these restrictions:

• The only data types int and long int are assumed to denote 32-bit and 64-bit 2complement integers. In particular there are no float types, no arrays and no pointers.

• There are no global variables but only local variables and arguments of the functions.

• The body of any function is restricted to a sequence of assignments, i.e. there is no compound instructions: no loop statements of any kind, neither if, nor switch and nor goto statements.

• The allowed expressions are the arithmetic expressions and the function calls.

This class of programs is simple for us because the corresponding assembly codes contain only general-purpose instructions, neither jump instructions nor any floating-point instructions.

Translation to Why

Now we will detail the translation of assembly code to Why. We present firstly how to translate operands to Why. Secondly, we will talk about how to translate annotations to Why. Finally, we present the translation of instructions.

Translation of 32-bit and 64-bit integers

Why has only unbounded mathematical integers built-in. Thus, 32-bit integers must be defined in Why. We follow here the same technique as what is done in the Jessie plug-in of Frama-C. The type int32 is an abstract type for an 32-bit integer.

type int32 logic integer_of_int32: int32 -> int

The logic function integer_of_int32 returns an integer value from an int32. We need a predicate is_int32 which verifies whether an integer is in the range of 32-bit word or not. predicate is_int32(x: int) = -2147483648 <= x and x <= 2147483647

The axiom axiom int32_coerce: forall x:int32, is_int32(integer_of_int32(x)) assures that all value x having type int32 always denotes an integer value in the range of 32-bit word.

A 64-bit integer has type int64. Like int32, we define a following logic function and a predicate for it: type int64 logic integer_of_int64: int64 -> int predicate is_int64(x: int) = -9223372036854775808 <= x and x <= 9223372036854775807

We also have an axiom

axiom int64_coerce: forall x:int64, is_int64(integer_of_int64(x))
Although the 8-, 16-integers are not considered here for simplicity, they are handled similarly, as in Jessie. We also consider signed integers and not unsigned integers, they are handled similarly as well.

Translation of operands

We want to translate operands being registers or memory references into Why variables. To do so, we make the following hypothesis: Assumption 7.1 (Separation Assumption) On a simple C program, the compiler generates an assembly code where syntactically distinct memory references denote disjoint memory locations.

For example, we assume that in any assembly code, the memory references -16(%rbp) and -8(%rax) are disjoint. Of course there is no reason that this is true in general, but we claim that for the "simple" C programs considered here, and our GNU compiler, this is true. Note that in Chapter 10, this assumption will not be made anymore.

The Separation Assumption allows us to translate each memory reference into a Why variable whose name is syntactically derived from it.

In assembly programs, %ax, %eax, %rax indicate the same register with different size:

0 16 32 64 %ax %eax %rax
As illustrated in the figure above, the size of %ax, %eax and %rax are 16, 32, 64 bits, respectively. When translated into Why, they have the same name: _rax. Once we have this, we do not need to cast from %ax and %eax to %rax and otherwise. With other general-purpose registers such as %bx, %cx, %di, %si, etc. we do similarly.

The following abstract type will be used in this chapter: type register Each register or memory reference used as an operand will be declared in Why as a variable with type register. In order to get a 32-bit and a 64-bit integer from a register, we need the two following logic functions:

logic sel_int32: register -> int32 logic sel_int64: register -> int64
The logic function sel_int32 and sel_int64 returns a 32-bit and a 64-bit integer from a register.

We distinguish two types of operands as follows:

• Immediate operand: begins by '$'. There is not any declaration here because this operand is a constant. We just delete the prefix '$'.

• Registers and memory references: We denote by op the variable in Why corresponding to the operand. Each register or memory reference will have a unique name in Why. In this document, we name the register or memory reference by replacing all special character (()+-%.) by '_'.

Example: -4(%rbp) = _4__rbp_. The Why variable to declare for this operand is parameter _4__rbp_: register ref

We denote by op int32 and by op int64 Why expressions (of Why type int) which denote the integer that op represents as a signed 32-bit or 64-bit integer. They are defined by:

$imm int32 = imm reg int32 = (integer_of_int32 (sel_int32 !reg)) mem int32 = (integer_of_int32 (sel_int32 !mem)) $imm int64 = imm reg int64 = (integer_of_int64 (sel_int64 !reg)) mem int64 = (integer_of_int64 (sel_int64 !mem))
For example, with the following instruction:

movl $5, -4(%rbp)
we interpret its operands into Why by

$5 int32 = 5 -4(%rbp) int32 = (integer_of_int32 (sel_int32 !_4__rbp_))
The instruction movq %rax, -8(%rbp) has two operands %rax and -8(%rbp). These operands are interpreted into Why as

%rax int64 = (integer_of_int64 (sel_int64 !_rax)) -8(%rbp) int64 = (integer_of_int64 (sel_int64 !_8__rbp_))

Annotations

In order to simplify the translation of annotations in this chapter, the annotations handled are only:

• preconditions • post-conditions • assertions
Thanks to the preparation step, all annotations in the original C program are kept in assembly code and all variables in annotations are replaced by registers/memory references.

We denote by A term the translation of a term (logic expression) into Why. The translation of annotations into Why is described as below: There are built-in constructors in ACSL such as \abs, \max, \min, etc. which are overloaded. This means that we can use these constructors for both integer and real type. Here, we decide to simplify it by using ACSL-style syntax. For example: instead of using \abs(e) which returns the absolute value of an integer or a real expression, we use \abs_int(e) and \abs_real(e).

A ==> B term = A term -> B term A == B term = A term = B term A && B term = A term and B term A || B term = A term or B term !A term = not( A term ) #int#v# term = v int32 #long#v# term = v int64 e 1 op
For example, the following annotation

(#int#-4(%rbp)# != 0 ==> #int#-4(%rbp)# == l_sign(\exact(#double#-24(%rbp)#))) && \abs_int(#int#-4(%rbp)#) <= 1 is interpreted into Why as (integer_of_int32(sel_int32(_4__rbp_)) <> 0 -> integer_of_int32(sel_int32(_4__rbp_)) = l_sign((sel_exact(_24__rbp_))
)) and abs_int(integer_of_int32(sel_int32(_4__rbp_))) <= 1

Translation of an instruction

The move instructions and addition/subtraction/multiplication/division instructions are translated thanks to the following abstract functions in Why program:

parameter set_int32_no_check: imm:int -> dest: register ref -> { } unit writes dest { integer_of_int32(sel_int32(dest)) = imm }
The previous abstract function will set a 32-bit integer to a register without verifying if this value overflows or not.

parameter set_int32: imm:int -> dest: register ref -> { is_int32(imm) } unit writes dest { integer_of_int32(sel_int32(dest)) = imm }
The post-condition of set_int32 is the same as set_int32_no_check. Its pre-condition verifies whether imm is a 32-bit integer.

We denote by ins i the Why translation of an instruction ins. Instructions in assembly code are interpreted to Why as follows:

movl src, dest i = set_int32_no_check src int32 dest addl src, dest i = set_int32 ( dest int32 + src int32 ) dest subl src, dest i = set_int32 ( dest int32 -src int32 ) dest imull src, dest i = set_int32 ( dest int32 * src int32 ) dest call label i = label _parameter()
With 64-bit instructions, we declare the following parameters similarly:

parameter set_int64_no_check: imm:int -> dest: register ref -> { } unit writes dest { integer_of_int64(sel_int64(dest)) = imm } parameter set_int64: imm:int -> dest: register ref -> { is_int64(imm) } unit writes dest { integer_of_int64(sel_int64(dest)) = imm }
The translation of 64-bit instructions are as follows: movq src, dest i = set_int64_no_check src int64 dest addq src,

dest i = set_int64 ( dest int64 + src int64 ) dest subq src, dest i = set_int64 ( dest int64 -src int64 ) dest imulq src, dest i = set_int64 ( dest int64 * src int64 ) dest
A special case for the translation is that each assertion is considered as an instruction. If we have an assertion A then its translation is /*@ assert A;*/ i = assert A term ; Notice that leave and ret are instructions in assembly language but they do not have any translation here. Notice also that when a function is called by another one, there are two instructions: pushq %rbp at the beginning of the function and popq %rbp at the end of the function before ret which push and restore the value of %rbp. The memory reference pointed by %rbp is used locally in each function. This means that the addresses of -4(%rbp) in two different functions are different. We do not translate pushq and popq to Why but we can assure that the value of -4(%rbp) in two functions are different.

The source src of the instructions movl(q) is either a constant, a register or a memory reference in 32(64) bits. Therefore, we do not need to verify if it overflows or not. However, for addition/subtraction/multiplication instructions, we need to assure that this computation does not overflow.

The case of division cannot be handled the same way as other operations since the divisor must be check non-null. We thus use idivl src, dest i = div_int32 src int32 dest idivq src, dest i = div_int64 src int64 dest with the special Why parameters:

parameter div_int32: imm: int -> dest: register ref -> { imm <> 0 and is_int32(computer_div(integer_of_int32(sel_int32(dest),imm))) } unit writes dest { integer_of_int32(sel_int32(dest)) = computer_div(integer_of_int32(sel_int32(dest@)),imm) } parameter div_int64: imm: int -> dest: register ref -> {
imm <> 0 and is_int64(computer_div(integer_of_int64(sel_int64(dest),imm))) } unit writes dest { integer_of_int64(sel_int64(dest)) = computer_div(integer_of_int64(sel_int64(dest@)),imm) }

The function computer_div is defined in the Why standard library and denotes the integer division which rounds the result towards 0, which corresponds to the usual convention for division in C and other programming languages.

Sequences and functions

Until now, we have described the translation of annotations and the one of each instruction. How do we translate a sequence of instructions or a function containing annotations into Why?

Assume that we have a function with preconditions, post-conditions and assertions. The translation of this function in assembly code into Why is illustrated in Figure 7.1. As we see in this figure, the post-condition becomes an assertion in Why.

For each function having pre-and post-condition, we define an interface of function in Why (at the bottom right hand side of Figure 7.1) where w is is a set of variables modified in the function.

Notice that the real semantic of assume P f:

-→ let f() = .cfi_startproc /*
[ {} unit reads w { P } ];
where w is a set of variables in P

Translation of annotations in presence of inline function

When the program is compiled with -O2, functions may be inlined. This means that the function called is integrated into the caller. Thanks to the inline assembly statement asm volatile, we can keep the annotations at the place they must be when integrating.

The question now is how to translate the annotations of inline function? Normally, when a function is called, in Why, we will use its interface and the precondition is demanded to prove. In case of inline function, the translation is specified below

• precondition_inline term = assert precondition_inline term • post-condition_inline term = assumes post-condition_inline term
The translation of assertions are not changed in this case.

Soundness of translation

The goal of this section is to demonstrate that if the verification conditions hold then the assembly program respects its annotations.

To express that a program respects its annotations, we use a blocking semantics [START_REF] Herms | A certified multi-prover verification condition generator[END_REF]: the execution of a program will block whenever an invalid annotation is met. Thanks to this definition, to prove that a program respects its annotations, we need to prove that there is a backward simulation [START_REF] Leroy | A formally verified compiler back-end[END_REF] of assembly code by the Why code. We denote (abuse of notation) S(t) the evaluation of a term t in S: it is defined by a natural induction:

S(t 1 + t 2 ) = S(t 1 ) + S(t 2 ) S(t 1 -t 2 ) = S(t 1 ) -S(t 2 ) S(t 1 * t 2 ) = S(t 1 ) * S(t 2 ) S(t 1 /t 2 ) = S(t 1 ) / S(t 2 ) S(!r) = S(r)
The execution of a program in Why is defined by an operational semantics. We use a bigstep-style semantics, that is we define the relation S, P ⇒ m S which means " in state S, the program P executes and terminates after m steps in the state S ".

It is defined by the following rules:

S, r := t ⇒ 1 S[r ← t] S, p 1 ⇒ m S S , p 2 ⇒ n S S, p 1 ; p 2 ⇒ (m+n+1) S A holds in S S, assert A ⇒ 1 S

Notice that if S(A) does not hold, then the program blocks.

There are two forms of Why functions that we use in this part:

• Case 1: Function (with let)

with body,

without parameters, the result is unit.

P re S holds S, body ⇒ n S P ost S holds S, f ⇒ n+1 S • Case 2: Function (with parameter) without body, with parameters: the pure ones x 1 , . . . , x m and the reference ones r 1 , . . . , r l , the result is unit reads u 1 , . . . , u m writes w 1 , . . . , w k , r i , . . . , r l (where w 1 , . . . , w k are global references). S := for all references r such that r / ∈ {w 1 , . . . , w l } then S(r) = S'(r).

P re S {r i =S(s i )} {x j =S(e i )} holds ∀i, j.i = j → s i = s j ∀i, j. s i = u j ∀i, j. s i = w j S ∼ w 1 ,...,
The reason why we need to assure that ∀i, j.i = j → s i = s j and ∀i, j.s i = w j is explained in subsection 7.3.2. Notice again that the execution blocks if a precondition or a postcondition is invalid.

Definition 7.3 A program respects its annotations if it executes without blocking.

In Why, there is a verification condition generator in which a set of Why functions produces a set of formulas (the verification condition). The soundness of this verification condition generator is expressed by the following Theorem:

Theorem 7.4 For all Why program P , if the generated proof obligations hold then P executes without blocking.

In P , the assertions, preconditions and post-conditions are verified when P is executed. Indeed, this theorem was proved by Filliâtre [START_REF] Filliâtre | Preuve de programmes impératifs en théorie des types[END_REF][START_REF] Filliâtre | Verification of non-functional programs using interpretations in type theory[END_REF]. Recently, Herms certified it in Coq [START_REF] Herms | A certified multi-prover verification condition generator[END_REF]. We define a global parameter x. In the parameter inc we modify both y and global reference x. In this example, s 1 = y and w 1 = x. If we call this function as follows:

About the condition in function call

let main() = {...} inc(x) { x = 1 and x = 2}
The Why VCGen rejects it with an error: "Application to x creates an alias". If it was accepted, we would be able to prove the post-condition above which is inconsistent. The reason why we also need to have s i = s j is illustrated by a similar example: • a value (a bit vector bv) from the name of a register, denoted by S(r)

• a value stored in memory at any address x, denoted by S(x).

We denote by bv_to_int32(bv) the value in two's complement integer of a bv, by int32_to_bv(const) the bitvector converted from a two's complement integer in 32 bits const. S 1 , i 2 ; . . . ; i m → S S, i 1 ; . . . ; i m ⇒ S P re S holds S, instr(P ) ⇒ S P ost S holds S, call P ⇒ S where instr(P) is a sequence of instructions from the address P until the instruction ret. For 64-bit arithmetic instructions, their execution is defined similarly. 7.3.4 Relation between the Why state and the assembly state (case of "simple" programs) Definition 7.6 Let bv be a bitvector and reg be a register. bv is congruent to reg (denotes by bv ∼ = reg) if

• bv_to_int32(bv) = reg int32
• bv_to_int64(bv) = reg int64 Definition 7.7 Let S be a Why state and S be an assembly state. S simulates S (denotes by S ∼ S) iff 1. For all register r, S(r) ∼ = S(r) where bv ∼ = reg with bv is bitvector and reg is register.

2. For all memory reference m = off(reg,...), either m is defined and then S(m) ∼ = S(m) or m is not defined.

It is important to notice that in the "simple" case must follow Assumption 7.1, that is two distinguished memory references (syntactically) imply two different addresses. This means that if address(m 1 ) = address(m 2 ) then m 1 = m 2 .

Lemma 1 For all assembly state S, Why state S, if S ∼ S and the translation of A into Why A annot is true in S then A is true in S.

Proof. This is a straightforward induction on the structure of A. 1. With m = 0: We have n = 0.

2. With m ≥ 1: We have n ≥ 1. From S, i 1 , . . . , i n i ⇒ m S
we have S, i 1 i ⇒ p S 1 and S 1 , i 2 , . . . , i n i ⇒ q S with p + q + 1 = m.

We prove below by case analysis on i 1 that there is a S 1 such that S, i 1 ⇒ S 1 where

S 1 ∼ S 1 .
By induction, because q < m, there exists S such that

• S 1 , i 2 , . . . , i n ⇒ S and

• S ∼ S
then S, i 1 , . . . , i n ⇒ S and S ∼ S (proved).

We now proceed our case analysis. We illustrate the proof by the three following cases: movl, addl and call instruction. With others instructions, we do similarly. The translation into Why:

addl src, dest i = set_int32 ( dest int32 + src int32 ) dest
The post-condition of sel_int32 says that integer_of_int32(sel_int32(S (dest))) = integer_of_int32(sel_int32(S(dest)))+integer_of_int32(sel_int32(S(src))) As S ∼ S, we have: integer_of_int32(sel_int32(S (dest)))

= bv_to_int32(S(dest)) + bv_to_int32(S(src)) = (bv_to_int32(S(dest)) + bv_to_int32(S(src))) mod 2 32 (because there is no overflow) = bv_to_int32(S (dest)) (proved).

parameter f _ p a r a m e t e r : _: u n i t-> { i n t e g e r _ o f _ i n t 3 2 ( s e l _ i n t 3 2 ( _ r d i ) ) >= 0 and i n t e g e r _ o f _ i n t 3 2 ( s e l _ i n t 3 2 ( _ r d i ) ) < 100 } u n i t r e a d s _ r d i { t r u e } l e t f_0 ( ) = _LFB0 :

( * #movl %e d i , -20(% r b p ) * ) move_reg32 ! _ r d i _20__rbp_ ; [ { }
u n i t r e a d s _20__rbp_ { i n t e g e r _ o f _ i n t 3 2 ( s e l _ i n t 3 2 (_20__rbp_ ) ) >= 0 and i n t e g e r _ o f _ i n t 3 2 ( s e l _ i n t 3 2 (_20__rbp_ ) ) < 1 0 0 } ] ;

( * #movl -20(% r b p ) , %e a x * ) move_reg32 ! _20__rbp_ _rax ; ( * #movl $100 , %edx * ) move_cte32 ( 1 0 0 ) ( 1 0 0 . 0 ) _rdx ; ( * #movl %edx , %e c x * ) move_reg32 ! _rdx _rcx ; ( * #s u b l %eax , %e c x * )
s e t _ r e g 3 2 ( ( i n t e g e r _ o f _ i n t 3 2 ( s e l _ i n t 3 2 ! _rcx ) ) -( i n t e g e r _ o f _ i n t 3 2 ( s e l _ i n t 3 2 ! _rax ) ) ) _rcx ; ( * #movl %ecx , %e a x * ) move_reg32 ! _rcx _rax ; ( * #movl %eax , -4(%r b p ) * ) move_reg32 ! _rax _4__rbp_ ; a s s e r t { i n t e g e r _ o f _ i n t 3 2 ( s e l _ i n t 3 2 (_4__rbp_ ) ) > 0 } ; a s s e r t { i n t e g e r _ o f _ i n t 3 2 ( s e l _ i n t 3 2 (_4__rbp_ ) ) <= / * @ r e q u i r e s 0 <= x <= 1000 ; @ e n s u r e s \ r e s u l t == x * x ; * / i n t s q u a r e ( i n t x ) { i n t tmp = x * x ; r e t u r n tmp ; } i n t main ( ) { i n t a = 5 ; i n t b = s q u a r e ( a ) ; //@ a s s e r t b== 2 5 ; r e t u r n 0 ; } parameter s q u a r e _ p a r a m e t e r : _: u n i t-> { 0 <= i n t e g e r _ o f _ i n t 3 2 ( s e l _ i n t 3 2 ( _ r d i ) ) <= 1000 } u n i t r e a d s _ r d i { i n t e g e r _ o f _ i n t 3 2 ( s e l _ i n t 3 2 ( _rax ) ) = i n t e g e r _ o f _ i n t 3 2 ( s e l _ i n t 3 2 ( _ r d i ) ) * i n t e g e r _ o f _ i n t 3 2 ( s e l _ i n t

Floating-point programs

Chapter 7 talks about the translation of the program containing only 32-bit and 64-bit integer type. In this chapter, we will extend it with the floating-point computations.

Assembly with floating-point arithmetic

Before entering into the translation, we give some basic knowledge about the different modes: SSE/SSE2, x87 and FMA and their instructions [2,[START_REF]Intel 64 and IA-32 Architectures Software Developer's Mannual[END_REF][START_REF]Intel 64 and IA-32 Architectures Software Developer's Mannual. Volume 2A: Instruction Set Reference[END_REF][START_REF]Intel 64 and IA-32 Architectures Software Developer's Mannual. Volume 2A: Instruction Set Reference[END_REF].

SSE/SSE2

The Intel MMX (MultiMedia eXtensions) technology introduced single-instruction multiple-data (SIMD) capacity into the IA-32 architecture, with the 64-bit mmx registers, 64-bit packed integer data types, and instructions that allowed SIMD operations to be performed on packed integers. SSE extensions expand the SIMD execution model by adding facilities for handling packed and scalar single-precision floating-point value contained in 128-bit registers.

The extension SSE2 is a major enhancement to SSE. It adds new math instructions for double-precision (64-bit) floating-point and also extends mmx instructions to operate on 128-bit xmm registers. The comisd and comiss instructions compare the double-precision (single-precision) floatingpoint values in the low quadwords (doublewords) of first operand and second operand, and sets the ZF, PF, and CF flags in the EFLAGS register according to the result (unordered, greater than, less than, or equal). The unordered result is returned if either source operand is a NaN. The OF, SF and AF flags in the EFLAGS register are set to 0.

Data Transfer Instruction

The comisd instruction differs from the ucomisd instruction in that it signals a SIMD floating-point invalid operation exception when a source operand is either a qNaN or sNaN. The ucomisd instruction signals an invalid numeric exception only if a source operand is an sNaN.

x87 Floating-point Unit

The x87 floating-point unit (FPU) instructions are executed by the processor's x87 FPU. These instructions operate on floating-point, integer and binary-coded decimal(BCD) operands.

x87 FPU registers

This FPU provides several registers. These registers are divided into three groups: data registers, control and status registers, and pointer registers.

The x87 FPU status register has 16 bits. It indicates the current state of the x87 FPU. The flags in the x87 FPU status register include the FPU busy flag, top-of-stack (TOS) pointer, condition code flags, error summary status flag, stack fault flag, and exception flags.

The x87 FPU has 8 floating-point registers to hold the floating-point operands. These registers supply the necessary operands to the floating-point instructions. Unlike the processor's general-purpose registers such as the eax and ebx registers, these registers are organized as a register stack.

R7 R6 R5 R4 R3 R2 R1 R0
Since these registers are organized as a register stack, these names are not statically assigned. That is, st(0) does not refer to a specific register. It refers to whichever register is acting as the top-of-stack (TOS) register. The next register is referred to as st(1), and so on; the last register as st (7). There is a 3-bit top-of-stack pointer in the status register to identify the TOS register. For example, if TOS points to R5 then st(0) indicates R5, st(1) indicates R6, etc. Load operations decrement TOS by one and load a value into the new top-of-stack register, and store operations store the value from the current TOS register in memory and then increment TOS by one.

Each data register can hold an extended-precision floating-point number. It uses the 80 bits format that we mentioned in Chapter 2.

x87 FPU instructions

Most floating-point instructions require one or two operands, located on the x87 FPU dataregister stack or in memory. When an operand is located in a data register, is referenced relative to the st(0) register, rather than by a physical register name. Often the st(0) is an implied operand.

These instructions are divided into the following groups: data transfer, load constants, and FPU control instructions.

Data Transfer Instructions

The data transfer instructions perform the following operations:

• Load a floating-point, integer, or packed BCD operand from memory into the st(0) register.

• Store the value in an st(0) register to memory in floating-point, integer, or packed BCD format.

• Move values between registers in the x87 FPU register stack.

Load Constant Instructions

The following instructions push commonly used constants onto the top st(0) of the x87 FPU register stack:

Each load instruction in the table above follows TOS ← TOS -1.

Basic Arithmetic Instructions The addition, subtraction, multiplication and division instructions operate on the following types of operands:

• Two x87 FPU data registers The src operand of instruction with one operand can be either 32-or 64-bit floating-point number in memory. With instructions with two operands src and dest, both src and dest must be FPU registers.

FMA

gcc uses AVX1 instructions when generating assembly code with option -mfma4. Before talking about FMA instructions we will present some AVX instructions [3].

AVX arithmetic instructions

FMA instructions

Now both AMD and Intel have specifications for FMA. In this document, FMA instructions generated thanks to gcc -mfma4 are specified by AMD [ The implementation of vfmaddsd is presented in Figure 8.1. The destination is a xmm register. When the result is written to the destination xmm register, the upper quadword of the destination register (bits 64 -127) and the upper 128-bits of the corresponding ymm register are cleared to zeros. The intermediate product is not rounded; the infinitely precise product is used in the addition. The result of the addition is rounded.

The implementation of vfnmaddss is in Figure 8.2. The destination is always a xmm register. When the result is written to the destination xmm register, the upper three doublewords of the destination register (bits 32 -127) and the upper 128-bits of the corresponding ymm register are cleared to zeros. The intermediate products are not rounded; the infinitely precise products are used in the addition. The results of the addition are rounded.

Definition of programs supported

This chapter is an extension of the previous one. The programs concerned are the ones in which the types int, long int, float double and long double are supported. An interesting point of this chapter is that we will show the different results obtained by compiling a floating-point program with different options of compiler and different architectures.

The hypothesis of the assembly programs handled in this chapter are

Translation to Why

Abstract functions

The abstract types single, double had already defined in the strict model and full mode and we also used them in Part 1. In the same way, we define abstract type binary80 for x87 mode. Similar to the functions sel_int32 and sel_int64 of Chapter 7, we declare the function for getting the content of a register: logic sel_single : register -> single logic sel_double : register -> double logic sel_80

: register -> binary80 logic sel_exact : register -> real

The logic functions sel_single, sel_double, sel_80 and sel_exact return a single, double, binary80 and an exact value, respectively, from a register.

We denote by opr an operand being register or memory reference; by opr int32 , opr single , opr double and opr binary80 the interpretation of an operand that returns a 32-bit integer value, a real in 32 bits, 64 bits and 80 bits respectively from a register. The translation of operands is specified as follows:

opr single = (single_value (sel_single !opr)) opr double = (double_value (sel_double !opr)) opr binary80 = (binary80_value (sel_80 !opr)) opr exact = (sel_exact !opr)
symbol single = (single_value (sel_single symbol)) symbol double = (double_value (sel_double symbol)) symbol binary80 = (binary80_value (sel_80 symbol)) symbol exact = (sel_exact symbol)

When constants are referenced by %rip

As specified in Intel document [START_REF]Intel 64 and IA-32 Architectures Software Developer's Mannual[END_REF], bytes, words and doublewords in the packed data types are stored in consecutive addresses. The least significant byte, word, or doubleword is stored at the lowest address and the most significant byte, word, or doubleword is stored at the high address.

The ordering of bytes, words, or doublewords in memory is always little endian. That is, the bytes with the low addresses are less significant than the bytes with high addresses.

In assembly language, a floating-point constant is often declared in data section and it is referenced by the special register %rip. For example:

.LC1:

.long 0 .long 1025507328 When a value is stored in memory, we do not know its type. In order to point out the values in Why, we define an axiom in which we presuppose all the types that we can access. Moreover, depending on the number of part stored in the memory, we will determine the size of the number pointed by the symbol.

With the example above, there are two 32-bit number following the symbol .LC1. Assume that the little endian architecture is used, we will consider that from .LC1(%rip), we can get both 32-bit and 64-bit numbers. Thus, in our axiom, we consider both forms: 32-bit integer and single; 64-bit integer and double as follow: logic _LC1__rip_: register axiom _LC1__rip__axiom: integer_of_int64(sel_int64(_LC1__rip_)) = 4404520435568345088 and integer_of_int32(sel_int32(_LC1__rip_)) = 0 and double_value(sel_double(_LC1__rip_))= 0x1p-45 and single_value(sel_single(_LC1__rip_))= 0.0 Note that the integer value 4404520435568345088 is obtained by 1025507328 32 (lower part is 0) and this value is converted to double, that is 2 -45 .

Modifying the translation of general-purpose instructions

As we have mentioned before, data transfer instructions copy data from one operand to another without knowing its type. For example, movl is a data transfer instruction but we do not know whether it transfers a 32-bit integer or a 32-bit floating-point number (float type in C).

In Chapter 7, with the simple program containing only integer value, the instructions movl and movq are considered to copy 32-bit and 64-bit integer value from source to destination. As in programs containing floating-point computations, the data of the mov instructions here may be a floating-point value and the translation of mov in Chapter 7 is not enough anymore.

Our idea is that we define a parameter move_cte32 which copies at the same time a 32-bit integer value, a single value and its exact value to dest and a parameter move_cte64 which copies a 64-bit integer value, a double value and its exact value to dest. We need to assure here that the integer value and the floating-point value are interpreted from the same bitvector. We cannot express this in Why but we did it directly in our translator when generating Why program. It is needed to say that if an instruction movq a b (movl a b) is followed by an instruction movsd b c (movss b c) then in the Why program, the condition integer_of_int64(sel_int(b)) (integer_of_int32(sel_int32(b))) is not used. It is just redundant information. This value is only used when the following instruction is a generalpurpose one.

x87 Floating-point Unit

Representation of the stack in Why

As we mentioned in 8.1.2, the stack has eight floating-point registers (R0 -R7) to hold the floating-point operands. There is a top-of-stack (TOS) pointer which identifies the TOS register.

To represent the stack, one solution is to use a Why array of type register. In order to prove floating-point programs, we use mostly Gappa tool. As Gappa can not use the axioms declared for Why array, we cannot prove Why programs by Gappa in x87 mode. We propose another solution which helps Gappa to prove Why programs in x87 mode: instead of using Why array, we define eight variables of type registers: st0, . . ., st7 to represent for the registers R0-R7 of the stack. When analyzing assembly program, we use a temporary variable st_ptr which simulates the top-of-stack TOS and translate "relative" stack registers to their physical ones. We insist that this temporary variable is not a variable in Why program, it is in our translator. In the generated Why program, there will be only physical stack registers. For translating from "relative" stack registers into physical ones, we set initial value of st_ptr as 8. Before each load instruction is executed, st_ptr decreases by 1. Otherwise, if the stack pops a value (for example fstp, faddp, ect.), st_ptr increases by 1 after this instruction is executed. The value of st_ptr is 0 ≤ st_ptr ≤ 8.

For illustration, Figure 8.3 shows us the changes of the stack after executing the instruction fldl. Assume that the current value of st_ptr is st_ptr = 3. This means that in Why, st3 is on the top of the stack. If we access %st(1), it will be st4 in Why. After executing the instruction fldl, the value of st_ptr decreases by 1 and its value is st_ptr = 2 and now st2 is on the top of the stack. If we want to access %st(i) then the physical register will be st(st_ptr + i).

We denote by st st the translation of st to Why. The translation of the operand -when the stack is used -is specified as follows:

st st
= sti where i = st_ptr st(i) st = sti where i = st_ptr + i and 0 ≤ st_ptr + i ≤ 7

Translation of instructions to Why

Similarly to the parameters for single and double types that we defined in Section 8.3.4, here are the parameters for binary80 type. fsts dest i = set_single st0 binary80 st0 exact dest fstl dest i = set_double st0 binary80 st0 exact dest

Arithmetic instructions

As before, the division is a special case as we have to verify in the precondition that the divisor is non-zero. The parameters for the division operation is declared below: The translation of some arithmetic instructions is specified as follows:

fadds src i = set_80 ( st0 binary80 + src single ) ( st0 exact + src exact ) st0 faddl src i = set_80 ( st0 binary80 + src double ) ( st0 exact + src exact ) st0 fsubs src i = set_80 ( st0 binary80 -src single ) ( st0 exact -src exact ) st0 fsubl src i = set_80 ( st0 binary80 -src double ) ( st0 exact -src exact ) st0 fsubrs src i = set_80 ( src single -st0 binary80 ) ( src exact -st0 exact ) st0 fsubrl src i = set_80 ( src double -st0 binary80 ) ( src exact -st0 exact ) st0 fmuls src i = set_80 ( st0 binary80 * src single ) ( st0 exact * src exact ) st0 fmull src i = set_80 ( st0 binary80 * src double ) ( st0 exact * src exact ) st0 fdivs src i = div_80 src single src exact st0 fdivl src i = div_80 src double src exact st0 fdivrs src i = divr_80 src single src exact st0 fdivrl src i = divr_80 src double src exact st0 fadd %st(i), %st(j) i = set_80 ( stj binary80 + sti binary80 ) ( stj exact + sti exact ) stj faddp %st(i), %st(j) i = set_80 ( stj binary80 + sti binary80 ) ( stj exact + sti exact ) stj fsub %st(i), %st(j) i = set_80 ( stj binary80 -sti binary80 ) ( stj exact -sti exact ) stj fsubp %st(i), %st(j) i = set_80 ( stj binary80 -sti binary80 ) ( stj exact -sti exact ) stj fsubr %st(i), %st(j) i = set_80 ( sti binary80 -stj binary80 ) ( sti exact -stj exact ) stj fsubrp %st(i), %st(j) i = set_80 ( sti binary80 -stj binary80 ) ( sti exact -stj exact ) stj fmul %st(i), %st(j) i
= set_80 ( stj binary80 * sti binary80 ) ( stj exact * sti exact ) stj fmulp %st(i), %st(j) i = set_80 ( stj binary80 * sti binary80 ) ( stj exact * sti exact ) stj fdiv %st(i), %st(j) i = div_80 sti binary80 sti exact stj fdivp %st(i), %st(j) i = div_80 sti binary80 sti exact stj fdivr %st(i), %st(j) i = divr_80 sti binary80 sti exact stj fdivrp %st(i), %st(j) i = divr_80 sti binary80 sti exact stj

AVX instructions

There are AVX instructions that are not FMA ones but they are generated when a program is compiled with gcc -mfma4. The translation of them is described in the following table. Indeed, the specification of these instructions are not different from the instructions in SSE/SSE2. The difference is that they have three operands and have the prefix 'v'. Here are the translations of some AVX instructions:

vaddss src2,

Translation of annotations to Why

Translation of annotations in presence of floating-point arithmetic

The translation of annotations was presented in Chapter 7 but it is only for 32-bit and 64-bit integer. In this chapter, we add some rules of the translation of annotations for floating-point types: float and double.

#float#v# term = single_value(sel_single(v)) #double#v# term = double_value(sel_doule(v)) \exact(#τ #v#) term = sel_exact(v) where τ ∈ {float, double}

When st(i) is in annotations

When a C program is compiled with optimization options in x87 mode, the variables in annotations in assembly code may be a x87 stack register. In this case, their type is no longer a double or float, we impose that its type is extended double.

#single#st(%i)# term = binary80_value(sel_80(sti)) #double#st(%i)# term = binary80_value(sel_80(sti))

For example, the following annotation:

Soundness of translation

The soundness of the translation in simple programs (which contain only general-purpose instructions) is presented in Section 7.3. In this section, we will extend it for floating-point numbers.

Definition of the execution of a assembly program

The 1. For all register r, S(r) ∼ = S(r) where bv ∼ = reg with bv is bitvector and reg is register.

2. For all memory reference m = off(reg,...), either m is defined and then S(m) ∼ = S(m) or m is not defined.

It is necessary to make sure that two distinguished memory references (syntactically) imply two different addresses (This follows Assumption 7.1). This means that if address(m 1 ) = address(m 2 ) then m 1 = m 2 .

We denote double(bv) = double_value(sel_double(reg)), this means that the translation of double(bv) to Why will be double_value(sel_double(reg)). We denote in the same way for single.

Here we will prove that Theorem 8.1 ( the stating of this theorem is the same as Theorem 7.8) is still correct with SSE/SSE2 instructions. Theorem 8.1 For all assembly state S, Why state S such that S ∼ S. For all sequence of assembly instructions i 1 , . . ., i n :

If S, i 1 ; . . . ; i n i ⇒ m S then ∃S such that S, i 1 ; . . . ; i n ⇒ S where S ∼ S .

As specified in parameters, each floating-point number has two parts: the rounded one and the real one (or the exact one). Here, S ∼ S does not depend on sel_exact. The proof of the rounded part is done below. The exact part will be mentioned in the next subsection. As before, we only detail a few instructions. Proof.

movsd instruction: movsd src, dest

In assembly program:

bv2double(S (dest)) = bv2double(S(src)) . Notice that in the translation of addsd instruction, the precondition must specify that the addition does not overflow. If the addition overflows in 64 bits then the program blocks.

The translation into

For proving other instructions in SSE/SSE2 mode and also in x87 and FMA, we do the same way as the proof of the two instructions above.

About exact value

In assembly, there does not exist the real value of an operand. As the assembly code generated contains annotations, we define here the semantic of the assembly in which each operand op (register or memory reference) value is a pair (bv, real ). This semantic is different from the one in Chapter 7 where each operand (register or memory reference) value is a bv. The modified semantic helps us to prove the properties about exact value in annotations.

We denote by S(op) the pair (S b (op), S e (op)) where S b (op) is the notation S(op) we used in subsection 8.5.1, S e (op) is the exact value of op. The execution of an instruction in SSE/SSE2 mode for exact value is defined below: The exact value is an additional information and this value may overflow. However, in order to assure that an instruction is executed without overflow in 32 and 64 bits, the precondition about single and double value of the result is needed. Now we will prove Theorem 8.1 for exact value. Proof.

S,

movsd instruction: movsd src, dest

In assembly program, we have:

S e (dest) = S e (src).
The translation into Why is: movsd src, dest i = set_double_no_check src double src exact dest From the post-condition of set_double_no_check, we have: sel_exact(S e (dest)) = sel_exact(S e (src))

Because S e (reg) = sel_exact(reg) and S ∼ S, we have:

S e (dest) = S e (src)
= S e (dest)(proved).

addsd instruction: addsd src, dest

In assembly program:

S e (dest) = S e (dest) + S e (src).

The translation into Why of this instruction is:

addsd src,dest i = set_double ( src double + dest double ) ( src exact + dest exact ) dest
From the post-condition of set_double we have:

sel_exact(S e (dest)) = sel_exact(S e (dest)) + sel_exact(S e (src))

As S e (reg) = sel_exact(reg) and S ∼ S, we have:

S e (dest) = S e (dest) + S e (src)
= S e (dest)(proved).

Notice that in the translation of addsd instruction, the precondition is necessary in order to make sure that the addition does not overflow. If it overflows then the program blocks.

Others instructions are handled similarly. //@ a s s e r t z == 1 . 0 ; r e t u r n z ; }

Illustrations

The two following examples are used for illustrating the translation of floating-point instructions to Why. These examples are proved with different results when compiling with different options of compiler.

Double rounding example

Let us go back to the example about double rounding in Figure 2.2. We modify it by adding an assertion //@ assert z == 1.0 (See Figure 8.4) which will be true in x87 case and won't be true in SSE2 case.

x87 mode

The assembly code in Figure 8.5 is generated by gcc -S -mfpmath=387.

• At line 1, the instruction movabsq copies the value 4607182418800017408 to %rax. The value 4607182418800017408 corresponds to 1.0 in double. Notice that in assembly code, the floating-point value is represented as an integer value.

• Then at line 2, the content of the register %rax is copied to the memory reference -32(%rbp) (corresponds to x).

• Similarly to line 1, at line 3, 4368493837572636672 (corresponds to 2 -53 + 2 -64 ) to %rax.

• Next, at line 4, -24(%rbp) (corresponds to y) receives the data from register %rax.

• At line 5, loads the the floating-point value at -32(%rbp) to the register st(0)(80 bits) of the stack by the instruction fld. • At line 7, the floating-point value in st(0) is rounded to 64 bits and copied to the memory reference -16(%rbp).

When this assembly code is fed into our translator to Why, and the result analyzed by Why, three proof obligations are produced. One is naturally for proving the assertion, the two others are required to prove the absence of overflow: once at line 6 of Figure 8.5, corresponding to the addition x+y in the source code, and once at line 7, which amounts to store the 80-bit value of the x87 stack into a 64-bit memory cell. These three obligations are proved valid using the Gappa automatic prover (See the screenshot in Figure 8.6).

Here are the VCs for proving the assertion in x87 mode:

_rax: register H1: double_value(sel_double(_rax)) = 0x1.p0 _32__rbp_: register H2: double_value(sel_double(_32__rbp_)) = double_value(sel_double(_rax)) _rax0: register H3: double_value(sel_double(_rax0)) = 0x1.002p-53 _24__rbp_: register H4: double_value(sel_double(_24__rbp_)) = double_value(sel_double(_rax0)) st7: register H5: binary80_value(sel_80(st7)) = double_value(sel_double(_32__rbp_)) H6: no_overflow_binary80(nearest_even, binary80_value(sel_80(st7)) + double_value(sel_double(_24__rbp_))) st7_0: register H7: binary80_value(sel_80(st7_0)) = round_binary80(nearest_even, (binary80_value(sel_80(st7)) + double_value(sel_double(_24__rbp_)))) H8: no_overflow_double(nearest_even, binary80_value(sel_80(st7_0))) _16__rbp_: register H9: double_value(sel_double(_16__rbp_)) = round_double(nearest_even, binary80_value(sel_80(st7_0))) ------------------------------------------------------------------------ Goal : double(_16__rbp_) = 1.0

If we do all possible substitutions, the goal reduces to:

• 64 (• 80 (0x1p0 + 0x1.002p -53)) = 1.0
This is proved automatically by Gappa.

SSE2 mode

The assembly code in Figure 8.7 is generated by gcc -S without any other options. In this code, 64-bit instructions (movabs, movq, movsd, addsd) are used. Thus, all the calculations are in 64 bits.

• From line 1 to line 4, the instructions generated are the same as the ones in Figure 8.5.

• At line 5, %xmm0 receives the floating-point value from -32(%rbp).

• At line 6, a floating-point addition in 64-bit is done by the instruction addsd, that is %xmm0 ← • 64 (%xmm0 + -24(%rbp)).

• Finally, the floating-point value of %xmm0 is rounded to 64 bits and then copied to -16(%rbp) (corresponds to z).

With this assembly code, the generated proof obligation corresponding to the assertion cannot be proved anymore. The modified assertion z == 1.0 + 0x1p-52 can be proved instead. The VCs generated corresponding to the assembly code in Figure 8. ------------------------------------------------------------------------ If we do all possible substitutions, the goal reduces to:

• 64 (0x1p0 + 0x1.002p -53) = 1.0 + 0x1p -52
This is proved by Gappa.

Overflow example

Monniaux [START_REF] Monniaux | The pitfalls of verifying floating-point computations[END_REF] considers the program in Figure 8.8 to illustrate differences between architectures with respect to overflows.

When compiled with non-optimized options

The overflow example is compiled by gcc -mfpmath=387 -O0. Excerpt of the generated assembly code are shown on Figure 8.9. • At line 2, the data in %rax is copied to the memory reference -8(%rbp) (corresponds to the variable x in C program).

• Next, at line 3, st(0) of the stack receives the value of -8(%rbp).

• Then, at line 4, a multiplication operation is done by the instruction fmul, that is st(0) ← • 80 (st(0) * -8(%rbp)).

• At line 5, the value in st(0) is rounded to 64 bits and copied to -16(%rbp) (corresponds to y).

• At line 6, st(0) receives the value of -16(%rbp).

• At line 7, st(0) ← • 80 (st(0) / -8(%rbp)).

• At line 8, -24(%rbp) ← • 64 (st(0)).

• Finally, at line 9, the value of -24(%rbp) is copied to the register %xmm0.

In brief, with the non-optimized assembly code in x87 mode, v * v is calculated in 80 bits (line 4) and then rounded in 64 bits (line 5) before doing the division operation (line 7). The return value of this version is Infinity.

For this version, 5 obligations are generated to check absence of overflow at lines 4, 5, 7 and 8 of assembly code of Figure 8.9, and to check that divisor is not zero at line 7. All are proved by Gappa except the overflow at line 5, where the content of the 80-bit register holding the result of the multiplication is moved into a 64-bit memory cell, which indeed overflows. The VCs of the obligation at line 5 which check the overflow are shown below: ---------------------------------------------------------------------------

When compiled with optimized options

The optimized version (Figure 8.10) is compiled with gcc -mfpmath=387 -O1. In this assembly code, the multiplication and division are done directly in 80 bits (line 3-4). Its result is v and it does not overflow. More precisely:

• At line 1, st(0) receives the value at .LC0(%rip) (the value • 64 (1e308) is put in data section determined by the symbol .LC0).

• At line 2, the stack loads the value at st(0), this means that after executing this instruction, st(1) and st(0) has the same value • 64 (1e308).

• At line 3, st(0) ← • 80 (st(0) * st(1)).

• Then at line 4, st(1) ← • 80 (st(0) / st(1)) and pops st(0).

• At line 5, the result of the division is rounded in 64 bits and put in -8(%rbp).

• The result at line 5 is copied to %xmm0 (at line 6).

For this version, 4 obligations are generated at lines 3, 4 and 5 of Figure 8.10 and all are proved by Gappa. Indeed there is no overflow in this version because the result of multiplication is not temporarily stored into a 64-bit register. The VCs generated for obligation at line 5 are illustrated below: H1: double_value(sel_double(_LC0__rip_)) = 0x1.1ccf385ebc8ap1023 st7: register H2: binary80_value(sel_80(st7)) = double_value(sel_double(_LC0__rip_)) st6: register H3: binary80_value(sel_80(st6)) = binary80_value(sel_80(st7)) H4: no_overflow_binary80(nearest_even, binary80_value(sel_80(st7)) * binary80_value(sel_80(st6))) st6_0: register H5: binary80_value(sel_80(st6_0)) = round_binary80(nearest_even, (binary80_value(sel_80(st7)) * binary80_value(sel_80(st6)))) H6: binary80_value(sel_80(st7)) <> 0.0 and no_overflow_binary80(nearest_even, binary80_value(sel_80(st6_0)) / binary80_value(sel_80(st7))) st7_0: register H7: binary80_value(sel_80(st7_0)) = binary80_value(sel_80(st6_0)) / binary80_value(sel_80(st7)) ----------------------------------------------------------------------no_overflow_double(nearest_even, binary80_value(sel_80(st7_0)))

The hypothesis H1 above comes from the axiom we define for the label .LC0 in data section (line 8-10 of Figure 8.10). The VCs above are explained as follows:

H1 : double(_LC0__rip_) = 0x1.1ccf 385ebc8ap1023 (• 64 (1e308)) H2 : binary80(st7) = double(_LC0__rip_) H3 : binary80(st6) = binary80(st7)
H4 : (binary80(st7) * binary80(st6)) does not overf low in 80 bits H5 : binary80(st6_0) = • 80 (binary80(st7) * binary80(st6)) H6 : binary80(st7) <> 0.0 and (binary80(st6_0)/binary80(st7)) does not overf low in 80 bits H7 : binary80(st7_0) = • 80 (binary80(st6_0)/binary80(st7))

Goal : binary80(st7_0) does not overf low in 80 bits

• 80 (• 80 (0x1.1ccf 385ebc8ap1023 * 0x1.1ccf 385ebc8ap1023)/0x1.1ccf 385ebc8ap1023)
does not overflow in 64 bits. This is proved by Gappa.

With SSE2 mode

We can also analyze the code compiled in the SSE2 mode, resulting in 3 obligations: overflows for the multiplication and division and check divisor is not null. As expected, it cannot be proved that the multiplication does not overflow.

Chapter 9

Handling Conditional and loop statements

The goal of this chapter is to present how to handle conditional and loop statements. In order to do that, we need to construct a control flow graph (CFG) and to translate a CFG to Why.

Conditional instructions in assembly

The conditional instructions are divided into two groups: Jump instructions and conditional move instructions. Both do their task based on the value of the status flags which are set by instructions such as test, cmp, comisd, ucomisd, etc. Jump instructions transfer program control to another point whereas conditional move instructions make a move operation.

Jump instructions

There are two types of jump instructions:

Unconditional jump instruction

The unconditional jump instruction transfers program control to a different point in the instruction stream without recording return information. The destination operand specifies the address of the instruction being jumped to. This operand may be an immediate value, a general-purpose register, or a memory reference. The syntax of unconditional jump instruction is: jmp target

Conditional jump instructions: Jcc

All the Jcc instructions are specified by Intel [START_REF]Intel 64 and IA-32 Architectures Software Developer's Mannual. Volume 2A: Instruction Set Reference[END_REF]. In this chapter, we present some of them which are often appeared in our examples.

We declare 6 status flags of EFLAGS register (See subsection 6.2.2) as follows:

Translation of cmp instruction

The cmp (CoMPare) instruction compares two operands (equal, not equal, and so on) and sets the status flags (See Section 6.2.2). It performs the same operation as the sub except that the result of subtraction is not saved. Thus, cmp does not change the source and destination operands. It is typically used in conjunction with a conditional jump instruction for decision making.

To translate the instruction cmp into Why, we declare the following parameter: The translation of the instruction cmp is as follows: cmpl op1, op2 i = cmp op1 int32 op2 int32 cmpq op1, op2 i = cmp op1 int64 op2 int64

Translation of floating-point comparison instructions

The floating-point comparison instructions are presented in subsection 8.1.1. To translate them into Why, we declare only one parameter comi as follows:

parameter comi: a:real->b:real-> { } unit writes CF, PF, AF, SF, ZF, OF { (b = a -> (ZF <> 0 and CF = 0 and PF = 0 and OF = 0 and SF = 0 and AF = 0)) and (b > a -> (ZF = 0 and CF = 0 and PF = 0 and OF = 0 and SF = 0 and AF = 0)) and (b < a -> (ZF = 0 and CF <> 0 and PF = 0 and OF = 0 and SF = 0 and AF = 0)) }

This parameter compares two real values and set status flags. The translation of this kind of instructions are as below: comisd op1, op2 i = comi op1 double op2 double comiss op1, op2 i = comi op1 single op2 single ucomisd op1, op2 i = comi op1 double op2 double ucomiss op1, op2 i = comi op1 single op2 single / * @ r e q u i r e s e1<= x-\exact ( x ) <= e2 ; @ e n s u r e s ( \ r e s u l t != 0 == > \ r e s u l t == l _ s i g n ( \ exact ( x ) ) ) && @ \ abs ( \ r e s u l t ) <= 1 ; @ * / i n t s i g n ( double x , double e1 , double e2 ) { i f ( x > e2 ) r e t u r n 1 ; i f ( x < e1 ) r e t u r n -1; r e t u r n 0 ; } By convention, the entry node is labeled 0. The exit nodes are those which are not unconditional jumps or do not have a normal successor node (the unconditional jumps do not have a normal successor node but they have a label of node to jump to). We construct one CFG for each function in assembly code.

Example with if statement

In order to understand how to construct a CFG with assembly code, we begin with a simple example in Figure 9.1 which is presented in Section 5.2. This example contains if statements. The assembly code of function sign generated by gcc -S is in Figure 9.2. From Definition 9.1, we construct a CFG for this code (See Figure 9.3). This program has only one function sign, we have thus only one CFG. The table in Figure 9.3 have 14 nodes, begins by node 0 and ends by node 13. Each node contains a label, its content and an optional successor node. Node 13 does not point to any node in order to show that it is the last one. The nodes that are neither the final ones nor unconditional jump ones will point to a successor node in the list. The unconditional jump nodes such as node 6, node 10 do not have a successor node but they have a label of node to jump to. The conditional jump nodes have both a successor node and a node to jump to. For example, node 3 contains the conditional jump instruction je .L2, its successor is node 5 and the node to jump to is node 4. We illustrate this table by a graph (below the table).

Example with do while statement

The goal of the example in Figure 9.4 is to show how to construct a CFG of a program containing a loop statement do while. This example has a function main with two variables x and i. Their initial values are x = 0 and i = 10. In the loop, we increase x by 1 and decrease i by 1. The condition to exit the loop is that i ≤ 0. The precondition of this function is x ≥ 0 and the post-condition ensures x = 10.

In assembly code, the jump condition instructions generated from the statement if in C and the one generated from the statement do while are not different. The one different thing is that the jump instruction of do while creates a cycle in CFG.

The CFG corresponding to the assembly code generated by gcc -S is presented in Figure 9.5. This CFG has 11 nodes from node 0 to node 10. Node 1 contains the preconditions, node 4 contains the invariants and node 9 contains the post-conditions. There is a cycle: 3 → 4 → 5 → 6 → 3.

Example with goto, do while and if statement

The example in Figure 9.6 finds the maximum value of a double array [9]. It contains many complex statements such as goto, if and do while. This CFG is thus a general one. The assembly instruction corresponding to goto statement is just a jmp instruction.

The CFG corresponding to the assembly code generated by gcc -S is presented in Figure 9.7. This CFG has 18 nodes: from node 0 to node 17 and has two cycles: i n t x ; i n t i ; / * @ r e q u i r e s x >= 0 ; @ e n s u r e s x == 1 0 ; * / v o i d main ( ) { x = 0 ; i = 1 0 ; do{ / * @ loop i n v a r i a n t x == 10 -i && 10 >= i > 0 ; @ loop v a r i a n t i ; * / x = x + 1 ; i = i -1 ; } w h i l e ( i >0); }

• 3 → 8 → 6 → 9 → 10 → 11 → 12 → 4 → 7 → 3 • 6 → 9 → 10 → 11 → 12 → 4 → 5 → 6
The node 10 contains invariants of the program. This node is important when we generate Why functions from this CFG. Notice that this example is not proved in this chapter because it contains arrays. It will be proved in Chapter 10.

What we have now is a Control Flow Graph for each function in assembly code. It facilitates the translation to Why which will be explained in the next section.

Translation from a CFG to Why

Our goal now is to build a Why program from a given CFG, so that the VCs generated from that Why program guarantee that the assembly program represented by the CFG satisfies its annotations.

Our construction of that Why program is inspired from other techniques proposed for dealing with unstructured programs in general [5,[START_REF] Filliâtre | Formal Verification of MIX Programs[END_REF].

We assume that on any cycle of the CFG, there is at least one invariant node.

Algorithm 9.2 We start from the initial node and traverse the CFG. The traversal stops whenever we meet a final node or an invariant node. From the assumption above, this traversal must terminate.

/ * @ r e q u i r e s n > 0 && \ v a l i d _ r a n g e ( t , 0 , n -1); @ e n s u r e s \ f o r a l l i n t e g e r k ; 0<=k<n == > \ r e s u l t >= t [ k ] ; @ * / double max_array ( double t [ ] , i n t n ) { double m; i n t i = 0 ; goto L ; The algorithm generateWhy(g:CFG) takes a CFG g and returns a list of Why functions corresponding to g. Assume that we have the following functions:

do{ i f ( t [ i ] > m) { L : m = t [ i ] ; } //@ a s s e r t m >=t [ i ] ; / * @ loop i n v a r i a n t 0 <= i < n && @ \ f o r a l l i n t e g e r k ; 0<=k<=i == > m >= t [ k ] ; @ * / i = i +1; } w h i l e ( i < n ) ;
• why_instr(i): translates a sequence of instructions i (except conditional instructions) into Why

• why_pred(p): translates a predicate p into Why

• why_cond(c): translates a conditional instruction c into Why

• proc_why_func: produces a Why function.

• prefix::str: concatenates a string str to prefix.

• n.succ : returns the successor node of the node n.

This algorithm is described by using the following sub-functions:

• Recursive function traverse_from(n:node):Why expression is invoked on each node n having type precondition or invariant of the CFG. We use the array done in order to ensure that each of such node is treated once. The following function explore will be called in this function.

• Function explore(n:node, pre: Why predicate) generates Why functions to encode a subgraph of the CFG from node n with precondition pre. We use an array visited to know whether a node is visited.

• Recursive function explore_rec(n:node, prefix: Why expression) traverses the CFG from node n and produces the Why function to encode the subgraph starting from n, assuming that prefix is the list of statements encoding the path which arrives to n. The implementation of this function is described as follows:

Firstly, if this node is visited then fails else we set visited[n] by true.

Secondly, we consider the type of node n:

-Instruction i: prefix::why_instr(i) Each e i is the corresponding Why program of p i . This means that if the VCs for the function e i in E are valid then p i satisfies its annotations. Now let us go back to the examples that we have presented in the previous section. By applying Algorithm 9.2, the corresponding Why functions for each CFG is as below:

-
• The CFG in Figure 9.3 does not contains any cycles, so the Why functions generated are:

-0 → 1 → 2 → 3 → 4 → 8 → 9 → 10 → 7 → 12 → 13 -0 → 1 → 2 → 3 → 4 → 8 → 11 → 7 → 12 → 13 -0 → 1 → 2 → 3 → 5 → 6 → 7 → 12 → 13
• For the CFG in Figure 9.5, there is a cycle in which node 4 contains the invariant. The Why functions generated are:

-0 → 1 → 2 → 3 → 4 -4 → 5 → 6 → 3 → 4 -4 → 5 → 7 → 8 → 9 → 10
• The CFG in Figure 9.7 has two cycles with the node 10 containing the invariant. There are 4 Why functions corresponding to this CFG:

-0 → 1 → 2 → 3 → 8 → 6 → 9 → 10 -10 → 11 → 12 → 4 → 5 → 6 → 9 → 10 -10 → 11 → 12 → 4 → 7 → 3 → 8 → 6 → 9 → 10 -10 → 11 → 13 → 14 → 15 → 16 → 17

Examples

Clock drift

This example (See Figure 9.8) is presented in Section 5.4. We prove it with our translator by analyzing its assembly code. The goal of this example is firstly to construct CFG from assembly code generated from a program containing for statement; secondly to apply Algorithm 9.2 to generate Why functions and thirdly to illustrate the value of st_ptr on assembly code generated by gcc -mfpmath=387. We also demonstrate in this example that when compiling with different options of gcc, the bound of the rounding error is changed. Now let us see the assembly code generated by gcc -S -mfpmath=387 -O2 in Figure 9.9. With this assembly code, we look at only instructions related to x87 stack: fldz, flds, etc. which will change the value of top-of-stack. The value of top-of-stack in our model is shown below: \ round_single ( \ nearest_even , 0 . 1 ) == 0 x1 . 9 9 9 9 9 ap -4; * / / * @ r e q u i r e s 0 <= n <= 1 0 0 ; @ e n s u r e s \ abs ( \ r e s u l tn * 0 . 1 ) <= n * C ; * / f l o a t f _ s i n g l e ( i n t n ) { f l o a t t = 0 . 0 f ; i n t i ; f o r ( i =0; i < n ; i ++) { / * @ loop i n v a r i a n t 0 <= i < n <=100 && @ \ abs ( t-i * 0 . 1 ) <= i * C ; @ loop v a r i a n t n-i ; * / L :

Line

//@ a s s e r t 0 . 0 <= t <= 1 0 0 . 0 * ( 0 . 1 +C) ; t = t + 0 . 1 f ; //@ a s s e r t \ abs ( t -(t@L+\round_single ( \ nearest_even ,0.1))) <=B ; } r e t u r n t ; } Remind that we use the intermediate variable st_ptr to store the value of top-of-stack. At the entry of the function, we set st_ptr = 8. For each load instruction, the value of st_ptr decreases by 1. For each store and pop instruction, the value of st_ptr increases by 1. As expected, at the entry and at the exit of the function f_single, st_ptr = 8 indicates that the stack is empty. There is a cycle (line 11 -16) in this assembly code but there are not any instructions which change the value of st_ptr. Obviously, at each node in the CFG corresponding to this function, the value of st_ptr is unique whatever the path of the CFG.

An excerpt of assembly code generated by gcc -mfpmath=387 is in At the entry and at the exit of the cycle, st_ptr = 8 indicates that the stack is empty. In this case, we can also make sure that the value of st_ptr is unique whatever the path.

f _ s i n g l e : .LFB0 :

. c f i _ s t a r t p r o c . . . . f _ s i n g l e : .LFB0 :

t
. c f i _ s t a r t p r o c . . . . . L 4 :

; ; i n v a r i a n t i s h e r e . L 3 : By analyzing assembly code using our translator, we construct a CFG for the function f_single. Then, from this CFG, we have four Why programs by using Algorithm 9.2. The obligations of these four Why functions (which are generated from assembly code compiled with gcc -S and gcc -S -mfpmath=387) are proved completely and automatically by the combination of Gappa, Alt-Ergo and CVC3.

f l d s -12(%rbp) f l d s .LC1 ( %rip ) faddp %st, %st ( 1 )
Thanks to Gappa, we have different values of the bound B when compiling the C program with different options of gcc. The following table show us the value of B depending on the mode and the optimization level. When the program is compiled in SSE2 mode or x87 mode without optimization, the value of B is 4.76838e-07. However, the value of B is much smaller than the previous ones, B = 4.33681e-19, when the C program is compiled with gcc -mfpmath=387 -O2. There, all the value of t at each step is stored in 80 bits stack registers.

KB3D

This example (See Figure 9.11) illustrates the handling of conditional statements, the handling of function calls, and the way we express properties on rounding errors across functions. The presentation has been done in Section 5.2. Our goal here is to analyze what should be the value of E depending on the architecture. Feeding the program in our assembly analyzer in SSE2 mode, the VCs are automatically proved valid using a combination of Gappa and SMT solvers (Alt-Ergo and CVC3). The bound E is indeed in that case exactly the same as the one found in Section 5.2 in a strict IEEE-754 mode. At least on this example, this shows that SSE2 assembly conforms strictly to the standard. The table below shows the value of E that are proved correct using various architecture-dependent settings.

Architecture E SSE2 0x1p-45 x87 -O0 0x1.004p-46 x87 -O2 0x1.004p-46 FMA -O2 0x1.8p-46
As expected, using FMA improves over SSE2 (25% less) since less rounding occur. The extended precision of x87 is even better (around 50% less). The result is the best one with x87 and optimization -O2 where all the values are stored in 80 bits. This new result is 1/1025 the one without optimization. Of course, all these bounds are smaller than the one we found in Section 5.2, which was 0x1.90641p-45 ≈ 3203 × 2 -56 , that is more than 50% higher than the SSE2 one. If we prove this program with Frama-C/Jessie, the value E we can prove is exactly the same as the case SSE2 in the table above.

Discussion

The weakest-precondition computation helps to prove assembly code which is an unstructured programs. However, this approach still have some limits. Compared to the traditional WP, this WP need more information to prove a program, especially when the program has a cycle with a loop invariant. For example, if a constant is stored in a variable and is used in the loop then this #d e f i n e E 0 x1p -45 //@ l o g i c i n t e g e r l _ s i g n ( r e a l x ) = ( x >= 0 . 0 ) ? 1 : -1; / * @ r e q u i r e s e1<= x-\exact ( x ) <= e2 ; @ e n s u r e s ( \ r e s u l t != 0 == > \ r e s u l t == l _ s i g n ( \ exact ( x ) ) ) && @ \ abs ( \ r e s u l t ) <= 1 ; @ * / i n t s i g n ( double x , double e1 , double e2 ) { i f ( x > e2 ) r e t u r n 1 ; i f ( x < e1 ) r e t u r n -1; r e t u r n 0 ; } / * @ r e q u i r e s @ s x == \ exact ( s x ) && s y == \ exact ( s y ) && @ vx == \ exact ( vx ) && vy == \ exact ( vy ) && @ \ abs ( s x ) <= 1 0 0 . 0 && \ abs ( s y ) <= 1 0 0 . 0 && @ \ abs ( vx ) <= 1 . 0 && \ abs ( vy ) <= 1 . 0 ; @ e n s u r e s \ r e s u l t != 0 @ == > \ r e s u l t == l _ s i g n ( \ exact ( s x ) * \ exact ( vx )+\ exact ( s y ) * \ exact ( vy ) ) @ * l _ s i g n ( \ exact ( s x ) * \ exact ( vy )-\ exact ( s y ) * \ exact ( vx ) ) ; * / i n t e p s _ l i n e ( double sx , double sy , double vx , double vy ) { i n t s1 , s 2 ; s 1=s i g n ( s x * vx+s y * vy , -E , E ) ; s 2=s i g n ( s x * vy-s y * vx , -E , E ) ; r e t u r n s 1 * s 2 ; } Figure 9.11: Avionics program information need to be put in invariant. Or if the precondition is needed to prove annotations in the loop then we also put it in invariant.

To illustrate what we have said above, let us go back to example Clock drift in Figure 9.8. When we generate assembly code with gcc -S -mfpmath=387 -O2, the assembly code is in Figure 9.12. Here, the constant • 32 (0.1) is stored in st(1) (line 7 -8). In the cycle (line 10 -18), the addition is done by the instruction fadd %st(1), %st(0), that is %st(0) = • 80 (%st(0) + %st(1)) and the value of %st( 1) is set only once at line 8 (outside the cycle). Hence, the Why functions corresponding to the paths: inv I → . . . → inv I and inv I → . . . → P ostcondition do not have any information about %st [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Std[END_REF]. It is thus necessary to insert the condition %st(1) = • 32 (0.1) into the invariant I.

Still in this example, in the precondition, we have 0 ≤ n ≤ 100. This condition must be added one more time into the invariant I because the Why functions corresponding to the paths: inv I → . . . → inv I and inv I → . . . → P ostcondition do not know that 0 ≤ n ≤ 100.

f _ s i n g l e : .LFB0 :

. c f i _ s t a r t p r o c t e s t l %edi, %edi f l d z j l e . L 2

f l d s .LC1 ( %rip ) f x c h %st ( 1 ) x o r l %eax, %eax . L 4 : . L 3 : #APP / * i n v a r i a n t I * / #NO_APP fadd %st ( 1 )

Handling Arrays and Pointers

With programs containing pointers, the Separation Assumption 7.1 is not correct anymore. In order to prove such programs, we propose a new memory model. New rules for translating instructions and annotations into Why based on the memory model will be presented in the next sections.

Handled programs

The programs supported in this chapter have the following characteristics:

• Arrays and pointers with C types: int, long, float, double and long double are supported,

• Dynamic memory allocation is allowed,

• There are no structure type, no pointer of pointer and no cast in the program.

• The assembly code is generated on IA-64 architecture.

In the previous model, we translated registers and memory references into Why variables. However, there were no differences between the translation of a register and the one of a memory reference. In this model, all the memory references are represented and their translation is similar to the way the program access and store value in the memory.

New rules of translation for operands and instructions

The form of memory references in assembly language has been presented in chapter 6. Based on its form, in this section we will talk about the representation of memory in Why.

Representation of memory in Why

Memory space is divided into segments. Code, data, and stacks may be all contained in the same or different segments depending on the memory model. Each of them is pointed by different registers. For example, local variables are pointed by the register %rbp, a data stored in memory is pointed by %rip, etc. In order to transform the memory references into Why, we model the memory in supposing that the memory is continuous and has three parts:

• Memory references for local variables and function arguments are normally pointed by the register %rbp. This %rbp points to the last data item placed on the stack. The memory references pointed by %rbp are -4(%rbp), -8(%rbp), -16(%rbp), etc. This means that the value of disp is negative.

• Memory references for global variables and constants are pointed by the special register %rip. The displacement disp is usually a symbol or non-negative value.

• Dynamic allocated variables (which are not global variables) occupy a space in memory and their addresses do not change. We reserve a memory space for allocated variables (if exist) and this space is between %rip and %rbp (see the following figure). In other words, we make an assumption: in our model, the address value of allocated variables is less than %rip and greater than %rbp.

We insist that what we suppose here is not exactly the same as what the compiler do. What we assure in our memory model is that the three regions do not overlap. In order to specify the validity of an address, we use the ACSL built-in predicate \valid and \valid_range. In our translator, they are interpreted by: the address is not less than %rbp. The translation of this predicate into Why will be discussed in the subsection 10.2.4.

To model the memory, we declare two following pointers: logic _rbp: int logic _rip: int

We set: _rbp ≤ 0. As _rbp and _rip do not point at the same address in the memory, we will set _rip a positive integer. In our model, we thus set _rip = 4. All the memory references are based on these two pointers. Note that the address of dynamic allocated variables are always not less than _rbp. Two axioms associated with _rbp and _rip are: axiom rbp_axiom: _rbp <= 0 axiom _rip_axiom: _rip = 4 i n t ga1 [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Std[END_REF] , ga2 [ 2 ] , ga3 [ 3 ] ; / * @ e n s u r e s ga1 [ 0 ] == i ; @ e n s u r e s ga2 [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Std[END_REF] == i +1; @ e n s u r e s ga3 [ 2 ] == i +2; @ * / v o i d f g ( i n t i ) { ga1 [ 0 ] = i ; ga2 [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Std[END_REF] = i +1; ga3 [ 2 ] = i +2; } Notice that in assembly code, %rip and %rbp addressing are relative. In our model, %rip and %rbp are considered as constants. This assumption is correct because the proofs are done function by function. Also notice that as we declare _rbp as a constant, all the instructions in which %rbp is the destination operand will not be translated into Why. This happens only when

• at the beginning of the function with instruction movq %rsp, %rbp to copy the value of %rsp to %rbp.

• at the end of the function with popq %rbp to restore the previous value of %rbp.

Without these instructions, the proofs do not change anything because they are done function by function.

We will present the translation of memory references depending on their scope: global variables, local variables and variables being arguments of a function as follows:

Translation of arrays defined as global variables

Remind that the directive .comm symbol , length declares a common symbol. By observing assembly code, we see that the declaration of global variables begins by the assembly directive .comm where symbol is the name of the variable and length is the total size of the variable. If it is an array, length is the memory allocated for this array (the number of elements × size of each element). These symbols are interpreted as constants like %rbp and %rip.

To illustrate the translation of an array being a global variable, we have a small C program in Figure 10.1. Its assembly code compiled with gcc -S is in Figure 10.2.

In this example, we have three integer arrays: ga1 having one element, ga2 having two elements and ga3 having three elements. The size of each element in each array is 4.

In the assembly code, they are declared by using the directive .comm (line 1-3 of Figure 10.2). We consider only the first and the second argument after .comm (the third one is the alignment of the symbol, it is an optional argument). The first one is the symbol name and the second one is the length (in bytes). In others words, the second argument is the total number of bytes allocated for each symbol. These arrays are defined in Why as follows: logic ga1: int logic ga2: int logic ga3: int An important question is: when we have many arrays in the same program, how can we make sure that the memory allocated for each array does not overlap the others?

There exists a model that separates the memory for pointers in C program [START_REF] Filliâtre | Multi-prover verification of C programs[END_REF]. As defined in this separation model, each pointer is considered to have a memory variable. For the three arrays above, we have three different memory variables. This model works well in C code source. We tried to apply this model in assembly code but it does not work. This is because if we use separation model, we cannot define the relation between these memory variables. We thus use only one memory variable for all. This memory variable will be presented in subsection 10.2.3.

As we have explained in subsection 10. We know that the size of int is 4 (bytes) so the total size of the array ga1 is 1 × 4 = 4, the total size of ga2 is 2 × 4 = 8, the total size of ga3 is 3 × 4 = 12. We can also get the total size in the second argument of .comm and we use this value in our translator.

In order to assure that the memory references for these arrays do not overlap, we have some constraints: The order of the symbols in the memory is not important. The important thing is that the memory space for each symbol does not overlap the others.

In Why, we use the following axioms to specify these constraints: axiom ga1_axiom: ga1 = 4 axiom ga2_axiom: ga2 = 8 axiom ga3_axiom: ga3 = 16

Once the axioms above are defined, when we translate these variables into Why, we do not need the pointer _rip anymore. To access ga2 [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Std[END_REF] for example, it is simply written (ga2 + 1*4).

Translation of arrays defined as local variables

An example containing arrays as local variables is in Figure 10.4. Its assembly code is in Figure 10.5.

It is simpler to translate local arrays to Why because in assembly code, they are referenced by _rbp and the memory reference of a certain element of an array has been indicated by the compiler. For example, with the local variable int la2 [2], the first element of la2 in v o i d l g ( i n t i ) { i n t l a 1 [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Std[END_REF] , l a 2 [ 2 ] , l a 3 [ 3 ] ;

l a 1 [ 0 ] = i ; l a 2 [ 1 ] = i + 1 ; l a 3 [ 2 ] = i + 2 ;
//@ a s s e r t l a 1 [ 0 ] == i ; //@ a s s e r t l a 2 [ 1 ] == i + 1 ; //@ a s s e r t l a 3 [ 2 ] == i + 2 ; } assembly code is at -32(%rbp), so in Why, it will be (_rbp -32). la2 [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Std[END_REF] then will be (_rbp + (-32+1*4)) or (_rbp -28) where 4 is the size of int (corresponds to the second operand in line 11 of Figure 10.5).

/ * @ r e q u i r e s \ v a l i d _ r a n g e ( t , 0 , n -1); @ e n s u r e s \ f o r a l l i n t e g e r k ; 0<=k<n == > t [ k ] = 1 . 0 ; @ * / v o i d i n i t ( double t [ ] , i n t n ) { i n t i ; f o r ( i =0; i <n ; i ++) t [ i ] = 1 . 0 ; } Normally, general-purpose registers are used for storing the address of the arguments of a function. We consider their address is a 64-bit integer. In the example above, %rdi (line 5 of Figure 10.7) contains the address of the array t. In order to get this address, in Why, we write (integer_of_int64 (sel_int64 _rdi)).

In the C program, the variable t must be allocated dynamically in the memory before it is used. In the precondition, it is necessary to specify that this pointer is valid by using \valid_range(t,0,n-1). As discussed in Section 10.2.1, there is a memory space reserved for dynamic allocated variables. The translation of \valid_range(t,0,n-1) into Why will be: the address of all elements of t must be greater than or equal to _rbp (See Figure 10.3). We do not specify that t must be less than _rip because the array t may be a global variable.

Definition of memory model

We define here a memory model which will be used in the next section. This model uses the standard theory of arrays [START_REF] Stump | A decision procedure for an extensional theory of arrays[END_REF]. In order to access a value at an address in the memory, we need a variable having the following type: type 'v memory where 'v is a type which has already been defined such as int, real, etc. In our approach, 'v is register.

We also need the following abstract function to access the value in the memory:

logic select: 'v memory, int -> 'v
The logic function select returns a 'v value at an address specified by an integer in the memory 'v memory.

Notice that in the standard theory of arrays, there exists a function store which stores a value to memory. In the previous version, we used both select and store. However, the proof obligations are proved with a big value of timeout. The reason is that the automatic provers take time to find the relation between select and store. In addition, Gappa is unable to use directly the axiomatics about the relation. We found another way to improve this: instead of using store which make the proofs slower, we use only select. The idea here is to indicate in the post-condition of the parameter all the properties we need: express what is changed and what is not changed in the memory after each modification.

Translation of operands and instructions to Why

Type of memory

The address of a pointer or of the first element of an array being argument of a function will be stored in a general-purpose register. As we consider an address to be an integer, there is no difference between the operations (addition, subtraction, multiplication, etc.) on addresses and the ones on integers. The translation of each instruction will be detailed in this section.

The memory is modeled as a map of register. We define the following variable: Remind again that %rbp is a special case of the translation as in our model we defined it as a constant (See Section 10.2.1).

In our model, memory is a map of register and we assume that the size of this register is not fixed (maybe in 32 or in 64 bits depending on the instruction). Therefore, we need some constraints for what is not changed when we update the value at an address. Let us see the figure above. Assume that we update a register with size n bytes at address addr in the memory, there are some different cases:

• If the address is less than or equal to addr-8 (8 is the number of bytes of a double or 64-bit integer) or is greater than or equal to addr+n then the bitvector (not greater than 64 bits) obtained from this address is not changed if we update the value at addr.

• If the address is not greater than addr-4 (4 is the number of bytes of a single or 32-bit integer) then the single value and 32-bit integer value at this address are not changed if we update the value at addr.

• We do not take into account the case where the address is greater than addr-4 and is less than addr because the bitvector (32 or 64 bits) at this address may be changed.

From the idea above, we have the following predicate:

predicate unchanged_mem(MEM1: register memory,MEM2:register memory, addr:int,nb:int) = (forall i:int. i<= addr-8 or i>=addr+nb -> select(MEM1,i) = select(MEM2,i) ) and (forall i:int. i<=addr-4 -> integer_of_int32(sel_int32(select(MEM1, i))) = integer_of_int32(sel_int32(select(MEM2, i))) and single_value(sel_single(select(MEM1, i))) = single_value(sel_single(select(MEM2, i))) )

where MEM1 and MEM2 are the memory after and before the value at address addr is updated.

Here are the translation of some instructions:

Translation of annotations to Why

The translation of variables in annotations is specified as follows:

#int#reg# term = reg int32 #single#reg# term = reg single #double#reg# term = reg double #intpointer#mem# term = select(MEM, mem addr ) int32 #singlepointer#mem# term = select(MEM, mem addr ) single #doublepointer#mem# term = select(MEM, mem addr ) double #intpointer#reg#[V] term = select(MEM, reg int64 + V int32 *4) int32 #singlepointer#reg#[V] term = select(MEM, reg int64 + V int32 *4) single #doublepointer#reg#[V] term = select(MEM, reg int64 + V int32 *8) double #intpointer#symbol#[V] term = select(MEM,symbol + V int32 *4) int32 #singlepointer#symbol#[V] term = select(MEM,symbol + V int32 *4) single #doublepointer#symbol#[V] term = select(MEM,symbol + V int32 *8) double #intpointer#d(%rbp)#[V] term = select(MEM,_rbp + d + V int32 *4) int32 #singlepointer#d(%rbp)#[V] term = select(MEM,_rbp + d + V int32 *4) single #doublepointer#d(%rbp)#[V] term = select(MEM,_rbp + d + V int32 *8) double
where V is the index of an element in the array, V can be either a register, a symbol or a constant in 32 bits. If V is in 64 bits then we replace V int32 by V int64 . For many compilers, %rbp is used to point to local variables. However, in some cases, %rsp is used instead of %rbp. In this section, we will discuss about this case.

A The translation of the operand containing %rsp will be:

d(%rsp, i, s) addr = _rsp + d + s* i int64
The translation of annotations containing %rsp will be:

#intpointer#d(%rsp)#[V] term = select(MEM,_rsp + d + V int32 *4) int32 #singlepointer#d(%rsp)#[V] term = select(MEM,_rsp + d + V int32 *4) single #doublepointer#d(%rbp)#[V] term = select(MEM,_rsp + d + V int32 *8) double
where V is the index of an element in the array, V can be either a register, a symbol or a constant in 32 bits. If V is in 64 bits then we replace V int32 by V int64 . The figure below is the same as the one in subsection 10.2.1. Before executing the instruction subq $totalbyte, %rsp the position of %rsp and %rbp is the same.

After this subtraction instruction, %rbp = %rsp + totalbytes. This is the reason why in the translation of \valid_range, instead of using _rbp, we use _rsp + totalbytes.

Maximum of an array

The example in Figure 10.8 finds the maximum floating-point value in an array. It was presented in subsection 9.4.3. This first simple example has no floating-point computation but it has an array being input of the function. Therefore, we can illustrate how to apply our memory model to prove it. Moreover, we also show that we can handle many complex statements: goto, do while and if in the same program.

The assembly code of this example is in Figure 10.9. The address of the first element of the array t is stored in -40(%rbp) (line 5). The statement m = t[i] in the source code corresponds to line 27 -32 in assembly code. There, -4(%rbp) (corresponds to i) is copied to %eax (line 27). In line 28, the instruction cltq extends %eax in 32 bits to %rax in 64 bits. Then in line 29, the instruction salq shifts %rax to the right by 3 bits, this means that %rax = %rax * 2 3 where 2 3 = 8 is the number of bytes of a double. In line 30, %rax receives the address of t[i], that is the address of the first element of t plus i*8. Then in line 31, the new value of %rax is the value obtained at the address stored in %rax (old value of %rax). Finally, in the line 32, the value of %rax is copied to -24(%rbp)(corresponds to m).

The CFG of the assembly code in Figure 10.9 was presented in Figure 9.7. By using Algorithm 9.2, the four Why functions are generated: The SSE2 mode, supposed to be strictly compliant with the standard, is worse than FMA and x87, because the roundings in FMA and x87 are slightly more precise.

• 0 → 1 → 2 → 3 → 8 → 6 → 9 → 10 • 10 → 11 → 12 → 4 → 5 → 6 → 9 → 10 • 10 → 11 → 12 → 4 → 7 → 3 → 8 → 6 → 9 → 10
The x87 without optimization is better than the SSE2 mode. The reason is that the value p is calculated in 80 bits and then it is rounded to 64 bits before entering the next step of the loop for.

The case of FMA is even better than the two previous modes because p + x[i] × y[i] is computed with a single rounding.

The improvement with x87 with optimization is impressive: around 2 11 2000 times better. The reason is that optimization makes the value of p stored into the x87 stack thus with extended 80-bits precision for the complete execution of the loop: no intermediate rounding to 64-bit is done.

#d e f i n e NMAX 1000 #d e f i n e NMAXR 1 0 0 0 . 0 #d e f i n e B 0 x1 . f57d5p -44 / * @ lemma bound_int_to_ r e a l : @ \ f o r a l l i n t e g e r i ; i <= NMAX == > i <= NMAXR; * / / * @ lemma t r i a n g u l a r _ i n e q u a l i t y : @ \ f o r a l l r e a l x , y , z ; \ abs ( x-z ) <= \ abs ( x-y ) + \ abs ( y-z ) ; * / / * @ r e q u i r e s 0 <= n <= NMAX; @ r e q u i r e s \ v a l i d _ r a n g e ( x , 0 , n-1) && \ v a l i d _ r a n g e ( y , 0 , n-1) ; @ r e q u i r e s \ f o r a l l i n t e g e r i ;

0 <= i < n == > @ \ abs ( x [ i ] ) <= 1 . 0 && \ abs ( y [ i ] ) <= 1 . 0 && @ x [ i ] == \ exact ( x [ i ] ) && y [ i ] ==\exact ( y [ i ] )
; @ e n s u r e s \ abs ( \ r e s u l t -\ exact ( \ r e s u l t ) ) <= n * B ; * / double s c a l a r _ p r o d u c t ( double x [ ] , double y [ ] , i n t n ) { double p = 0 . 0 ; / * @ loop i n v a r i a n t 0 <= i <= n ; @ loop i n v a r i a n t \ abs ( \ exact ( p ) ) <= i ; @ loop i n v a r i a n t \ abs ( p

-\ exact ( p ) ) <= i * B ; @ loop v a r i a n t n-i ; * / f o r ( i n t i =0; i < n ; i ++) { //

Bit-level reasoning

There are numerical programs containing bit-level operations such as negation by XOR, conversion of an integer to a double which were not handled in previous chapters. Finding a model which we can integrate into our translator to prove such programs automatically is a difficult task. Hence, this is a future work. What we can do now is to write manually such programs directly in Why 31 [10]. The goal of this step is to create a bit-level model. Firstly, we will explain why we handle these examples. Secondly, we will present the Why3. Then we will detail how to construct theories in Why 3. Finally, there are examples for illustrating our theories.

Motivations

Before talking about the goal of this chapter, we present firstly some examples which are not proved by our hardware-dependent approach. These examples containing bit-level operations and they are what we want to prove in this chapter.

Examples of the chapter

Negation by XOR

This example is the computation of the negation of a floating-point number using XOR operation which is often found in assembly code. A C program is shown in Figure 11.1 and its assembly code corresponding to the negation is in Figure 11 In this assembly code, %xmm1 stores the value of x (line 1), %xmm0 contains the value 0x8000000000000000 (line 2) which is represented in 64-bit binary format as 100000. . .0. The XOR operation is executed by the instruction xorpd (line 3). Return to our hardware-dependent approach, as we do not handle instructions in bit level, this program cannot be proved. In order to compute the negation of a floating-point number x:

• x XOR 0x80000000 (if x is a single) and • x XOR 0x8000000000000000 (if x is a double).
It is illustrated in Figure 11.3.

Conversion of an integer to double

This example is presented by a C program in Figure 11.4. Notice that this example is tested in a little endian architecture. The definition of endianness and the difference between little endian and big endian will be talked in the next section. The idea of this example is to convert an integer to a double, shown in Figure 11.5. We have two integers: j = 0x43300000 and j' = 0x80000000. The value of Const is obtained by concatenating j and j'. The value of Const is Const = 2 52 + 2 31 . Var is the concatenation of j and (j' XOR x). The value of Var is Var = 2 52 + 2 31 +x. Then, Var -Const is the value of x in double.

We cannot prove the examples above with our hardware-dependent approach because this approach did not handle bit-level operations. If we change the translation in our approach in order to prove them, this means that the translation of each instructions must be done in bitvectors, all the examples from Chapter 7 to Chapter 10 will not be proved anymore.

Goals of the chapter

The goal of this chapter is to do an experience in Why 3: We will implement examples directly in Why 3 and see what is possible to prove with bitvectors.

Also notice that in this chapter, we will use Why 3 instead of Why 22 as in Why 3, it is easy to prove a part of obligations in Coq and with automatic provers. Moreover, in Why 3, we can use or clone a theory, this helps to write Why 3 code clearer.

There does not exist the Why 3 standard library for bitvector. Writing a bitvector library in Why 3 is also a goal of this chapter.

About endianness

Endianness3 is the format indicating how multi-byte data is stored in computer memory. It describes the location of the most significant byte (MSB) and least significant byte (LSB) of an address in memory.

There are two types of endianness:

• Little endian: stores the LSB at the lowest address. The architectures supporting the little endian are Intel 80x86, DEC Alpha, etc.

• Big endian: stores the MSB at the lowest address. The architectures using this type of endianness are Motorola 68k, IBM Power, etc.

Besides, there are architectures such as ARM, PowerPC, etc. which feature a setting which allow to switch to either big or little endian (Bi-Endian) by setting a processor register.

Here is an example illustrating the difference between little endian and big endian: consider a 32-bit integer (in hex) 0xabcdef12. It consists of 4 bytes: ab, cd, ef and 12. Assuming that we store this 32-bit integer at the memory address starting 1000. Then, on little and big endian system, the memory will be: In big endian architecture, we also swap the value of i[0] and i [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Std[END_REF] of the variable Var.

We do not intend to handle the endianness in this chapter. However, it is necessary to know a little about it because it is related to the the bitvector. We consider in the next sections only little endian. The big endian is handled similarly.

Why 3

Why 3 is the next generation of Why platform. It clearly separates the purely logical specification part from generation of verification conditions for programs. In this section, we present logic declarations and theories in Why 3 that we will use in the next sections.

Logic

The logic in Why 3 is a first-order logic with polymorphic types with several extensions: recursive definitions, algebraic data types and inductive predicates.

Theory BitVector

This theory consists of elements and functions of a bitvector. The detail of this theory is presented in Figure 11.8. It contains:

• the size of a bitvector

• an abstract type bv

• the functions about logic operations: bitwise and, bitwise or, bitwise xor, bitwise not, logical shift right/left;

• the function to_nat_sub b j i: returns the non-negative integer value of the bitvector b from index i to j;

• the function from_int n: returns the bitvector representing an integer n on size bits;

• the function from_int2c: returns the bitvector of an 2-complement integer.

We present below how to define the three functions that we will use in our examples in Section 11.5.

Function to_nat_sub

The function to_nat_sub(b:bv, j:int, i:int):int is defined recursively by: to_nat_sub The function from_int(n:int):bv returns the bitvector representing the non-negative integer n on size bits. It is defined by: For all integers n, i such as 0 ≤ i < size

b j i =    to_nat_sub b (j -1) i if 0 ≤ i ≤ j < size and nth b j = False 2 j-i + to_nat_sub b (j -1) i if 0 ≤ i ≤ j <
• If (n/2 i )%2 = 0 then nth (from_int n) i = False • If (n/2 i )%2 <> 0 then nth (from_int n) i = True
The corresponding Why function is shown below:

function from_int (n:int) : bv axiom nth_from_int_high_even: forall n i:int. size > i >= 0 /\ mod (div n (pow2 i)) 2 = 0 -> nth (from_int n) i = False axiom nth_from_int_high_odd: forall n i:int. size > i >= 0 /\ mod (div n (pow2 i)) 2 <> 0 -> nth (from_int n) i = True

Function from_int2c

The function from_int2c (n:int):bv takes a signed integer n as input and returns a bitvector bv with two's complement representation. It is defined by:

• For sign bit: For all integer n:

-If n ≥ 0 then nth (from_int2c n) (size-1) = False -If n < 0 then nth (from_int2c n) (size-1) = True
• For all integer n, i such as

0 ≤ i < size -1 -If (n/2 i )%2 = 0 then nth (from_int n) i = False -If (n/2 i )%2 <> 0 then nth (from_int n) i = True
The declaration of this function in Why is below:

n i:int. size-1 > i >= 0 /\ mod (div n (pow2 i)) 2 = 0 -> nth (from_int2c n) i = False axiom nth_from_int2c_high_odd: forall n i:int. size-1 > i >= 0 /\ mod (div n (pow2 i)) 2 <> 0 -> nth (from_int2c n) i = True
The logic functions div and mod are the integer division (which rounds down) and the modulo functions. The function pow2 i, which is declared in theory Pow2int, return the integer value of 2 i .

Theory BV32 and BV64

The theory BV64, presented in Figure 11.9 and the theory BV32 is the similar one. The theory BV32(64) declares a theory for bitvector in 32(64) bits. It is a clone version of the theory BitVector with size = 32 [START_REF] Pavlova | Vérification de bytecode et ses application[END_REF]. f o r a l l i : i n t . 0<=i <32 -> BV64 . n t h ( c o n c a t b1 b2 ) i = BV32 . n t h b2 i axiom c o n c a t _ h i g h : f o r a l l b1 b2 : BV32 . bv . f o r a l l i : i n t . 32<=i <64 -> BV64 . n t h ( c o n c a t b1 b2 ) i = BV32 . n t h b1 ( i -32) end

Theory BV_double

The theory BV_double contains declarations for converting a 64-bit bitvector to a double. It uses the theory BV64 (see Figure 11.11). In this theory, we do not take into account special

Conversion of an integer to a double

This example was presented in subsection 11.1.1. We implement it directly in Why 3. In order to do that, we have three parts:

• Declaration of constants j and j'

• Proof of Const = 2 52 + 2 31

• Proof of Var(x) = 2 52 + 2 31 + x

• Final goal: Const -Var(x) = (double)x

The detail of this program is shown below: lemma nth_0_30 : f o r a l l x : i n t . f o r a l l i : i n t . i s _ i n t 3 2 ( x ) /\ 0<=i <=30 -> BV32 . n t h ( BV32 . bw_xor j ' ( BV32 . f r o m _ i n t 2 c x ) ) i = BV32 . n t h ( BV32 . f r o m _ i n t 2 c x ) i lemma nth_jpxor_0_30 :

f o r a l l x : i n t . f o r a l l i : i n t . i s _ i n t 3 2 ( x ) /\ 0<=i <=30 -> BV32 . n t h ( j p x o r x ) i = BV32 . n t h ( BV32 . f r o m _ i n t 2 c x ) i lemma nth_var31 : f o r a l l x : i n t .

Discussion

We had tried to define a bitvector model in Why 2 and tried to prove the two examples above. As Why 3 allows proving a part of obligations in Coq and in addition, it allows reusing existing theories, we choose Why 3 for this experience. What we obtained in this chapter are that firstly we defined a bitvector library in Why 3. Secondly, we are successful in proving programs containing bit-level operation in Why 3. With this experience, we had to prove many lemmas in Coq. It is thus too difficult to prove such program only with automatic provers. For this reason, in order to integrate the bitvector model into hardware-dependent approach, the improvement of bitvector theory is needed.

Chapter 12

Conclusion and Future works

Summary

This thesis demonstrates that the formal proofs of numerical programs in considering architecture and compiler aspect is possible. This possibility is shown by two proposed approaches:

Hardware-independent approach

This approach gives correct rounding errors whatever the architecture and allowing many choices to the compiler. This is implemented in the Jessie plugin of the Frama-C framework for all basic operations: addition, subtraction (with possible reordering), multiplication, division, square root, negation, absolute value.

Moreover, it handles both rounding according to 64-bit rounding in IEEE-754 double precision, 80-bit rounding in x87, double rounding in IA-32 architecture, and FMA in Itanium and PowerPC processors and all possible reorganizations of additions and subtractions.

A drawback of this approach is that we may only prove rounding errors. There is no way to prove, for example, that a computation is correct (even if it would be correct in all possible roundings and compilations). This means that some subtle floating-point properties may be lost but bounding the final rounding error is usually what is wanted by engineers and this does not appear to be a big flaw.

Note that we only consider double precision numbers as they are the most used. This is easily applied to single precision computations the same way (with single rounding, 80-bit rounding or double rounding). The idea would be to give similar formulas and to provide the basic operations with those post-conditions.

We only handle rounding-to-nearest (ties to even and ties away from zero). The reason is that directed roundings do not suffer from these problems: double rounding gives the correct answer and if some intermediate computations are done in 80-bit precision, the final result is more accurate, but still correct as it is always rounded in the correct direction. When additions are reordered, we may have different results, but all are smaller than the exact result, so the wanted property still holds whatever the order. Only rounding-to-nearest causes discrepancies.

This work is at the boundary between software and hardware for floating-point programs and this aspect of formal verification is very important. Moreover, this work deals both with normal and subnormal numbers, the latter ones being usually dismissed.

Another interesting point is that our error bounds may be used by other tools. As shown here, considering a slightly bigger error bound for each operation suffices to give a correct final error. This means that if Fluctuat [START_REF] Delmas | Towards an industrial use of FLUCTUAT on safety-critical avionics software[END_REF] for example would use them, it would also consider all our cases of hardware and of compilation.

Hardware-dependent approach

Former work on the verification of assembly code are mainly in the context of the so-called proof-carrying-code, where proof obligations for safety (of memory dereferencing, absence of overflow, etc.) are generated on the object code. However these do not consider any behavioral specification language to specify deeper properties than safety. Thus, we believe that what we have presented is the first method being able to prove architecture-and compiler-dependent behavioral properties of floating-point programs. Our approach and our prototype show that proving complex behaviors of floating-point programs is possible. Our approach is implemented by modifying GAS (≈ 10k lines of C code).

The advantage of this method is that we can prove programs based on their assembly code. For this reason, firstly, it is possible to consider the accuracy of the result depending on the architecture and compiler. Secondly, the precision of each operations is also considered in the calculations. Thirdly, the optimizations are taken into account. Finally, there is no reasoning at the bit-level representation. More importantly, we can prove that the same program may give different result when it is compiled by different architecture/options of compiler.

We handled both single, double and extended rounding. In addition, we also handled FMA instructions. Operations on 32-bit and 64-bit integer were also taken into account. We were successful in proving C programs with complex statement such as for, goto, do while, etc. We supported inlining functions which allow us to prove programs compiled with optimization at -O2 level. By supporting this, we could find better results of programs.

We had two different models:

• the simple one that ignores the access of memory, dedicates to proving programs without arrays and pointers.

• the memory model which supports for programs containing pointers and arrays.

We also considered a special class of programs containing bit operations such as the negation of a floating-point number by using XOR, the bit-level conversion (integer number from/to floating-point number). This work is still in progress as it is only an experience in Why 3 for finding a way to prove such programs automatically.

However, it is clearly not mature enough for a non-expert user, because there is a lot of open issues. First, some languages features are not supported, like pointer casts at the C level, and also at the assembly level. Second, we are not able to interpret all the compiler optimizations. Besides, there are still constraints for each class of programs we handle.

Comparison with related works

Before comparing with previous works, it is important to talk about the role of Why. Why is an excellent tool which supports many automatic and interactive provers. Thanks to Why, we can prove a program by using a combination of many provers. It is very difficult to use only one prover to prove all proof obligations and none of provers can support many features. For example, Gappa is a good prover which can prove obligations about floating point behaviors while others cannot do that. Otherwise, many tools such as Alt-Ergo, CVC3, Z3, etc. helps us to prove non-floating-point properties. For the reasons above, using Why is an advantage of our work compared to previous works which support only one prover. About the first approach, we compared it firstly with the strict model and full model. We based on these two models to construct our approach. The difference is that our model consider architecture and compiler aspects while these models support strictly the IEEE-754 standard. Although the bound of each rounding error is higher than the one in previous works, it is true for multiple architecture. In addition, we consider rounding in 64 bits while the strict model and full model support both 32-bit and 64-bit rounding. Moreover, the full model supports also for special values while we prove that they do not happen like in the strict mode.

Fluctuat is a tool aiming at analyzing the numerical precision and stability of complex algorithms. Its aim is to detect automatically a possible catastrophic loss of precision and its source, or else prove its absence. It relies on abstract domains for the estimation of values and errors, based on interval and affine arithmetic. This tool follows only the IEEE-754 standard while we consider the rounding error for multiple architecture and compiler. About the second approach, as said before, we are the first ones who can prove numerical programs considering architecture and compiler issues. This is thus the first difference compared to the related works.

The work of Burdy and Pavlova [START_REF] Burdy | Java bytecode specification and verification[END_REF] defined a specification language for bytecode called BML (Bytecode Modeling Language) and proposed a verification condition generator for Java bytecode which is completely independent from the source code. They also defined a compiler from JML to BML which allows Java to benefit from the JML source specification (requires not to optimize). Our work is a similar one, dedicates for proving C numerical program, and we use ACSL-style specifications and translate the annotations into assembly without defining a compiler for it. We do not need to define another specification language for assembly. More precisely, we simplify what this work did and we focus on floating-point aspect. Moreover, thanks to Why, we can prove assembly program (compiler optimizations are allowed) with automatic provers while the work of Burdy and Pavlova dedicate to the interactive prover Coq.

Another work that we want to compare with is the one of Myreen in 2008. This work proposed to translate assembly code into recursive functions in the HOL4 system and it did not support for floating-point computation. Our work is a similar one but we translate assembly code of floating-point program into Why intermediate language and the proof obligations can be proved with many different provers.

Compared with the default model of Frama-c/Jessie, proving a program on its assembly code is more complicated than proving it at source level. When we prove a C program on its assembly code, in some cases, we need more information in order to prove it with automatic provers. For example, in the case when we use memory model or when we have a loop invariant in the program.

Difficulties

The difficulty I met is my limited knowledge about assembly language. Indeed, there are a lot of assemblers with different syntaxes. It took me a large amount of time to understand different syntaxes and finally I chose AT&T syntax which is supported by gcc. I also took a lot of time to try to write a parser in OCaml for assembly and finally I found that it is a difficult task and not the main goal of the thesis. For this reason, I modified GAS which already had a parser for assembly.

12.2 Future works 12.2.1 Hardware-independent approach Among the many future works are the numerous other possible compiler optimizations. We have looked a little into multiplication reordering, but, due to underflow, the natural formulas are either wrong or unusable. It is very difficult to only modify the one operation formula to handle all possible underflows in a sequence of multiplications. If we had dismissed underflows, this would have been easy, but we are still trying to find a correct and useful solution. We are also interested in distributivity, meaning a ⊗ (b ⊕ c) ←→ (a ⊗ b) ⊕ (a ⊗ c), and in the replacing of the division by the multiplication by the inverse: a b ←→ a ⊗ (1 b) (this is known to be incorrect, but may speed up a lot of computations, such as Gaussian elimination [START_REF] Cosnard | Computation costs versus communication costs in parallel Gaussian elimination[END_REF]). The interaction of all those optimizations with one another should be carefully studied.

Hardware-dependent approach

A simple future work is to consider special values such as NaNs, infinities, etc. A similar work was done by Marché and Ayad [START_REF] Ayad | Multi-prover verification of floating-point programs[END_REF]. We handled only rounding-to-nearest mode which is a default of gcc. Another perspective is to consider other rounding modes.

In this thesis, we had two models: simple one (Chapter 7 to 9) and memory one (Chapter 10). The goal is to use only memory model as the main model for this work. However, we need to find a way to prove program automatically more easily. Until now, we need redundant assertions to help automatic provers to prove it.

An ongoing work is bit-level reasoning. As discussed in Chapter 11, we proved some examples containing bit-level operations in Why 3 such as the negation of a floating-point number, the conversion of an integer number from/to a floating-point number. However, to prove such programs automatically, it is still a hard work. A possible future work is to find a way to do it automatically.

One another work that we would like to do is to see how the assembly code is changed when it is compiled with gcc -funsafe-math-optimizations where the operations are re-organized and may be rewritten and we would like to study what happens when a processor allows to vectorize [START_REF] Monniaux | The pitfalls of verifying floating-point computations[END_REF]. Our approach may work when the program is compiled with this kind of optimization.

An important limitation is that the translation of annotations by using inline assembly may change assembly code such as the compiler may add more move instructions. Thus it is not a safe way. We wonder if the assert statement in C program is a better choice.

Our translator is tested only on Intel-64 architecture with gcc compiler. One of the future work is try to test it in other architectures. For example, instead of generating FMA instructions by gcc on Intel processors, we can generate assembly code directly by a PowerPC or an Intel Itanium and use our translator to prove it. In order to do that, the translator must be modified because the instructions of PowerPC are different from Intel-64. This work is tedious but direct from our work.

A future work that is not related to floating-point arithmetic is to improve the weakestprecondition algorithm for unstructured programs. We took some time to find out another algorithm but we were not successful. We reused proposed techniques [START_REF] Filliâtre | Formal Verification of MIX Programs[END_REF]5] and in some cases we have to add more information for loop invariants.

For longterm perspectives, our ambition is to integrate this approach into a certified compiler. A similar work related to certified compiler was presented by Leroy [START_REF] Leroy | Formal certification of a compiler back-end or: programming a compiler with a proof assistant[END_REF]. Author proposed the formal certification of a compiler from Cminor (a Clike imperative language) to PowerPC assembly code, using the Coq proof assistant both for programming the compiler and for proving its correctness. A future work is to see whether we can either integrate this approach into such a work or reuse it, not dedicate to interactive language like Coq but our goal is for automatic provers. We also intend to improve the robustness and reduce the trusted code base.

  .1) where s is 0 or 1, e is any integer such that emin ≤ e ≤ emax, m is a number represented by a digit string of the form d 0 • d 1 d 2 . . . d p-1 where d i is an integer digit 0 ≤ d i < b (therefore 0 ≤ m < b).

w bits t = p - 1

 1 1)-bit trailing significand field T = d 1 d 2 . . . d p-1 . The leading bit of the significand, d 0 , is implicitly encoded in the biased exponent E. bits E 0 . . . . . . . . . E w-1 d 1 . . . . . . . . . . . . . . . . . . . . . . . . d p-1

Figure 2 . 1 :

 21 Figure 2.1: Bad case for double rounding

Figure 2 . 3 :

 23 Figure 2.3: A more complex program giving different answers depending on the architecture.

Figure 2 . 4 :

 24 Figure 2.4: A program giving different answers depending on the optimization.

Figure 2 . 5 :

 25 Figure 2.5: Why platform

Figure 2 . 6 :

 26 Figure 2.6: Result of a Why program

Figure 2 . 8 :

 28 Figure 2.8: A C program annotated in ACSL

Figure 2 . 9 :

 29 Figure 2.9: Proof of numerical programs in Frama-C/Jessie

Figure 3 . 1 : 1 3. 1 . 2

 31112 Figure 3.1: Rounding error in 64-bit rounding vs. Theorem 3.1

Figure 3 . 2 :

 32 Figure 3.2: Rounding error in 80-bit rounding vs. Theorem 3.1

  n s _ c o r r e c t n e s s : f o r a l l x f , Rabs x <= powerRZ 2 ( 3 5 0 0 0 ) -> ( f = g ap p a_ ro u nd i ng ( r o u n d i n g _ f l o a t roundNE 53 ( 1 0 7 4 ) ) x \/ f = ga p pa _ ro un d in g ( r o u n d i n g _ f l o a t roundNE 64 ( 1 6 4 4 5 ) ) x \/ f = ga p pa _ ro un d in g ( r o u n d i n g _ f l o a t roundNE 53 ( 1 0 7 4 ) ) ( ga p pa _r o un d in g ( r o u n d i n g _ f l o a t roundNE 64 ( 1 6 4 4 5 ) ) x ) \/ f = x ) -> ( powerRZ 2 ( -1022) <= Rabs x -> ( Rabs ( ( f -x ) / x ) <= 2050 * powerRZ 2 ( -64) /\ powerRZ 2 ( -1022) <= Rabs f ) ) /\ ( Rabs x <= powerRZ 2 ( -1022) -> ( Rabs ( f -x ) <= 2049 * powerRZ 2 ( -1086) /\ Rabs f <= powerRZ 2 ( -1 0 2 2 ) ) ) .

Figure 3 . 3 :

 33 Figure 3.3: Coq theorem certifying the correctness of Theorem 3.1

  parameter add_double : m:mode -> x:double -> y:double -> { no_overflow_double(m,double_value(x) + double_value(y)) } double { add_double_post(m,x,y,result) }

  strict IEEE-754 standard only one possible result: a ⊕ b = • 64 (a + b) ↓ FMA and extended registers a ⊕ b is some real such that |a ⊕ b -(a + b)| ≤ ε|a + b| + η ↓ FMA, extended registers and addition reorganization a ⊕ b is some real such that |a ⊕ b -(a + b)| ≤ ε (|a| + |b|) + η

  Theorem 4.2 is the combination of Theorem 3.1 and Theorem 4.1. It is proved by the previous results.

Theorem 4 . 2

 42 Let ε = 2051 • 2 -60 and ε = 2050 × 2 -64 . Let η = 2049 × 2 -1082 and η = 2049 × 2 -1086 . If we define each operation result as any real such that

#Figure 5 . 1 :

 51 Figure 5.1: Double rounding example with ACSL annotation.

Figure 5 . 2 :

 52 Figure 5.2: Avionics program

Figure 5 . 3 :

 53 Figure 5.3: Result of Figure 5.2 program

#Figure 5 . 4 :

 54 Figure 5.4: Summation program

5 5 1 1 4 e - 18 /

 18 /@ lemma r o u n d 0 1 : \ abs ( ( d o u b l e ) 0 . 1 -0 . 1 ) <= A ; #d e f i n e B 1 . 1 2 2 4 2 1 e -15 // B i s a bound o f r o u n d _ e r r o r ( t +( d o u b l e ) 0 . 1 ) f o r 0 <= t <=NMAXR+ 0 . 0 1 #d e f i n e C (B + A)

Figure 5

 5 Figure 5.5: Clock drift program

40 #

 40 pragma J e s s i e F l o a t M o d e l ( m u l t i r o u n d i n g ) #d e f i n e NMAX 1000 #d e f i n e NMAXR 1 0 0 0 . 0 #d e f i n e B 0 x1 . f57d5p -44

Figure 5 . 6 :

 56 Figure 5.6: Scalar product program

Figure 6 . 1 .Figure 6 . 1 :

 6161 Figure 6.1: Step-by-step from C program to Why proof obligations

Figure 6 . 3 :

 63 Figure 6.3: The program of Figure 6.2 after passing the preparation step

Figure 6 . 4 :

 64 Figure 6.4: Assembly code of the example of Figure 6.3 (compiled by gcc -S)

Figure 6 . 5 :

 65 Figure 6.5: Square program

Figure 6 . 6 :

 66 Figure 6.6: Example in Figure 6.5 after preparation step

Figure 6 . 7 :

 67 Figure 6.7: Assembly code of the example of Figure 6.6 (compiled by gcc -S)

Figure 7 .Figure 7 . 2 :

 772 Figure 7.2: Soundness of the translation

  parameter x: int ref parameter inc: y:int ref -> {...} unit writes x, y { y = 1 and x = 2 }

  parameter inc: x:int ref -> y:int ref -> {...} unit writes x, y { y = 1 and x = 2 } parameter z: int ref let main() = {...} inc(z,z); { z = 1 and z = 2} 7.3.3 Definition of the execution of an assembly program Definition 7.5 (Memory state in assembly program) A memory state S is a map which returns

  For example: bv_to_int32(0x00000001) = 1 bv_to_int32(0xffffffff) = -1 The execution of an instruction i is denoted as: S, i ⇒ S . It is defined as follows: S, movl src, dest ⇒ S[dest ← S(src)] bv_to_int32(S(dest)) + bv_to_int32(S(src)) does not overf low S, addl src, dest ⇒ S[dest ← int32_to_bv(bv_to_int32(S(dest)) + bv_to_int32(S(src)))] bv_to_int32(S(dest)) -bv_to_int32(S(src)) does not overf low S, subl src, dest ⇒ S[dest ← int32_to_bv(bv_to_int32(S(dest)) -bv_to_int32(S(src)))] bv_to_int32(S(dest)) * bv_to_int32(S(src)) does not overf low S, imull src, dest ⇒ S[dest ← int32_to_bv(bv_to_int32(S(dest)) * bv_to_int32(S(src)))] bv_to_int32(S(dest)) <> 0 bv_to_int32(S(dest))/bv_to_int32(S(src)) does not overf low S, idivl src, dest ⇒ S[dest ← int32_to_bv(bv_to_int32(S(dest))/bv_to_int32(S(src)))] S, leave ⇒ S S, ret ⇒ S S, i 1 → S 1

  (a) movl instruction: movl src, dest In assembly program, what we have is: bv_to_int32(S (dest)) = bv_to_int32(S(src)). The translation into Why is: movl src, dest i = set_int32_no_check src int32 dest From the post-condition of set_int32_no_check we have: integer_of_int32(sel_int32(S (dest))) = integer_of_int32(sel_int32(S (src))) As bv_to_int32(bv) = integer_of_int32(sel_int32(reg)) and S ∼ S, we have: bv_to_int32(S (dest)) = bv_to_int32(S(src)) = bv_to_int32(S (dest)) (proved). (b) addl instruction: add src, dest In assembly program: bv_to_int32(S (dest)) = (bv_to_int32(S(src)) + bv_to_int32(S(dest))) mod 2 32

Figure 7 . 4 :

 74 Figure 7.4: Why program of Figure 6.4

Figure 7 . 5 :

 75 Figure 7.5: Result of Figure 7.4 program

Figure 7 . 6 :

 76 Figure 7.6: Square program

3 2 (Figure 7 . 7 :

 277 Figure 7.7: Why program of Figure 6.7

Figure 7 . 8 :

 78 Figure 7.8: Result of Figure 7.7 program

  movsd xmm1 xmm2/m64 Move scalar double-precision floating-point value from xmm1 register to xmm2/m64 movsd xmm2/m64 xmm1 Move scalar double-precision floating-point value from xmm2/m64 to xmm1 register movss xmm1 xmm2/m32 Move scalar single-precision floating-point value from xmm1 register to xmm2/m32 movss xmm2/m32 xmm1 Move scalar single-precision floating-point value from xmm2/m32 to xmm1 register These instructions move a scalar double-precision (single-precision) floating-point value from the source operand (first operand) to the destination operand (second operand). The source and destination operands can be xmm registers or 64-bit (32-bit) memory locations. Packed Arithmetic Instructions addsd xmm2/m64, xmm1 Add the low double-precision floating-point value from xmm2/m64 to xmm1 addss xmm2/m32, xmm1 Add the low single-precision floating-point value from xmm2/m32 to xmm1 subsd xmm2/m64, xmm1 Subtracts the low double-precision floating-point values in xmm2/mem64 from xmm1 subss xmm2/m32, xmm1 Subtracts the low single-precision floating-point values in xmm2/mem32 from xmm1 mulsd xmm2/m64, xmm1 Multiply the low double-precision floating-point value in xmm2/mem64 by low double-precision floating-point value in xmm1 mulss xmm2/m32, xmm1 Multiply the low single-precision floating-point value in xmm2/mem32 by low single-precision floating-point value in xmm1 divsd xmm2/m64, xmm1 Divide low double-precision floating-point value in xmm1 by low double-precision floating-point value in xmm2/mem64 divss xmm2/m32, xmm1 Divide low single-precision floating-point value in xmm1 by low single-precision floating-point value in xmm2/mem32 Comparison Instructions comisd xmm2/m64, xmm1 comiss xmm2/m32, xmm1 ucomisd xmm2/m64, xmm1 ucomiss xmm2/m32, xmm1

  ← • 80 (π) fldl2t st(0) ← • 80 (log 2 10) fldl2e st(0) ← • 80 (log 2 e) fldlg2 st(0) ← • 80 (log 10 2) fldln2 st(0) ← • 80 (log e 2)

  vaddss xmm3/mem32, xmm2, xmm1 Add Scalar Single-Precision floating-point vaddsd xmm3/mem64, xmm2, xmm1 Add Scalar Double-Precision floating-point vsubss xmm3/mem32, xmm2, xmm1 Subtract Scalar Single-Precision floating-point vsubsd xmm3/mem64, xmm2, xmm1 Subtract Scalar Double-Precision floating-point vmulss xmm3/mem32, xmm2, xmm1 Multiply Scalar Single-Precision floating-point vmulsd xmm3/mem64, xmm2, xmm1 Multiply Scalar Double-Precision floating-point vdivss xmm3/mem32, xmm2, xmm1 Divide Scalar Single-Precision floating-point vdivsd xmm3/mem64, xmm2, xmm1 Divide Scalar Double-Precision floating-point The first source operand is either an xmm register or a 64-bit memory location and the second source operand is a xmm register. The destination is a third xmm register. Bits [127:64] of the second source operand are copied to bits [127:64] of the destination. Bits [255:128] of the ymm register that correspond to the destination are cleared.

  parameter move_cte32:a:int-> b:real->bexact:real->c:register ref-> { } unit writes c { integer_of_int32(sel_int32(c)) = a and single_value(sel_single(c)) = b and sel_exact(c) = bexact } parameter move_cte64:a:int->b:real->bexact:real->c:register ref-> { } unit writes c { integer_of_int64(sel_int64(c)) = a and double_value(sel_doule(c)) = b and sel_exact(c) = bexact } The translation of movl and movq is shown as below: movl src, dest i = move_cte32 src int32 src single src exact dest movq src, dest i = move_cte64 src int64 src double src exact dest unit writes c { single_value(sel_single(c)) = round_single(nearest_even, single_value(sel_single(b))/single_value(sel_single(a))) and sel_exact(c) = sel_exact(b)/sel_exact(a) } parameter div_double: a:register -> b:register -> c:register ref -> { double_value(sel_double(a))<>0 and no_overflow_double(nearest_even, double_value(sel_double(b))/double_value(sel_double(a))) } unit writes c { double_value(sel_double(c)) = round_double(nearest_even, double_value(sel_double(b))/double_value(sel_double(a))) and sel_exact(c) = sel_exact(b)/sel_exact(a) } The translation of data instructions and arithmetic ones are presented as follows: Data Transfer Instructions

Figure 8 . 3 :

 83 Figure 8.3: Illustration of the stack with instruction fldl

  parameter set_80_no_check: a:real -> aexact:real -> b:register ref -> { } unit writes b { binary80_value(sel_binary80(b)) = a and sel_exact(b) = aexact } parameter set_80: a:real -> aexact:real -> b:register ref -> { no_overflow_binary80(\nearest_even,a) } unit writes b { binary80_value(sel_binary80(b)) = round_binary80(nearest_even,a) and sel_exact(b) = aexact } The interpretation of instructions is presented as follows: Data transfer instructions flds src i = set_80 src single src exact st0 fldl src i = set_80 src double src exact st0 fldlz i = set_80 (0.0) (0.0) st0 fldl1 i = set_80 (1.0) (1.0) st0

  parameter div_80: a:real -> aexact:real -> b:register ref -> { a<>0 and no_overflow_binary80(nearest_even, binary80_value(sel_80(b))/a) } unit writes b { binary80_value(sel_80(b)) = round_binary80(nearest_even, binary80_value(sel_80(b@))/a) and sel_exact(b) = sel_exact(b@)/aexact } parameter divr_80: a:real -> aexact:real-> b:register ref -> { binary80_value(sel_80(b))<>0.0 and no_overflow_binary80(nearest_even, a/binary80_value(sel_80(b))) } unit writes b { binary80_value(sel_80(b)) = round_binary80(nearest_even, a/binary80_value(sel_80(b@))) and sel_exact(b) = aexact/sel_exact(b@) }

  \abs_real(#double#%st#) <= 1000.0 * (1.0 + 0x1.004p-44) is translated into Why as abs_real(binary80_value(sel_80(st0))) <= 1000.0*(1.0+0x1.004p-44)

  Why: movsd src, dest i = set_double_no_check src double src exact dest From the post-condition of set_double_no_check we have: double_value(sel_double(S (dest)))=double_value(sel_double(S(src))) Because bv2double(bv) = double_value(sel_double(reg)) and S ∼ S, we have: bv2double(S (dest)) = bv2double(S(src)) = bv2double(S (dest))(proved).

2 .

 2 addsd instruction: addsd src, dest In assembly program: bv2double(S (dest)) = • 64 (bv2double(S(dest)) + bv2double(S(src))) if the addition does not overflow. The translation into Why of this instruction is: addsd src,dest i = set_double ( src double + dest double ) ( src exact + dest exact ) dest From the post-condition of set_double we have: double_value(sel_double(S (dest))) = round_double(nearest_even,double_value(sel_double(S(dest)))+ double_value(sel_double(S(src)))) As bv2double(bv) = double_value(sel_double(reg)) and S ∼ S, we have: bv2double(S (dest)) = • 64 (bv2double(S(dest)) + bv2double(S(src))) = bv2double(S (dest)) (proved).

double d o u

  b l e r o u n d i n g ( ) { double x = 1 . 0 ; double y = 0 x1p -53 + 0 x1p -64; double z = x + y ;

Figure 8 . 4 :Figure 8 . 5 :

 8485 Figure 8.4: A simple floating-point program

Figure 8 . 6 :

 86 Figure 8.6: Result of Figure 8.4 program

movabsqFigure 8 . 7 :

 87 Figure 8.7: Assembly code in SSE2 mode of Figure 8.4 example

  7 are: _rax: register H1: double_value(sel_double(_rax)) = 0x1.p0 _32__rbp_: register H2: double_value(sel_double(_32__rbp_)) = double_value(sel_double(_rax)) _rax0: register H3: double_value(sel_double(_rax0)) = 0x1.002p-53 _24__rbp_: register H4: double_value(sel_double(_24__rbp_)) = double_value(sel_double(_rax0)) _xmm0: register H5: double_value(sel_double(_xmm0)) = double_value(sel_double(_32__rbp_)) H6: no_overflow_double(nearest_even, double_value(sel_double(_xmm0)) + double_value(sel_double(_24__rbp_))) _xmm0_0: register H7: double_value(sel_double(_xmm0_0)) = round_double(nearest_even, (double_value(sel_double(_xmm0)) + double_value(sel_double(_24__rbp_)))) _16__rbp_: register H8: double_value(sel_double(_16__rbp_)) = double_value(sel_double(_xmm0_0))

  -double_value(sel_double(_16__rbp_)) = 1.0 + 0x1.p-52 The VCs above say that: H1 : double(_rax) = 0x1p0 H2 : double(_32__rbp_) = double(_rax) H3 : double(_rax0) = 0x1.002p -53 H4 : double(_24__rbp_) = double(_rax0) H5 : double(_xmm0) = double(_32__rbp_) H6 : (double(_xmm0) + double(_24__rbp_)) does not overf low in 64 bits H7 : double(_xmm0_0) = • 64 (double(_xmm0) + double(_24__rbp_)) H8 : double(_16__rbp_) = double(_xmm0_0) Goal : double(_16__rbp_) = 1.0 + 0x1p -52

Figure 8

 8 Figure 8.8: Overflow example

Figure 8 . 9 :

 89 Figure 8.9: Non-optimized assembly code of overflow example

Figure 8 . 10 :

 810 Figure 8.10: Optimized assembly of overflow example

  -no_overflow_double(nearest_even, binary80_value(sel_80(st7_0))) They are explained as follows: H1 : double(_rax) = 0x1.1ccf 385ebc8ap1023 (• 64 (1e308)) H2 : double(_8__rbp_) = double(_rax) H3 : binary80(st7) = double(_8__rbp_) H4 : binary80(st7) * double(_8__rbp_) does not overf low in 80 bits H5 : binary80(st7_0) = • 80 (binary80(st7) * double(_8__rbp_)) Goal : binary80(st7_0) does not overf low in 64 bits If we do all possible substitutions, the goal reduces to: • 80 (0x1.1ccf 385ebc8ap1023 * 0x1.1ccf 385ebc8ap1023) does not overflow in 64 bits This obligation is not proved because • 80 (0x1.1ccf 385ebc8ap1023 * 0x1.1ccf 385ebc8ap1023) overflows in 64 bits.

  parameter CF: int ref parameter PF: int ref parameter AF: int ref parameter ZF: int ref parameter SF: int ref parameter OF: int ref Notice that the value of these status flags in Why is either equal to or different 0 (instead of equal to 1).

  parameter cmp: a:int ->b:int -> { } unit writes OF, ZF, SF { (b = a -> (ZF <> 0 and OF = 0 and SF = 0)) and (b > a -> ZF = 0 and OF = 0 and SF = 0) and (b < a -> ZF = 0 and OF = 0 and SF <> 0) }

Figure 9 . 1 :

 91 Figure 9.1: Example with if

Figure 9 . 2 :

 92 Figure 9.2: Assembly code of program in Figure 9.1

Figure 9 Figure 9 . 5 :

 995 Figure 9.4: Program with loop statement

Figure 9 . 6 :Figure 9 . 7 :

 9697 Figure 9.6: Program with loop and goto statement

---

  Assertion p: prefix::assert why_pred(p) Precondition p: prefix::assume why_pred(p) -Post-condition p: prefix::assert why_pred(p) Invariant p: * produce_why_func(prefix::assert why_pred(p)) * traverse_from(n) Unconditional jump instruction with label l: explore_rec (l, prefix) -Jump instruction with condition c and label l: * explore_rec (l, prefix:: assume why_cond(c))

8 #

 8 d e f i n e A 1 . 4 9 0 1 2 e -09 // A i s a bound o f ( f l o a t ) 0 . 1 -0 . 1 #d e f i n e B 4 . 7 6 8 3 8 e -07 // B i s a bound o f r o u n d _ e r r o r ( t +( f l o a t ) 0 . 1 ) f o r 0 <= t <= 1 0 . 0 1 #d e f i n e C (B + A) / * @ lemma r o u n d 0 1 :\ abs ( \ round_single ( \ nearest_even , 0 . 1 ) -0 . 1 ) <= A ; * / / * @ lemma round_01 :

Figure 9 . 8 :

 98 Figure 9.8: Clock drift program

Figure 9 .

 9 10. All the instructions that change the value of top-of-stack is in the cycle. The change of st_ptr is shown as below

Figure 9 . 9 :

 99 Figure 9.9: Assembly code of program in Figure 9.8 (generated by gcc -S -mfpmath=387 -O2)

Figure 9 . 10 :

 910 Figure 9.10: Assembly code of program in Figure 9.8 (generated by gcc -S -mfpmath=387)

Figure 9 . 12 :

 912 Figure 9.12: Assembly code of program in Figure 9.8 (generated by gcc -S -mfpmath=387 -O2)

  dynamic allocated variables are in this region

Figure 10 . 1 :

 101 Figure 10.1: An example containing arrays as global values.

  e n s u r e s #i n t p o i n t e r#%edx#[START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Std[END_REF] == #i n t #-4(%rbp)#+1; e n s u r e s #i n t p o i n t e r#%ecx#[2] == #i n t #-4(%rbp)#+2; * / #NO_APP r e t . c f i _ e n d p r o c . . . .

Figure 10 . 2 :

 102 Figure 10.2: Assembly code in SSE2 mode of Figure 10.1 example

2 . 1 ,

 21 memory references accessed by %rbp and %rip have two different directions. Memory references pointed by these two pointers are always disjoint. Now let us back to the example in Figure 10.1. The three arrays represented in this model are illustrated in Figure 10.3.

Figure 10 . 3 :

 103 Figure 10.3: Memory model

Figure 10 . 4 :

 104 Figure 10.4: C code of a program with arrays defined as local variables

Figure 10 . 5 :

 105 Figure 10.5: Assembly code in SSE2 mode of Figure 10.4 example

Figure 10 . 6 :Figure 10 . 7 :

 106107 Figure 10.6: A C program with an array defined as an argument of a function

  parameter MEM: register memory refRemind that we have already defined logic functions to get a value from a register in Chapter 7 and Chapter 8, we will use it again.A memory reference mem = disp(base, index, scale) is interpreted as the integer address mem addr = base + disp + index × scale. The translation of memory reference is as follows:d(b, i, s) addr = b int64 +d+s* i int64 d(%rbp, i, s) addr = _rbp + d + s* i int64 symbol(%rip) addr = symbol mem int32 = integer_of_int32(sel_int32(select(MEM, mem addr ))) mem int64 = integer_of_int64(sel_int64(select(MEM, mem addr ))) mem single = single_value(sel_single(select(MEM, mem addr ))) mem double = double_value(sel_double(select(MEM, mem addr ))) mem exact = sel_exact(select(MEM, mem addr ))

  parameter store_double: a:real -> aexact:real -> b:int-> { no_overflow_double(nearest_even,a) } unit writes MEM { double_value(sel_double(select(MEM, b))) = round_double(nearest_even,a) and sel_exact(select(MEM,b)) = aexact and unchanged_mem(MEM, MEM@, b, 8) } The translation of arithmetic instructions are as follows.

3 10. 3

 33 \valid_range(#intpointer#p#,i,j) term = forall k:int. i<=k<=j -> integer_of_int64(sel_int64(p)) + k*4 >= _rbp \valid_range(#singlepointer#p#,i,j) term = forall k:int. i<=k<=j -> integer_of_int64(sel_int64(p)) + k*4 >= _rbp \valid_range(#doublepointer#p#,i,j) term = forall k:int. i<=k<=j -> integer_of_int64(sel_int64(p)) + k*8 >= _rbp For example, the translation of the annotation #int#a#[2] == 3 into Why is as below: #int#a#[2] == 3 term = integer_of_int32(sel_int32(select (MEM, a + 2*4))) = When local variables are pointed by %rsp

• 10 →Figure 10 . 9 :

 10109 Figure 10.9: Assembly code in SSE2 mode of Figure 10.8 example

  .

2 .Figure 11 . 1 :

 2111 Figure 11.1: Negation of a floating-point number

Figure 11 . 2 :

 112 Figure 11.2: Excerpt of assembly code in SSE mode of Figure 11.1 example

Figure 11 . 3 :Figure 11 . 4 :

 113114 Figure 11.3: Negation of a double by XOR

= 2 52 + 2 Figure 11 . 5 :

 522115 Figure 11.5: Convert an integer to a double

Figure 11 . 7 :

 117 Figure 11.6: A Why 3 theory example

Figure 11 .

 11 Figure 11.7 contains bitvector theories that we define for this chapter. It illustrates the relations between the theory BitVector and the others. Besides, we also define other theories: Pow2int Pow2real for the exponent of 2 in integer and in real value.

  size and nth b j = True 0 if i > j For example, if we have a 8-bit bitvector b = 10101010 then to_nat_sub b 3 0 = 2 3 + to_nat_sub b 2 0 (nth b 3 = True) = 2 3 + to_nat_sub b 1 0(nth b 2 = False) = 2 3 + 2 1 + to_nat_sub b 0 0 (nth b 1 = True) = 2 3 + 2 1 + to_nat_sub b -1 0 (nth b 0 = False) = 2 3 + 2 1 (0 > -1) = 10The corresponding function4 in Why 3 is presented as follows:functionto_nat_sub bv int int : int (* (to_nat_sub b j i) returns the non-negative integer whose binary repr is b[j..i] *) axiom to_nat_sub_zero : forall b:bv, j i:int. 0 <= i <= j < size -> nth b j = False -> to_nat_sub b j i = to_nat_sub b (j-1) i axiom to_nat_sub_one : forall b:bv, j i:int. 0 <= i <= j < size -> nth b j = True -> to_nat_sub b j i = (pow2 (j-i)) + to_nat_sub b (j-1) i axiom to_nat_sub_high : forall b:bv, j i:int. i > j -> to_nat_sub b j i = 0 Function from_int

  axiom nth_sign_positive: forall n :int. n>=0 -> nth (from_int2c n) (size-1) = False axiom nth_sign_negative: forall n:int. n<0 -> nth (from_int2c n) (size-1) = True axiom nth_from_int2c_high_even: forall

  Figure 11.9: BV64 theory

Figure 11 . 1 x

 111 Figure 11.10: Bv32_64 theory

t h e o r yFigure 11 . 11 :Figure 11 . 12 :

 11111112 Figure 11.11: Bv_double theory

f u n c

  t i o n j : BV32 . bv = BV32 . f r o m _ i n t 0 x43300000 f u n c t i o n j ' : BV32 . bv = BV32 . f r o m _ i n t 0 x80000000 ( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ) ( * d e f i n i t i o n s : * ) ( * c o n s t : bv64 = c o n c a t j j ' * ) ( * c o n s t _ a s _ d o u b l e : r e a l = double_of_bv64 ( c o n s t ) * ) ( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ) f u n c t i o n c o n s t : BV64 . bv = BV32_64 . c o n c a t j j ' f u n c t i o n c o n s t _ a s _ d o u b l e : r e a l = double_of_bv64 c o n s t ( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ) ( * n e x t lemma : c o n s t _ a s _ d o u b l e = 2^52 + 2^31 * ) ( * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

31 Figure 11 . 14 : 31 (

 31111431 Figure 11.14: Declaration and proofs of Const

BV32 . n tFigure 11 . 16 :

 1116 Figure 11.16: Proofs of lemma1

  Why as the output of our translator. This means that we will translate the assembly code into programs in Why. These Why programs are then the inputs of VC generator to generate VCs.

		ACSL-annotated
		C program
	Frama-C	inline assembly
		translator
		C program
		+ inline assembly
		multirounding
		model
		Assembly code
	Jessie plug-in
		modified version
		of GAS
		Why
		VC generator
	Why platform	verification conditions
	Automatic provers	Interactive provers
	(Alt-Ergo, CVC3,	(Coq, PVS, etc.)
	Gappa, etc.)
		Figure 2.7: Use of Why

The use of Why is illustrated in Figure 2.7. In Part 1, we use Frama-C/Jessie which is a front-end of the Why platform for deductive program verification for proving a numerical program. Here, we define a model written in Why called multirounding. The Jessie plug-in uses this model to interpret each operation in a C program. In Part 2, we use In Chapter 11, Why 3 -the next generation of Why (also called Why 2) -is used. We consider in this chapter that Why 3 is a specification language and we write the program directly in Why 3. The VCs are generated from this Why 3 program by Why 3 platform.

and for an ordering o 1 of the additions, we are able to deduce that |S o 1

	n -n 0 a i | ≤ I.
	Now if we set the addition post-condition as: x ⊕ y is any real number r such that
	|r -(x + y)| ≤ ε • |x + y| + η.
	Then, whatever the ordering o 2 of the additions, we have |S o 2 n -n 0 a i | ≤ I.

  w l ,s 1 ,...,s k

S P ost S {r

i =S (s i )} {x j =S(e i )} {r i @=S(s i )} holds S, f (e 1 , . . . , e n , s 1 , . . . , s k ) ⇒ n+1 S

where S ∼ w 1 ,...,w l

  Theorem 7.8 For all assembly state S, Why state S such that S ∼ S. For all sequence of assembly instructions i 1 , . . ., i This says that if i 1 , . . . , i n executes (without blocking) and we obtain the Why state S then i 1 , . . . , i n executes (without blocking) and we have the final memory state S of the execution of i 1 , . . . , i n so that S ∼ S . Proof. By recurrence on m:

n : If S, i 1 ; . . . ; i n i ⇒ m S then ∃S such that S, i 1 ; . . . ; i n ⇒ S where S ∼ S .

  movsd src, dest i = set_double_no_check src double src exact dest movss src, dest i = set_single_no_check src single src exact dest Arithmetic Instructions addsd src, dest i = set_double ( src double + dest double ) ( src exact + dest exact ) dest addss src, dest i = set_single ( src single + dest single ) ( src exact + dest exact ) dest subsd src, dest i = set_double ( dest double -src double ) ( dest exact -src exact ) dest subss src, dest i = set_single ( dest single -src single ) ( dest exact -src exact ) dest mulsd src, dest i = set_double ( dest double * src double ) ( dest exact * src exact ) dest mulss src, dest i = set_single ( dest single * src single ) ( dest exact * src exact ) dest

	divsd src, dest i = div_double !src !dest dest
	divss src, dest i	= div_single !src !dest dest

  src1, dest i = set_single ( src1 single + src2 single ) ( src1 exact + src2 exact ) dest vaddsd src2, src1, dest i = set_double ( src1 double + src2 double ) ( src1 exact + src2 exact ) dest vsubss src2, src1, dest i = set_single ( src1 single -src2 single ) ( src1 exact -src2 exact ) dest vsubsd src2, src1, dest i = set_double ( src1 double -src2 double ) ( src1 exact -src2 exact ) dest vmulss src2, src1, dest i = set_single ( src1 single * src2 single ) ( src1 exact * src2 exact ) dest vmulsd src2, src1, dest i = set_double ( src1 double * src2 double ) ( src1 exact * src2 exact ) dest vdivss src2, src1, dest i = div_single !src1 !src2 dest vdivsd src2, src1, dest i = div_double !src1 !src2 dest The translation of FMA instructions is specified as follows: vfmaddss src3,src2,src1,dest i = set_single ( src1 single * src2 single + src3 single ) ( src1 exact * src2 exact + src3 exact ) dest vfmaddsd src3,src2,src1,dest i = set_double ( src1 double * src2 double + src3 double ) ( src1 exact * src2 exact + src3 exact ) dest vfmsubss src3,src2,src1,dest i = set_single ( src1 single * src2 single -src3 single ) ( src1 exact * src2 exact -src3 exact ) dest vfmsubsd src3,src2,src1,dest i = set_double ( src1 double * src2 double -src3 double ) ( src1 exact * src2 exact -src3 exact ) dest vfnmaddss src3,src2,src1,dest i = set_single (-( src1 single * src2 single )+ src3 single ) (-( src1 exact * src2 exact )+ src3 exact ) dest vfnmaddsd src3,src2,src1,dest i = set_double (-( src1 double * src2 double )+ src3 double ) (-( src1 exact * src2 exact )+ src3 exact ) dest vfnmsubss src3,src2,src1,dest i = set_single (-( src1 single * src2 single )-src3 single ) (-( src1 exact * src2 exact )-src3 exact ) dest vfnmsubsd src3,src2,src1,dest i = set_double (-( src1 double * src2 double )-src3 double )

(-( src1 exact * src2 exact )-src3 exact ) dest

  definition of memory state in assembly program in Section 7.3 is reused. We denote by bv2single(bv) and bv2double(bv) the representation in single and double of bv, by single2bv(x) and double2bv(x) the bitvector of a single and double value x, by • 32 (x) and • 64 (x) the rounding value in round-to-nearest mode in 32 and 64 bits of x. The execution of an instruction in SSE/SSE2 mode is : S, i ⇒ S . It is defined as follows: Relation between Why state and assembly state (case of floating-point programs) Remind what we have said in Section 7.3 that if S is a Why state and S is assembly state then S simulates S (denotes by S ∼ S) iff

	bv2double(S(dest)) + bv2double(S(src)) does not overf low in 64 bits
	S, addsd src, dest ⇒ S[dest ← double2bv(• 64 (bv2double(S(dest)) + bv2double(S(src))))]
	bv2single(S(dest)) -bv2single(S(src)) does not overf low in 32 bits
	S, subss src, dest ⇒ S[dest ← single2bv(• 32 (bv2single(S(dest)) -bv2single(S(src))))]
	bv2double(S(dest)) -bv2double(S(src)) does not overf low in 64 bits
	S, subsd src, dest ⇒ S[dest ← double2bv(• 64 (bv2double(S(dest)) -bv2double(S(src))))]
	bv2single(S(dest)) * bv2single(S(src)) does not overf low in 32 bits
	S, mulss src, dest ⇒ S[dest ← single2bv(• 32 (bv2single(S(dest)) * bv2single(S(src))))]
	bv2double(S(dest)) * bv2double(S(src)) does not overf low in 64 bits
	S, mulsd src, dest ⇒ S[dest ← double2bv(• 64 (bv2double(S(dest)) * bv2double(S(src))))]
	bv2single(S(src)) <> 0 bv2single(S(dest))/bv2single(S(src)) does not overf low in 32 bits
	S, divss src, dest ⇒ S[dest ← single2bv(• 32 (bv2single(S(dest))/bv2single(S(src))))]

S, movss src, dest ⇒ S[dest ← single(S(src))] S, movsd src, dest ⇒ S[dest ← double(S(src))] bv2single(S(dest)) + bv2single(S(src)) does not overf low in 32 bits S, addss src, dest ⇒ S[dest ← single2bv(• 32 (bv2single(S(dest)) + bv2single(S(src))))] bv2double(S(src)) <> 0 bv2double(S(dest))/bv2double(S(src)) does not overf low in 64 bits S, divsd src, dest ⇒ S[dest ← double2bv(• 64 (bv2double(S(dest))/bv2double(S(src))))] 8.5.2

  movss src, dest ⇒ S e (dest) ← S e (src) S, movsd src, dest ⇒ S e (dest) ← S e (src) bv2single(S b (dest)) + bv2single(S b (src)) does not overf low in 32 bits S, addss src, dest ⇒ S bv2single(S b (dest)) * bv2single(S b (src)) does not overf low in 32 bits S, mulss src, dest ⇒ S e (dest) ← S e (dest) * S e (src) bv2double(S b (dest)) * bv2double(S b (src)) does not overf low in 64 bits S, mulsd src, dest ⇒ S e (dest) ← S e (dest) * S e (src) bv2single(S b (src)) <> 0 bv2single(S b (dest))/bv2single(S b (src)) does not overf low in 32 bits S, divss src, dest ⇒ S e (dest) ← S e (dest)/S e (src) bv2double(S b (src)) <> 0 bv2double(S b (dest))/bv2double(S b (src)) does not overf low in 64 bits S, divsd src, dest ⇒ S e (dest) ← S e (dest)/S e (src)

e (dest) ← S e (dest) + S e (src) bv2double(S b (dest)) + bv2double(S b (src)) does not overf low in 64 bits S, addsd src, dest ⇒ S e (dest) ← S e (dest) + S e (src) bv2single(S b (dest)) -bv2single(S b (src)) does not overf low in 32 bits S, subss src, dest ⇒ S e (dest) ← S e (dest) -S e (src) bv2double(S b (dest)) -bv2double(S b (src)) does not overf low in 64 bits S, subsd src, dest ⇒ S e (dest) ← S e (dest) -S e (src)

  We present here some generalpurpose and SSE/SSE2 instructions:addl reg, mem i = store_imm32 ( mem int32 + reg int32 ) mem addr addl src, reg i = set_int32 ( reg int32 + src int32 ) reg addq src, reg i = set_int64 ( reg int64 + src int64 ) reg addq reg, mem i = store_imm64 ( mem int64 + reg int64 ) mem addr subl reg, mem i = store_imm32 ( mem int32 -reg int32 ) mem addr subl src, reg i = set_int32 ( reg int32 -src int32 ) reg subq src, reg i = set_int64 ( reg int32 -src int64 ) reg subq reg, mem i = store_imm64 ( mem int64 -reg int64 )mem addr addss reg, mem i = store_single ( mem single + reg single ) ( mem exact + reg exact ) mem addr addss mem, reg i = set_single ( reg single + mem single ) ( reg exact + mem exact ) reg addsd reg, mem i = store_double( mem double + reg double ) ( mem exact + reg exact ) mem addr addsd mem, reg i = set_double ( reg double + mem double ) ( reg exact + mem exact ) reg subss reg, mem i = store_single( mem single -reg single ) ( mem exact -reg exact ) mem addr subss mem, reg i = set_single ( reg single -mem single ) ( reg exact -mem exact ) reg subsd reg, mem i = store_double( mem double -reg double ) ( mem exact -reg exact ) mem addr subsd mem, reg i = set_double ( reg double -mem double ) ( reg exact -mem exact ) reg

  piece of assembly function f which uses %rsp to point to local variables looks like:At the beginning of f, %rsp is decreased by n = 32 bytes where n is the total bytes of local variables of the function. After subtracting n from %rsp, the displacement of memory reference pointed by %rsp is a non-negative integer. Then, the value of parameters of the function is copied to local variables pointed by %rsp. The pre-condition generated is put right after these data transfer instructions. At the end of the function f, the value of %rsp is restored by adding with n. Now let us go back to the memory model, we will declare the Why variable corresponding to %rsp as follows:parameter _rsp: int ref and add to each Why function generated from this function a condition assume _rsp <= 0.

	f:
	subq $32, %rsp
	....
	movsd %xmm0, (%rsp)
	movsd %xmm1, 8(%rsp)
	movsd %xmm2, 16(%rsp)
	movsd %xmm3, 24(%rsp)
	(preconditions)
	....
	addq $32, %rsp

  Why variable. It is imm in the instruction subq $imm, %rsp of the current function. If the function does not have this instruction then totalbytes = 0. Also notice that in the same program, the value of totalbytes is changed from function to function. Our translator will replace totalbytes by a constant in the translation \valid_range into Why.

	\valid_range(#intpointer#p#,i,j) term	= forall k:int. i<=k<=j ->
	integer_of_int64(sel_int64(p)) + k*4 >= _rsp + totalbytes
	\valid_range(#singlepointer#p#,i,j) term = forall k:int. i<=k<=j ->
	integer_of_int64(sel_int64(p)) + k*4 >= _rsp + totalbytes
	\valid_range(#doublepointer#p#,i,j) term = forall k:int. i<=k<=j ->
	integer_of_int64(sel_int64(p)) + k*8 >= _rsp) + totalbytes
	where totalbytes is not a The Why function corresponding to the function f above is:
	let f =
	assume _rsp <=0
	subq $32, %rsp i
	....
	movsd %xmm0, (%rsp) i
	movsd %xmm1, 8(%rsp) i
	movsd %xmm2, 16(%rsp) i
	movsd %xmm3, 24(%rsp) i
	assume preconditions term
	....
	addq $32, %rsp i
	void

  Let us go back to example in Figure11.4, now we assume that the memory address of Const starting at 1000 and with both little and big endian, the memory of Const will be:With little endian, the value of Const = 2 52 + 231 . However, the value of Const with big endian does not. In order to get the same value Const = 2 52 + 231 with big endian, we just swap the value of i[0] and i[START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Std[END_REF], the memory is as below:

		Address Little endian Big endian
		1000		12		ab
		1001		ef		cd
		1002		cd		ef
		1003		ab		12
			i[1]			i[0]
	Little endian	43	30	00	00	80	00	00	00
	Big endian	00	00	30	43	00	00	00	80
	Address	1007 1006 1005 1004 1003 1002 1001 1000
			i[1]				i[0]
	Big endian	00	00	00	80	00	00	30	43
	Address	1007 1006 1005 1004 1003 1002 1001 1000

  * * ) lemma n t h _ c o n s t 1 : f o r a l l i : i n t . 0 <= i <= 30 -> BV64 . n t h c o n s t i = F a l s e lemma n t h _ c o n s t 2 : BV64 . n t h c o n s t 31 = True lemma n t h _ c o n s t 3 : f o r a l l i : i n t . 32 <= i <= 51 -> BV64 . n t h c o n s t i = F a l s e lemma n t h _ c o n s t 4 : f o r a l l i : i n t . 52 <= i <= 53 -> BV64 . n t h c o n s t i = True lemma n t h _ c o n s t 5 : f o r a l l i : i n t . 54 <= i <= 55 -> BV64 . n t h c o n s t i = F a l s e lemma n t h _ c o n s t 6 : f o r a l l i : i n t . 56 <= i <= 57 -> BV64 . n t h c o n s t i = True lemma n t h _ c o n s t 7 : f o r a l l i : i n t . 58 <= i <= 61 -> BV64 . n t h c o n s t i = F a l s e lemma n t h _ c o n s t 8 : BV64 . n t h c o n s t 62 = True

SRI stands for Système de Référence Inertielle or Inertial Reference System.

http://java.sun.com/docs/books/jls/third_edition/html/expressions.html#249198

http://frama-c.com/

This approach is integrated in the version 2.30 of Why, available at http://why.lri.fr/download/why-2. 30.tar.gz

These tools are available at http://www.lri.fr/~nguyen/contrib/

http://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html

http://why.lri.fr/

http://frama-c.com/acsl.html

A monotonic function f is a function such that, for all x and y, x ≤ y implies f (x) ≤ f (y).

The Gappa Coq Library adds a Gappa tactic to the Coq Proof Assistant. This tactic invokes the Gappa tool to solve properties about floating-point or fixed-point arithmetic. It can also solve simple inequalities over real numbers.

The proof is available at http://www.lri.fr/~nguyen/research/rnd_64_80_post.html

See also http://research.nianet.org/fm-at-nia/KB3D/.

http://autarkaw.wordpress.com/2008/06/02/round-off-errors-and-the-patriot-missile/

This example can also be found at http://hisseo.saclay.inria.fr/drift.html

Netwide Assembler (NASM): http://www.nasm.us/

Microsoft Macro Assembler (MASM): http://www.masm32.com/

GNU Assembler (GAS): http://sourceware.org/binutils/docs-2.21/as/index.html

In Intel 64 architecture, ESP, EBP and EIP are replaced by RSP, RBP and RIP, respectively

See http://gcc.gnu.org/onlinedocs/gcc/Simple-Constraints.html

http://www.cs.berkeley.edu/~necula/cil/

Advanced Vector Extensions (AVX) is an extension to the x86 instruction set architecture for microprocessors from Intel and AMD proposed by Intel in March

http://why3.lri.fr/

The Why syntax that we use from Chapter 2 to Chapter 10 is called Why 2

http://en.wikipedia.org/wiki/Endianness http://www.intel.com/design/intarch/papers/endian.pdf

Notice that Why 3 does not allow defining functions recursively on integers, this is the reason why we use axioms instead of defining recursive functions.

Acknowledgements

This thesis is the result of more than three years of work under the supervision of Sylvie Boldo and Claude Marché, to who I would like to express my sincere gratitude. They gave me a chance to work in a big team. I would like to thank to Sylvie for her helps from the beginning of my PhD. She was always available to discuss with me each time I needed her. I would like to thank to Claude for his guidance. Spending time working with him gave me more ideas to finish my thesis.

/ * @ r e q u i r e s n >= 0 && n < 1 0 0 ; * / i n t f ( i n t n ) { i n t tmp = 100 -n ; //@ a s s e r t tmp > 0 ; //@ a s s e r t tmp <= 1 0 0 ; r e t u r n tmp ; } It is important to notice that in this translation, the precondition is necessary as the result of addition may overflow and the destination is in 32 bits. We need to assure that the addition does not overflow. If it overflows then the Why program blocks.

(c) call instruction: call proc

We have i 1 i = f(): invokes a function in Why.

Because S, f () i ⇒ p S 1 , the rule of the function call implies

• P re S holds, where P re = P re annot • S, body ⇒ p-1 S 1 , where body = instr(P) i

• P ost S 1 holds Lemma 1 says that P re S holds.

By induction, because p -1 < m and there exists S 1 such that

• S, instr(f ) ⇒ S 1 and

From ( 1), (2) and with Lemma 1, we have P ost S 1 holds.

Consequently, the rule of call instruction in assembly applies: S, call f ⇒ S 1 .

Examples

Simple example

This example is presented in Section 6.4 and its assembly code is in Figure 6.4. The C source code is also shown in Figure 7.3. The Why program corresponding to the assembly code is presented in Figure 7.4. The assembly instruction is put in the Why comments under the form (*....*) and the translation of this instruction follows each instruction.

The function f has a precondition and no post-condition. Thus, in the parameter f_parameter the post-condition is true. The translation of this example is very simple. There are only two instructions we need to translate: movl and subl and two assertions. What we prove are these two assertions. They are proved by both Gappa, Alt-Ergo and CVC3 (See • For SSE/SSE2 mode, there is no difference between ucomiss/comiss and ucomisd/comisd.

• For x87 mode:

-We translate only floating-point operations, the integer instructions (fimul, fiadd, etc.) will not be considered here.

-We consider the case where the stack is empty at the entrance and the exit of a function (the top-of-stack pointer will point to the top register of the stack).

Translation of SSE/SSE2 instructions

We need to set a floating-point value to a register. To do that, we define the following parameter functions: Setting a single has two cases:

• Case 1: We do not need to check if the input value is overflow or not. We use it when transferring data from src to dest in the movss instructions. Pay attention that in this case, the value a is not rounded because we already know that it is a 32-bit floating-point number.

• Case 2: We have to check the input value. This parameter is used when we set a value of a computation (addition, subtraction, etc.) to a register.

We do similarly with double floating-point value. Division is a special case. We ensure that the divisor is not equal to 0.

Definition of programs supported

In this chapter, we continue to extend the model from the previous chapters. More precisely, the assembly programs may contain condition instructions, corresponding to the complex statements in C language: if then else, switch, for, do while, goto, etc. We add here some other constraints:

• It is possible that the compiler generates the calculated jump instruction for C statements such as switch. In this chapter, we do not handle this kind of jump instruction. We only consider the jump to a near relative address in which its operand is a label.

• Remind that in subsection 8.3.5, we try to translate "relative" stack registers into physical stack register statically. When compiling with x87 mode, one thing we want to make sure is that the value we calculate for top-of-stack register must be unique whatever is the path of the control-flow graph to reach the instruction. In all the programs that we handled in the previous chapter, this hypothesis was satisfied and we think that a compiler never generates the assembly code in which the top-of-stack value is not unique at each point.

Translation of comparison instructions

As conditional instructions depend on the status flags which are set by comparison instructions such as cmp, comisd, ucomisd, etc., we present here how to translate such instructions into Why. -Conditional move instruction i with condition c and label l: * explore_rec (l, prefix:: assume why_cond(c)::why_instr(i)) * explore_rec (n.succ, prefix:: assume not why_cond(c))

Finally, if n is the last node of the subgraph then produce_why_fun (prefix)

This algorithm fails only when there exists a cycle without any invariant. Its complexity is an exponential number of the size of the CFG. Such complexity is introduced by a following example: if c0 then t0 else e0; if c1 then t1 else e1; if c2 then t2 else e2; ... if ck then tk else ek;

The number of paths from n0 to s of this CFG is 2 k , so our algorithm will generate 2 k Why functions. However, this exponent is not a big problem. If we meet such case, it is possible to solve the problem manually by inserting invariants. With the example above, if we put manually invariants at each node nk then the number of Why functions will be 2k.

Theorem 9.3 (Soundness) Let p be the assembly code of a function, g its CFG, and E the set of Why programs generated by generateWhy(g). If the VCs for the Why programs in E are valid, then p satisfies its annotations.

Proof. From any state satisfying the precondition of p, a given execution corresponds to a path P in the graph g, from its initial node to a final node. If that path contains at least a cycle (with invariant I inside) then it can be split into subpaths p 1 , . . . , p n , such that p 1 goes from the initial node to a node with an invariant I; p 2 , . . . , p n-1 goes from a node containing I to a node with an invariant I; and p n goes from a node containing I to a final node. Each p i corresponds to an execution of one of the generated Why function, in a state satisfying its precondition (straightforward induction on n) and thus its satisfies the precondition of p.

Data transfer instructions

We define two parameters move_mem_to_reg64 and move_reg_to_mem32 which move the data (in 64 bits and in 32 bits) of a register to memory: The translation of the mov instructions are specified as follows: movl mem, reg i = move_cte32 mem int32 mem single mem exact reg movl reg, mem i = move_reg_to_mem32 !reg mem addr movq mem, reg i = move_cte64 mem int64 mem double mem exact reg movq reg, mem i = move_reg_to_mem64 !reg mem addr

The instruction lea loads effective address. It computes the effective address of the source operand and stores it in the destination operand. The source operand is a memory address (offset part), the destination operand is a general-purpose register. leaq mem, reg i = set_int64 mem int64 reg leal mem, reg i = set_int32 mem int32 reg

The instruction cltq extends the 32-bit register %eax to the 64-bit register %rax. In order to translate this instruction into Why, we declare a parameter:

and integer_of_int32(sel_int32(_rax)) = integer_of_int32(sel_int32(_rax@)) }

Remind that %eax = _rax. The translation fo cltq is

Notice that for all instructions, when dest is a register, we reuse the Why parameter of Chapter 7 and Chapter 8.

Arithmetic instructions

The instruction salq imm, dest multiplies the destination dest by 2, imm times. In other words, we multiply dest by 2 imm . We need the following parameter:

The translation of salq is as follows: salq imm, dest i = salq 2pimm dest where 2pimm = 2 imm . Notice that this is a bit-level instruction and we can not translate it into Why. Instead, the value 2pimm is calculated directly by our translator, not in generated Why program.

To store a 32-and 64-bit integer to memory, we need the following parameters: / * @ r e q u i r e s n > 0 && \ v a l i d _ r a n g e ( t , 0 , n -1); @ e n s u r e s \ f o r a l l i n t e g e r k ; 0<=k<n == > \ r e s u l t >= t [ k ] ; @ * / double max_array ( double t [ ] , i n t n ) { double m; i n t i = 0 ; goto L ;

Scalar Product

This example illustrates how we combine floating-point analysis with other features such as loops and arrays. It is presented in Section 5.5 and the same code source is in Figure 10.10.

In this example, the role of the preconditions \valid_range(x,0,n-1) and \valid_range(y,0,n-1) is important as these preconditions help to separate x[i] and p and also to separate y[i] and p. This works because:

The table below displays the value of B (found by Gappa tool) in function of NMAX and the architecture-dependent settings. The predicate eq defines that the value of the n th bit of v1 is equal to the n th bit of v2 for all values of n such that 0 ≤ n < size.

Theories

A Why 3 input is a file organized as a list of theories. A theory is a list of declarations. Declarations introduce new types, functions and predicates, state axioms, lemmas and goals. These declarations can be directly written in the theory or taken from existing theories.

Excerpt of the theory BitVector in Figure 11.6 illustrates how to define a theory in Why 3. The detail of this theory will be talked in the next section.

This theory has a list of declarations:

• the size of a bitvector,

• a type bv representing a bitvector,

• the function nth: returns the value of the bitvector bv at a given position,

• functions bvzero and bvone and the corresponding axioms Nth_zero and Nth_one: express that bvzero is a bitvector with all values being False and bvone is a bitvector with all value being True,

• a predicate eq.

use and clone

When we need to reuse a theory, we use the keyword use or clone. If we use a theory without duplicate it, we use the keyword use. On the contrary, when we clone a theory with the keyword clone, we create a local copy of every cloned declaration, and the newly created symbols, despite having the same names, are different from their originals. More precisely: If a theory T1 is used in another theory T2 then

• the symbols of T1 are shared . i ] * ) f u n c t i o n to_nat_sub bv i n t i n t : i n t ( * ( to_nat_sub b j i ) r e t u r n s t h e non-n e g a t i v e i n t e g e r whose b i n a r y r e p r i s b [ j . . i ] * ) . . . .

( * ( f r o m _ i n t n ) r e t u r n s t h e b i t v e c t o r r e p r e s e n t i n g t h e i n t e g e r n on s i z e b i t s . * ) f u n c t i o n f r o m _ i n t ( n : i n t ) : bv . . . . ( * f r o m _ i n t 2 c : i n t -> bv t a k e an i n t e g e r a s i n p u t and r e t u r n s a bv w i t h 2-complement * ) f u n c t i o n f r o m _ i n t 2 c ( n : i n t ) : bv . . . . end The declaration of the variable Const and the proofs related to Const are shown in Figure 11.14. We declare the function const that returns a 64-bit bitvector by concatenating j and j' (line 8). The function const_as_double returns the real value of bitvector const (line 10). We need to prove that Const = 2 52 + 2 31 . The corresponding lemma const_value is at line 39. All the lemmas from line 15 -37 help to prove const_value. All the lemmas for proving Const = 2 52 + 2 31 are proved automatically by Alt-Ergo, CVC3 and Z3 except the lemma exp_const which is proved using Coq (30 lines). There are 5 other lemmas that help to prove the lemmas above by automatic provers.

Proof of Var(x) = 2 52 + 2 31 The lemma lemma1 helps to prove the lemma lemma2. lemma2 proves that the fraction value of var(x) is 2 31 + x. The lemmas lemma2, lemma3, lemma4 helps to find the double value of var(x).

The lemma lemma1 is the most difficult one to prove. The proofs of lemma2, lemma3, lemma4 are similar to the proofs of fraction, exponent and sign in the case of const. The Why code of lemma1 and others lemmas needed to prove lemma1 is shown in Figure 11. [START_REF] Burdy | Java bytecode specification and verification[END_REF].

These lemmas are proved by the combination of automatic provers and interactive prover Coq. They are shown in the following table:

Abstract

On some recently developed architectures, a numerical program may give different answers depending on the execution hardware and the compilation. These discrepancies of the results come from the fact that each floating-point computation is calculated with different precisions. The goal of this thesis is to formally prove properties about numerical programs while taking the architecture and the compiler into account. In order to do that, we propose two different approaches. The first approach is to prove properties of floating-point programs that are true for multiple architectures and compilers. This approach states the rounding error of each floatingpoint computation whatever the environment and the compiler choices. It is implemented in the Frama-C platform for static analysis of C code. The second approach is to prove behavioral properties of numerical programs by analyzing their compiled assembly code. We focus on the issues and traps that may arise on floating-point computations. Direct analysis of the assembly code allows us to take into account architecture-or compiler-dependent features such as the possible use of extended precision registers. It is implemented above the Why platform for deductive verification.

Keywords: Floating-point arithmetic, Numerical programs, Static analysis, Compile-time optimizations, the Why platform, the Frama-C platform.

Résumé

Sur des architectures récentes, un programme numérique peut donner des réponses différentes en fonction du hardware et du compilateur. Ces incohérences des résultats viennent du fait que chaque calcul en virgule flottante est effectué avec des précisions différentes. Le but de cette thèse est de prouver formellement des propriétés des programmes opérant sur des nombres flottants en prenant en compte l'architecture et le compilateur. Pour le faire, nous avons proposé deux approches différentes. La première approche est de prouver des propriétés des programmes en virgule flottante qui sont vraies sur plusieurs architectures et compilateurs. Cette approche ne considère que les erreurs d'arrondi qui doivent être validées quels que soient l'environnement matériel et le choix du compilateur. Elle est implantée dans la plate-forme Frama-C pour l'analyse statique de code C. La deuxième approche consiste à prouver des propriétés des programmes en analysant leur code assembleur. Nous nous concentrons sur des problèmes et des pièges qui apparaissent sur des calculs en virgule flottante. L'analyse directe du code assembleur nous permet de considérer des caratéristiques dépendant de l'architecture ou du compilateur telle que l'utilisation des registres en précision étendue. Cette approche est implantée comme une sur-couche de la plate-forme Why pour la vérification déductive.

Mots clés: Arithmétique en virgule flottante, Programmes numériques, Analyse statique, Optimisations à la compilation, Plate-forme Why, Plate-forme Frama-C.