Correct Approximation of IEEE 754 Floating-Point Arithmetic for Program Verification

Verification of programs using floating-point arithmetic is challenging on several accounts. One of the difficulties of reasoning about such programs is due to the peculiarities of floating-point arithmetic: rounding errors, infinities, non-numeric objects (NaNs), signed zeroes, denormal numbers, different rounding modes, etc. One possibility to reason about floating-point arithmetic is to model a program computation path by means of a set of ternary constraints of the form z = x op y and use constraint propagation techniques to infer new information on the variables' possible values. In this setting, we define and prove the correctness of algorithms to precisely bound the value of one of the variables x, y or z, starting from the bounds known for the other two. We do this for each of the operations and for each rounding mode defined by the IEEE 754 binary floating-point standard, even in the case the rounding mode in effect is only partially known. This is the first time that such so-called filtering algorithms are defined and their correctness is formally proved. This is an important slab for paving the way to formal verification of programs that use floating-point arithmetics.


Introduction
Programs using floating-point numbers are notoriously difficult to reason about [33]. Many factors complicate the task: 1. compilers may transform the code in a way that does not preserve the semantics of floating-point computations; 2. floating-point formats are an implementation-defined aspect of most programming languages; 3. there are different, incompatible implementations of the operations for the same floating-point format; 4. mathematical libraries often come with little or no guarantee about what is actually computed; 5. programmers have a hard time predicting and avoiding phenomena caused by the limited range and precision of floating-point numbers (overflow, absorption, cancellation, underflow, etc.); moreover, devices that modern floating-point formats possess in order to support better handling of such phenomena (infinities, signed zeroes, denormal numbers, non-numeric objects a.k.a. NaNs) come with their share of issues; 6. rounding is a source of confusion in itself; moreover, there are several possible rounding modes and programs can change the rounding mode any time.
As a result of these difficulties, the verification of floating-point programs in industry relies, almost exclusively, on informal methods, mainly testing, or on the evaluation of the numerical accuracy of computations, which only allows to determine conservative (but often too loose) bounds on the propagated error [19]. The satisfactory formal treatment of programs engaging in floating-point computations requires an equally satisfactory solution to the difficulties summarized in the above enumeration. Progress has been made, but more remains to be done. Let us review each point: 1. Some compilers provide options to refrain from rearranging floating-point computations. When these are not available or cannot be used, the only possibility is to verify the generated machine code or some intermediate code whose semantics is guaranteed to be preserved by the compiler back-end. 2. Even though the used floating-point formats are implementation-defined aspects of, say, C and C++ 1 the wide adoption of the IEEE 754 standard for binary floating-point arithmetic [24] has improved things considerably. 3. The IEEE 754 standard does provide some strong guarantees, e.g., that the results of individual additions, subtractions, multiplications, divisions and square roots are correctly rounded, that is, it is as if the results were computed in the reals and then rounded as per the rounding mode in effect. But it does not provide guarantees on the results of other operations and on other aspects, such as, e.g., when the underflow exception is signaled [17]. 2 4. A pragmatic, yet effective approach to support formal reasoning on commonly used implementation of mathematical functions has been recently proposed in [6]. The proposed techniques exploit the fact that the floating-point implementation of mathematical functions preserve, not completely but to a great extent, the piecewise monotonicity nature of the approximated functions over the reals. 5. A static analysis for detecting floating-point exceptions based on abstract interpretation has been presented in [32]. A few attempts at this task have been made using other techniques [7,38] but, as we argue in Sections 1.4 and 5, they present precision and soundness issues. 6. Most verification approaches in the literature assume the round-to-nearest rounding mode [10], or over-approximate by always considering worst-case rounding modes [32]. Analyses based on SMT solvers [13] can treat each rounding mode precisely, but only if the rounding mode in use is known exactly. As we show in Section 5, some SMT solvers also suffer from soundness issues.
The contribution of this paper is in areas 5 and 6. In particular, concerning point 5, by defining and formally proving the correctness of constraint propagation algorithms for IEEE 754 arithmetic constraints, we enable the use of formal methods for a broad range of programs. Such methods, i.e., abstract interpretation and symbolic model checking, allow for proving that a number of generally unwanted phenomena (e.g., generation of NaNs and infinities, absorption, cancellation, instability, etc.) do not happen or, in case they do happen, allow the generation of a test vector to reproduce the issue. Regarding point 6, handling of all IEEE 754 rounding modes, and being resilient to uncertainty about the rounding mode in effect, is another original contribution of this paper.
While the round-to-nearest rounding mode is, by far, the most frequently used one, it must be taken into account that: -the possibility of programmatically changing the rounding mode is granted by IEEE 754 and is offered by most of its implementations (e.g., in the C programming language, via the fesetround() standard function); -such possibility is exploited by interval libraries and by numerical calculus algorithms (see, e.g., [35,36]); -setting the rounding mode to something different from round-to-nearest can be done by third parties in a way that was not anticipated by programmers: this may cause unwanted non-determinism in video games [20] and there is nothing preventing the abuse of this feature for more malicious ends, denial-of-service being only the least dangerous in the range of possibilities. Leaving malware aside, there are graphic and printer drivers and sound libraries that are known to change the rounding mode and may fail to set it back [37].
As a possible way of tackling the difficulties described until now, and enabling sound formal verification of floating-point computations, this paper introduces new algorithms for the propagation of arithmetic constraints over floating-point numbers. These algorithms are called filtering algorithms as their purpose is to prune the domains of possible variable values by filtering out those values that cannot be part of the solution of a system of constraints. Algorithms of this kind must be employed in constraint solvers that are required in several different areas, such as automated test-case generation, exception detection or the detection of subnormal computations. In this paper we propose fully detailed, provably correct filtering algorithms for floating-point constraints. Such algorithms handle all values, including symbolic values (NaNs, infinities and signed zeros), and rounding modes defined by IEEE 754. Note that filtering techniques used in solvers over the reals do not preserve all solutions of constraints over floating-point numbers [30,31], and therefore they cannot be used to prune floating-point variable domains reliably. This leads to the need of filtering algorithms such as those we hereby introduce.
The choice of the IEEE 754 Standard for floating-point numbers as the target representation for our algorithms is due to their ubiquity in modern computing platforms. Indeed, although some programming languages leave the floating-point format as an implementation-defined aspect, all widely-used hardware platforms -e.g., x86 [25] and ARM [3]-only implement the IEEE 754 Standard, while older formats are considered legacy.
Before defining our filtering algorithms in a detailed and formal way, we provide a more comprehensive context on the propagation of floating-point constraints and its practical applications (Sections 1.1 and 1.2), and justify their use in formal program analysis and verification (Section 1.3). We also give a more in-depth view of related work in Section 1.4, and clarify our contribution in Section 1.5.

From programs to floating-point constraints
Independently from the application, program analysis starts with parsing, the generation of an abstract syntax tree and the generation of various kinds of intermediate program representations. An important intermediate representation is called three-address code (TAC). In this representation, complex arithmetic expressions and assignments are decomposed into sequences of assignment instructions of the form result := operand 1 operator operand 2 .
A further refinement is the computation of the static single assignment form (SSA) [2] whereby, labeling each assigned variable with a fresh name, assignments can be considered as if they were equality constraints. For example, the TAC form of the floating-point assignment z := z * z + z is t := z * z; z := t + z, which in an SSA form becomes t 1 := z 1 * z 1 ; z 2 := t 1 + z 1 . These, in turn, can be regarded as the conjunction of the constraints t 1 = z 1 z 1 and z 2 = t 1 z 1 , where by and we denote the multiplication and addition operations on floating-point numbers, respectively. The Boolean comparison expressions that appear in the guards of if statements and loops can be translated into constraints similarly. This way, a C/C++ program translated into an SSA-based intermediate representation can be represented as a set of constraints on its variables. In particular, a constraint set arises form each execution path in the program. For this reason, this approach to program modeling can be viewed as symbolic execution [15,26]. Constraints can be added or removed from such a set in order to obtain a constraint system that describes a particular behavior of the program (e.g., the execution of a certain instruction, the occurrence of an overflow in a computation, etc.). Once such a constraint system has been solved, the variable domains only contain values that cause the desired behavior. If one of the domains is empty, then that behavior can be ruled out. For more details on the symbolic execution of floating-point computations, we refer the reader to [4,10].

Constraint propagation
Once constraints have been generated, they are amenable to constraint propagation: under this name goes any technique that entails considering a subset of the constraints at a time, explicitly removing elements from the set of values that are candidate to be assigned to the constrained variables. The values that can be removed are those that cannot possibly participate in a solution for the selected set of constraints. For instance, if a set of floating-point constraints contains the constraint x x = x, then any value outside the set {NaN, +0, 1, +∞} can be removed from further consideration. The degree up to which this removal can actually take place depends on the data-structure used to record the possible values for x, intervals and multi-intervals being typical choices for numerical constraints. For the example above, if intervals are used, the removal can only be partial (negative floating-point numbers are removed from the domain of x). With multi-intervals more precision is possible, but any approach based on multi-intervals must take measures to avoid combinatorial explosion.
In this paper, we only focus on interval-based constraint propagation: the algorithms we present for intervals can be rather easily generalized to the case of multi-intervals. We make the further assumption that the floating-point formats available to the analyzed program are also available to the analyzer: this is indeed quite common due to the wide adoption of the IEEE 754 formats.
Interval-based floating-point constraint propagation consists of iteratively narrowing the intervals associated to each variable: this process is called filtering. A projection is a function that, given a constraint and the intervals associated to two of the variables occurring in it, computes a possibly refined interval for the third variable (the projection is said to be over the third variable). Taking z 2 = t 1 z 1 as an example, the projection over z 2 is called direct projection (it goes in the same sense of the TAC assignment it comes from), while the projections over t 1 and z 1 are called indirect projections.

Applications of constraint propagation to program analysis
When integrated in a complete program verification framework, the constraint propagation techniques presented in this paper enable activities such as abstract interpretation, automatic test-input generation and symbolic model checking. In particular, symbolic model checking means exhaustively proving that a certain property, called specification, is satisfied by the system in exam, which in this case is a computer program. A model checker can either prove that the given specification is satisfied, or provide a useful counterexample whenever it is not.
For programs involving floating-point computations, some of the most significant properties that can be checked consist of ruling out certain undesired exceptional behaviors such as overflows, underflows and the generation of NaNs, and numerical pitfalls such as absorption and cancellation. In more detail, we call a numeric-to-NaN transition a floating-point arithmetic computation that returns a NaN despite its operands being non-NaN. We call a finite-to-infinite transition the event of a floating-point operation returning an infinity when executed on finite operands, which occurs if the operation overflows. An underflow occurs when the output of a computation is too small to be represented in the machine floatingpoint format without a significant loss in accuracy. Specifically, we divide underflows into three categories, depending on their severity: Gradual underflow: an operation performed on normalized numbers results in a subnormal number. In other words, a subnormal has been generated out of normalized numbers: enabling gradual underflow is indeed the very reason for the existence of subnormals in IEEE 754. However, as subnormals come with their share of problems, generating them is better avoided. Hard underflow: an operation performed on normalized numbers results in a zero, whereas the result computed on the reals is nonzero. This is called hard because the relative error is 100%, gradual overflow does not help (the output is zero, not a subnormal), and, as neither input is a subnormal, this operation may constitute a problem per se. Soft underflow: an operation with at least one subnormal operand results in a zero, whereas the result computed on the reals is nonzero. The relative error is still 100% but, as one of the operands is a subnormal, this operation may not be the root cause of the problem.
Absorption occurs when the result of an arithmetic operation is equal to one of the operands, even if the other one is not the neutral element of that operation. For example, absorption occurs when summing a number with another one that has a relatively very small exponent. If the precision of the floating-point format in use is not enough to represent them, the additional digits that would appear in the mantissa of the result are rounded out. In this section, we show how symbolic model checking can be used to either rule out or pinpoint the presence of these run-time anomalies in a software program by means of a simple but meaningful practical example. Floating-point constraint propagation has been fully implemented with the techniques presented in this paper in the commercial tool ECLAIR, 3 developed and commercialized by BUGSENG. ECLAIR is a generic platform for the formal verification of C/C++ and Java source code, as well as Java bytecode. The filtering algorithms described in the present paper are used in the C/C++ modules of ECLAIR that are responsible for semantic analysis based on abstract interpretation [16], automatic generation of test-cases, and symbolic model checking. The latter two are based on symbolic execution and constraint satisfaction problems [22,23], whose solution is based on multiinterval refinement and is driven by labeling and backtracking search. Indeed, the choice of ECLAIR as our target verification platform is mainly due to its use of constraint propagation for solving constraints generated by symbolic execution, which makes it easier to integrate the algorithms presented in this paper. However, such techniques are general, and could be used to solve the constraints generated by any symbolic execution engine.
Constraints arising from the use of mathematical functions provided by C/C++ standard libraries are also supported. Unfortunately, most implementations of such libraries are not correctly rounded, which makes the realization of filtering algorithms for them rather challenging. In ECLAIR, propagation for such constraints is performed by exploiting the piecewise monotonicity properties of those functions, which are partially retained by all implementations we know of [6].
To demonstrate the capabilities of the techniques presented in this paper, we applied them to the C code excerpt of Fig. 1. It is part of the implementation of the Bessel functions in the GNU Scientific Library, 4 a widely adopted library for numerical computations. In particular, it computes the scaled regular modified cylindrical Bessel function of first order, exp(−|x|)I 1 (x), where x is a purely imaginary argument. The function stores the Fig. 1 Function extracted from the GNU Scientific Library (GSL), version 2.5. The possible numerical exceptions detected by ECLAIR are marked by the raised letters next to the operators causing them. h, s and g stand for hard, soft and gradual underflow, respectively; a for absorption; i for finite-to-infinity; n for numeric-to-NaN computed result in the val field of the data structure result, together with an estimate of the absolute error (result->err). Additionally, the function returns an int status code, which reports to the user the occurrence of certain exceptional conditions, such as overflows and underflows. In particular, this function only reports an underflow when the argument is smaller than a constant. We analyzed this program fragment with ECLAIR's symbolic model checking engine, setting it up to detect overflow (finite-to-infinite transitions), underflow and absorption events, and NaN generation (numeric-to-NaN transitions). Thus, we found out the underflow guarded against by the if statement of line 12 is by far not the only numerical anomaly affecting this function. In total, we found a numeric-to-NaN transition, two possible finite-to-infinite transitions, two hard underflows, 5 gradual underflows and 6 soft underflows. The code locations in which they occur are all reported in Fig. 1.
For each one of these events, ECLAIR yields an input value causing it. Also, it optionally produces an instrumented version of the original code, and runs it on every input it reports, checking whether it actually triggers the expected behavior or not. Hence, the produced input values are validated automatically. For example, the hard underflow of line 17 is triggered by the input x = -0x1.8p-1021 ≈ −6.6752 × 10 −308 . If the function is executed with x = -0x1p+1023 ≈ −8.9885 × 10 307 , the multiplication of line 29 yields a negative infinity. Since ax = |x|, we know x = 0x1p+1023 would also cause the overflow. The same value of x causes an overflow in line 30 as well. The division in the same line produces a NaN if the function is executed with x = −∞.
The context in which the events we found occur determines whether they could cause significant issues. For example, even in the event of absorption, the output of the overall computation could be correctly rounded. Whether or not this is acceptable must be assessed depending on the application. Indeed, the capability of ECLAIR of detecting absorption can be a valuable tool to decide if a floating-point format with a higher precision is needed. Nevertheless, some of such events are certainly problematic. The structure of the function suggests that no underflow should occur if control flow reaches past the if guard of line 12. On the contrary, several underflows may occur afterwards, some of which are even hard. Moreover, the generation of infinities or NaNs should certainly either be avoided, or signaled by returning a suitable error code (and not GSL SUCCESS). The input values reported by ECLAIR could be helpful for the developer in fixing the problems detected in the function of Fig. 1. Furthermore, the algorithms presented in this paper are provably correct. For this reason, it is possible to state that this code excerpt presents no other issues besides those we reported above. Notice, however, that due to the way the standard C mathematical library functions are treated, the results above only hold with respect to the implementation of the exp function in use. In particular, the machine we used for the analysis is equipped with the x86 64 version of EGLIBC 2.19, running on Ubuntu 14.04.1.

Filtering algorithms
In [30] C. Michel proposed a framework for filtering constraints over floating-point numbers. He considered monotonic functions over one argument and devised exact direct and correct indirect projections for each possible rounding mode. Extending this approach to binary arithmetic operators is not an easy task. In [10], the authors extended the approach of [30] by proposing filtering algorithms for the four basic binary arithmetic operators when only the round-to-nearest tails-to-even rounding mode is available. They also provided tables for indirect function projections when zeros and infinities are considered with this rounding mode. In our approach, we generalize the initial work of [10] by providing extended interval reasoning. The algorithms and tables we present in this paper consider all rounding modes, and contain all details and special cases, allowing the interested reader to write an implementation of interval-based filtering code.
Recently, [21] presented optimal inverse projections for addition under the round-tonearest rounding mode. The proposed algorithms combine classical filtering based on the properties of addition with filtering based on the properties of subtraction constraints on floating-points as introduced by Marre and Michel [29]. The authors are able to prove the optimality of the lower bounds computed by their algorithms. However, [21] only covers addition in the round-to-nearest rounding mode, leaving other arithmetic operations (subtraction, multiplication and division) and rounding modes to future work. Special values (infinities and NaNs) are also not handled. Conversely, this paper presents filtering algorithms covering all such cases.
It is worth noting that the filtering algorithms on intervals presented in [29] have been corrected for addition/subtraction constraints and extended to multiplication and division under the round-to-nearest rounding mode by some of these authors (see [4,5]). In this paper we discuss the cases in which the filtering algorithms in [4,5,29] should be used in combination with our filters for arithmetic constraints. However, the main aim of this paper is to provide an exhaustive and provably correct treatment of filtering algorithms supporting all special cases for all arithmetic constraints under all rounding modes.

SMT solvers
Satisfiability Modulo Theories (SMT) is the problem of deciding satisfiability of first-order logic formulas containing terms from different, pre-defined theories. Examples of such theories are integer or real arithmetic, bit-vectors, arrays and uninterpreted functions. Recently, SMT solvers have been widely employed as backends for different software verification techniques, such as model checking and symbolic execution [9]. The need for verifying floating-point programs lead to the introduction of a floating-point theory [13] in SMT-LIB, a library defining a common input language for SMT solvers. Since then, the theory has been implemented in different ways into several solvers. CVC4 [8,12], MathSAT [14] and Z3 [18] use bit-blasting, i.e., they convert floating-point constraints to bit-vector formulae, which are then solved as Boolean SAT problems. Some tools, instead, use methods based on interval reasoning. MathSAT also supports Abstract Conflict Driven Learning (ACDL) for solving floating-point constraints based on interval domains [11]. Colibri [28] uses constraint programming techniques, with filtering algorithms such as those in [5,10] and those presented in this paper. However, [28] does not report such filters in detail, nor proves their correctness. This leads to serious soundness issues, as we shall see in Section 5.2. An experimental comparison of such tools can be found in [12].
Note that SMT-LIB, the input language of all such tools, only allows to specify one single rounding mode for each floating-point operation. Thus, the only way of dealing with uncertainty of the rounding mode in use is to solve the same constraint system with all possible rounding mode combinations, which is quite unpractical. Our filtering algorithms are instead capable of working with a set of possible rounding modes, and retain soundness by always choosing the worst-case one.

Floating-point program verification
Several program analyses for automatic detection of floating-point exceptions were proposed in the literature.
Relational abstract domains for the analysis of floating-point computations through abstract interpretation have been presented in [32] and implemented in the tool Astrée. 5 Such domains, however, over-approximate rounding operations by always assuming worst cases, namely rounding toward plus and minus infinity, which may cause precision issues (i.e., false positives) if only round-to-nearest is used. Also, [32] does not offer a treatment of symbolic values (e.g., infinities) as exhaustive as the one we offer in this paper.
In [7] the authors proposed a symbolic execution system for detecting floatingpoint exceptions. It is based on the following steps: each numerical program is transformed to directly check each exception-triggering condition, the transformed program is symbolically-executed in real arithmetic to find a (real) candidate input that triggers the exception, the real candidate is converted into a floating-point number, which is finally tested against the original program. Since approximating floating-point arithmetic with real arithmetic does not preserve the feasibility of execution paths and outputs in any sense, they cannot guarantee that once a real candidate has been selected, a floating-point number raising the same exception can be found. Even more importantly, even if the transformed program over the reals is exception-free, the original program using floating-point arithmetic may not be actually exception-free.
Symbolic execution is also the basis of the analysis proposed in [38], that aims at detecting floating-point exceptions by combining it with value range analysis. The value range of each variable is updated with the appropriate path conditions by leveraging interval constraint-propagation techniques. Since the projections used in that paper have not been proved to yield correct approximations, it can be the case that the obtained value ranges do not contain all possible floating-point values for each variable. Indeed, valid values may be removed from value ranges, which leads to false negatives. In Section 6, the tool for floatingpoint exception detection presented in [38] is compared with the same analysis based on our propagation algorithms. As expected, no false positives were detected among the results of our analysis.

Contribution
This paper improves the state of the art in several directions: 1. all rounding modes are treated and there is no assumption that the rounding mode in effect is known and unchangeable (increased generality); 2. utilization, to a large extent, of machine floating-point arithmetic in the analyzer with few rounding mode changes (increased performance); 3. accurate treatment of round half to even -the default rounding mode of IEEE 754-(increased precision); 4. explicit and complete treatment of intervals containing symbolic values (i.e., infinities and signed zeros); 5. application of floating-point constraint propagation techniques to enable detection of program anomalies such as overflows, underflows, absorption, generation of NaNs.

Plan of the paper
The rest of the paper is structured as follows: Section 2 recalls the required notions and introduces the notation used throughout the paper; Section 3 presents some results on the treatment of uncertainty on the rounding mode in effect and on the quantification of the rounding errors committed in floating-point arithmetic operations; Section 4 contains the complete treatment of addition and division constraints on intervals, by showing detailed special values tables and the refinement algorithms; Section 5 reports the results of experiments aimed at evaluating the soundness of existing tools. Section 6 concludes the main part of the paper. Appendix A contains the complete treatment of subtraction and multiplication constraints. The proofs of results not reported in the main text of the paper can be found in Appendix B.

Preliminaries
We will denote by R + and R − the sets of strictly positive and strictly negative real numbers, respectively. The set of affinely extended reals, R ∪ {−∞, +∞}, is denoted by R.
Definition 2 (IEEE 754 binary floating-point numbers) A set of IEEE 754 binary floating-point numbers [24] is uniquely identified by: p ∈ N, the number of significant digits (precision); e max ∈ N, the maximum exponent, the minimum exponent being e min def = 1 − e max . The set of binary floating-point numbers F(p, e max , e min ) includes: -all signed zero and non-zero numbers of the form (−1) s · 2 e · m, where s is the sign bit; -the exponent e is any integer such that e min ≤ e ≤ e max ; -the mantissa m, with 0 ≤ m < 2, is a number represented by a string of p binary digits with a "binary point" after the first digit: -the infinities +∞ and −∞; the NaNs: qNaN (quiet NaN) and sNaN (signaling NaN).
Numbers such that d 0 = 1 are called normal. The smallest positive normal floating-point number is f nor min def = 2 e min and the largest is f max The non-zero floatingpoint numbers such that d 0 = 0 are called subnormal: their absolute value is less than 2 e min , and they always have fewer than p significant digits. Every finite floating-point number is an integral multiple of the smallest subnormal magnitude f min def = 2 e min +1−p . Note that the signed zeroes +0 and −0 are distinct floating-point numbers. For a non-zero number x, we will write even(x) (resp., odd(x)) to signify that the least significant digit of x's mantissa, d p−1 , is 0 (resp., 1).
In the sequel we will only be concerned with IEEE 754 binary floating-point numbers and we will write simply F for F(p, e max , e min ) when there is no risk of confusion.
Definition 3 (Floating-point symbolic order) Let F be any IEEE 754 floating-point format. The relation ≺⊆ F × F is such that, for each x, y ∈ F, x ≺ y if and only if both x and y are not NaNs and either: x = −∞ and y = −∞, or x = +∞ and y = +∞, or x = −0 and y ∈ {+0} ∪ R + , or x ∈ R − ∪ {−0} and y = +0, or x, y ∈ R and x < y. The partial order ⊆ F × F is such that, for each x, y ∈ F, x y if and only if both x and y are not NaNs and either x ≺ y or x = y.
Note that F without the NaNs is linearly ordered with respect to '≺'. For x ∈ F that is not a NaN, we will often abuse the notation by interchangeably using the floating-point number or the extended real number it represents. The floats −0 and +0 both correspond to the real number 0. Thus, when we write, e.g., x < y we mean that x is numerically less than y (for example, we have −0 ≺ +0 though −0 ≮ +0, but note that x y implies x ≤ y). Numerical equivalence will be denoted by '≡' so that

Definition 4 (Floating-point predecessors and successors) The partial function
otherwise.
The partial function F F is defined by reversing the ordering, so that, for each x ∈ F, pred(x) = −succ(−x) whenever succ(x) is defined.
The rounding functions are defined as follows. Note that they are not defined for 0: the IEEE 754 standard, in fact, for operations whose exact result is 0, bases the choice between +0 and −0 on the operation itself and on the sign of the arguments [24, Section 6.3].
Note that, when the result of an operation has magnitude lower than f nor min , it is rounded to a subnormal number, by adjusting it to the form (−1) s · 2 e min · m, and truncating its mantissa m, which now starts with at least one 0, to the first p digits. This phenomenon is called gradual underflow, and while it is preferred to hard underflow, which truncates a number to 0, it still may cause precision issues due to the reduced number of significant digits of subnormal numbers.
The rounding modes ↓ and ↑ are the most "extreme", while n and 0 are always contained between them. We formalize this observation as follows: Moreover, In this paper, we use intervals of floating-point numbers in F that are not NaNs.
Definition 6 (Floating-point intervals) Let F be any IEEE 754 floating-point format. The set I F of floating-point intervals with boundaries in F is given by The set I F is a bounded meetsemilattice with least element ∅, greatest element [−∞, +∞], and the meet operation, which is induced by set-intersection, will be simply denoted by ∩.
Floating-point intervals with boundaries in F allow to capture the extended numbers in F: NaNs should be tracked separately.

Rounding modes and rounding errors
The IEEE 754 standard for floating-point arithmetic introduces different rounding operators, among which the user can choose on compliant platforms. The rounding mode in use affects the results of the floating-point computations performed, and it must be therefore taken into account during constraint propagation. In this section, we present some abstractions aimed at facilitating the treatment of rounding modes in our constraint projection algorithms.

Dealing with uncertainty on the rounding mode in effect
Even if programs that change the rounding mode in effect are quite rare, whenever this happens, the rounding mode in effect at each program point cannot be known precisely. So, for a completely general treatment of the problem, such as the one we are proposing, our choice is to consider a set of possible rounding modes. To this aim, in this section we define two auxilliary functions that, given a set of rounding modes possibly in effect, select a worst-case rounding mode that ensures soundness of interval propagation. Soundness is guaranteed even if the rounding mode used in the actual computation differs from the one selected, as far as the former is contained in the set. Of course, if a program never changes the rounding mode, the set of possible rounding modes boils down to be a singleton.
The functions presented in the first definition select the rounding modes that can be used to compute the lower (function r ) and upper (function r u ) bounds of an operation in case of direct projections.
Definition 7 (Rounding mode selectors for direct projections) Let F be any IEEE 754 floating-point format and S ⊆ R be a set of rounding modes. Let also y, z ∈ F and be such that either . Then The following functions select the rounding modes that will be used for the lower (functionsr r andr ) and upper (functionsr r u andr u ) bounds of an operation when computing inverse projections. Note that there are different functions depending on which one of the two operands is being projected:r r andr r u for the right one,r andr u for the left one.
Definition 8 (Rounding mode selectors for inverse projections) Let F be any IEEE 754 floating-point format and S ⊆ R be a set of rounding modes. Let also a, b ∈ F and First, we define Secondly, we define the following selectors: The usefulness in interval propagation of the functions presented above will be clearer after considering Proposition 2. Moreover, it is worth noting that, if the set of possible rounding modes is composed by a unique rounding mode, then all the previously defined functions return such rounding mode itself. In that case, the claims of Proposition 2 trivially hold.
Proposition 2 Let F, S, y, z and ' ' be as in Definition 7. Let also and . Then, for each r ∈ S Moreover, there exist r , r ∈ S such that (10) Now, consider with x, z ∈ F and r ∈ S. Let and , according to Definition 8. Moreover, letŷ be the minimum y ∈ F such that , and letỹ be the maximum y ∈ F such that . Then, the following inequalities hold:ŷ y ỹ .
The same result holds if , with and .
Proof Here we only prove the claims for direct projections (namely, (9) and (10) In order to prove inequality (9), it is now sufficient to consider all possible sets S ⊆ R and use the relations above. For claim (10), observe that for any combination of rounding modes in S except for one case: that is when y • z > 0, and 0 ∈ S but ↓ / ∈ S. In this case, however, by Definition 5, . Similarly, except for the case when y • z ≤ 0, 0 ∈ S but ↑ / ∈ S. First, assume that y • z < 0: in this case, by Definition 5, . For the remaining case, that is y • z = 0, we observe that for multiplication and division the result is the same for all rounding modes [24, Section 6.3], while for addition or subtraction we have .
Thanks to Proposition 2 we need not be concerned with sets of rounding modes, as any such set S ⊆ R can always be mapped to a pair of "worst-case rounding modes" which, in addition, are never round-to-zero. Therefore, projection functions can act as if the only possible rounding mode in effect was the one returned by the selection functions, greatly simplifying their logic. For example, consider the constraint , meaning "x is obtained as the result of for some r ∈ S." Of course, implies and , which, by Proposition 2, imply and , where and . The results obtained by projection functions that only consider r and r u are consequently valid for any r ∈ S.

Rounding errors
For the precise treatment of all rounding modes it is useful to introduce a notation that expresses, for each floating-point number x, the maximum error that has been committed by approximating with x a real number under the different rounding modes (as shown in the previous section, we need not be concerned with round-to-zero).
Definition 9 (Rounding Error Functions) The partial functions ∇ ↑ : F R, ∇ ↓ : F R, ∇ n− 2 : F R and ∇ n+ 2 : F R are defined as follows, for each x ∈ F that is not a NaN: An interesting observation is that the values of the functions introduced in Definition 9 are always representable in F and thus their computation does not require extra-precision, something that, as we shall see, is exploited in the implementation. This is the reason why, for round-to-nearest, ∇ n− 2 and ∇ n+ 2 have been defined as twice the approximation error bounds: the absolute value of the bounds themselves, being f min /2, is not representable in F for each x ∈ F such that |x| ≤ f nor min . When the round-to-nearest rounding mode is in effect, Proposition 3 relates the bounds of a floating-point interval [x , x u ] with those of the corresponding interval of R it represents.
Proof (sketch) To prove (15), we separately consider the two cases defined by (13).
By monotonicity of 'pred', the minimum value of (x + pred(x))/2 occurs when x = succ(−f max ), and which proves (15) in this case. If, instead, x > −f max , applying monotonicity of 'pred' suffices. The full proof is in Appendix B.2, together with the one of (16), which is symmetric.

Real approximations of floating-point constraints
In this section we show how inequalities of the form and , with r ∈ {↓, ↑, n} can be reflected on the reals. Indeed, it is possible to algebraically manipulate constraints on the reals so as to numerically bound the values of floating-point quantities. The results of this and of the next section will be useful in designing inverse projections.

Proposition 4
The following implications hold, for each x, y, z ∈ F such that all the involved expressions do not evaluate to NaN, for each floating-point operation and the corresponding extended real operation • ∈ {+, −, ·, /}, where the entailed inequalities are to be interpreted over R: The proof of Proposition 4 is carried out by applying the inequalities of Proposition 1 to each rounded operation, resulting in a quite long case analysis. It can be found in Appendix B.2.

Floating-point approximations of constraints on the reals
In this section, we show how possibly complex constraints involving floating-point operations can be approximated directly using floating-point computations, without necessarily using infinite-precision arithmetic.
Without being too formal, let us consider the domain E F of abstract syntax trees with leafs labelled by constants in F and internal nodes labeled with a symbol in {+, −, ·, /} denoting an operation on the reals. While developing propagation algorithms, it is often necessary to deal with inequalities between real numbers and expressions described by such syntax trees. In order to successfully approximate them using the available floating-point arithmetic, we need two functions: [[·]] ↓ : E F → F and [[·]] ↑ : E F → F. These functions provide an abstraction of evaluation algorithms that: (a) respect the indicated approximation direction; and (b) are as precise as practical. Point (a) can always be achieved by substituting the real operations with the corresponding floating-point operations rounded in the right direction. For point (b), maximum precision can trivially be achieved whenever the expression involves only one operation; generally speaking, the possibility of efficiently computing a maximally precise (i.e., correctly rounded) result depends on the form of the expression (see, e.g., [27]).

Definition 10 (Evaluation functions) The two partial functions
The implications of Proposition 5 can be derived from Definition 10 and Proposition 1. Their proof is postponed to Appendix B.2.

Propagation for simple arithmetic constraints
In this section we present our propagation procedure for the solution of floating-point constraints obtained from the analysis of programs engaging in IEEE 754 computations.
The general propagation algorithm, which we already introduced in Section 1.2, consists in an iterative procedure that applies the direct and inverse filtering algorithms associated with each constraint, narrowing down the intervals associated with each variable. The process stops when fixed point is reached, i.e., when a further application of any filtering algorithm does not change the domain of any variable.

Propagation algorithms: definitions
Constraint propagation is a process that prunes the domains of program variables by deleting values that do not satisfy any of the constraints involving those variables. In this section, we will state these ideas more formally.
Let and S ⊆ R. Consider a constraint with x ∈ X = [x , x u ], y ∈ Y = [y , y u ] and z ∈ Z = [z , z u ]. Direct propagation aims at inferring a narrower interval for variable x, by considering the domains of y and z. It amounts to computing a possibly refined interval for x, X = [x , x u ] ⊆ X, such that (31) Property (31) is known as the direct propagation correctness property. Of course it is always possible to take X = X, but the objective of constraint propagation is to compute a "small", possibly the smallest X enjoying (31), compatibly with the available information. The smallest X that satisfies (31) is called optimal and is such that (32) Property (32) is called the direct propagation optimality property. Inverse propagation, on the other hand, uses the domain of the result x to deduct new domains for the operands, y or z. For the same constraint, , it means computing a possibly refined interval for y, Y = [y , y u ] ⊆ Y , such that (33) Property (33) is known as the inverse propagation correctness property. Again, taking Y = Y is always possible and sometimes unavoidable. The best we can hope for is to be able to determine the smallest such set, i.e., satisfying (34) Property (34) is called the inverse propagation optimality property. Satisfying this last property can be very difficult.

The Boolean domain for NaN
From now on, we will consider floating-point intervals with boundaries in F. They allow for capturing the extended numbers in F only: NaNs (quiet NaNs and signaling NaNs) should be tracked separately. To this purpose, a Boolean domain N def = { , ⊥}, where stands for "may be NaN" and ⊥ means "cannot be NaN", can be used and coupled with the arithmetic filtering algorithms.
Let be an arithmetic constraint over floating-point numbers, and (X, NaN x ), (Y, NaN y ) and (Z, NaN z ) be the variable domains of x, y and z respectively. In practice, the propagation process for such a constraint reaches a fixed point when the combination of refining domains (X , NaN x ), (Y , NaN y ) and (Z , NaN z ) remains the same obtained in the previous iteration. For each iteration of the algorithm we analyze the NaN domain of all the constraint variables in order to define the next propagator action.
The IEEE 754 Standard [24, Section 7.2] lists all combinations of operand values that yield a NaN result. For the arithmetic operations considered in this paper, NaN is returned if any of the operands is NaN. Moreover, addition and subtraction return NaN when infinities are subtracted (e.g., +∞ −∞ or +∞ +∞), and also , , and Thus, direct projections are such that if NaN y = or NaN z = , then also NaN x = ; indirect projections yield NaN y = NaN z = if NaN x = . Moreover, e.g., if , then the direct projection yields NaN x = also if ±∞ ∈ Y and ∓∞ ∈ Z, and the indirect one allows for ±∞ in Y and ∓∞ in Z only if NaN x = , and so on for the other operators.

Filtering algorithms for simple arithmetic constraints
Filtering algorithms for arithmetic constraints are the main focus of this paper. In the next sections, we will propose algorithms realizing optimal direct projections and correct inverse projections for the addition ( ) and division ( ) operations. The reader interested in implementing constraint propagation for all four operations can find the algorithms and results for the missing operations in Appendix A. We report the correctness proofs of the projections for addition in the main text, leaving those for the remaining operations to Appendix B.3.
The filtering algorithms we are about to present are capable of dealing with any set of rounding modes and are designed to distinguish between all different (special) cases in order to be as precise as possible, especially when the variable domains contain symbolic values. Much simpler projections can be designed whenever precision is not of particular concern. Indeed, the algorithms presented in this paper can be considered as the basis for finding a good trade-off between efficiency and the required precision.

Addition
Here we deal with constraints of the form Thanks to Proposition 2, any set of rounding modes S ⊆ R can be mapped to a pair of "worst-case rounding modes" which, in addition, are never round-to-zero. Therefore, the projection algorithms use the selectors presented in Definition 7 to choose the appropriate worst-case rounding mode, and then operate as if it was the only one in effect, yielding results implicitly valid for the entire set S. It can be proved that Algorithm 1 computes a correct and optimal direct projection, as stated by its postconditions.

Theorem 1 Algorithm 1 satisfies its contract.
Proof Given the constraint x = y S z with x ∈ X = [x , x u ], y ∈ Y = [y , y u ] and z ∈ Z = [z , z u ], Algorithm 1 sets X = [x , x u ] ∩ X; hence, we have X ⊆ X. Moreover, by Proposition 2, for each y ∈ Y , z ∈ Z and r ∈ S, we have y r z y r z y r u z, and because a b implies a ≤ b for any a, b ∈ F according to Definition 3, we know that y r z ≤ y r z ≤ y r u z. Thus, by monotonicity of , we have y r z ≤ y r z ≤ y r z ≤ y r u z ≤ y u r u z u . Therefore, we can focus on finding a lower bound for y r z and an upper bound for y u r u z u . Fig. 2 Direct projection of addition: the function da (resp., da u ); values for y (resp., y u ) on rows, values for z (resp., z u ) on columns Such bounds are given by the functions da and da u of Fig. 2. Almost all of the cases reported in the tables can be trivially derived from the definition of the addition operation in the IEEE 754 Standard [24]; just two cases need further explanation. Concerning the entry of da in which y = −∞ and z = +∞, note that z = +∞ implies z u = +∞. Then for any y > y = −∞, y +∞ = +∞. On the other hand, by the IEEE 754 Standard [24], −∞ +∞ is an invalid operation. For the symmetric case, i.e., the entry of da u in which y u = −∞ and z u = +∞, we can reason dually.
We are now left to prove that ∀X ⊂ X : ∃r ∈ S, y ∈ Y, z ∈ Z : y r z ∈ X . Let us focus on the lower bound x , proving that there always exists a r ∈ S such that y r z = x .
First, consider the cases in which y ∈ (R − ∪ R + ) or z ∈ (R − ∪ R + ). In these cases, a case analysis proves that da (y , z , r ) is equal to y r z . Indeed, if either of the operands (say y ) is −∞ and the other one (say z ) is not +∞, then according to the IEEE 754 Standard we have y r z = −∞ for any r ∈ R. Symmetrically, y r z = +∞ if one operand is +∞ and the other one is not −∞. If, w.l.o.g., y = +∞ and z = −∞, the set X is non-empty only if z u = −∞, and y r z u = +∞ for any r ∈ R.
For the cases in which y ∈ (R − ∪ R + ) and z ∈ (R − ∪ R + ) we have x = y r z , by definition of da of Fig. 2. Remember that, by Proposition 2, there exists r ∈ S such that y r z = y r z . Since y ∈ Y and z ∈ Z, we can conclude that for any X ⊆ X , x ∈ X implies y r z ∈ X .
An analogous argument allows us to conclude that there exists an r ∈ S for which the following holds: for any X ⊆ X , x u ∈ X implies y u r z u ∈ X .
The following example will better illustrate how the tables in Fig. 2 should be used to compute functions da and da u . All examples in this section refer to the IEEE 754 binary single precision format. 8], and that the selected rounding mode is r = r u =↓. In order to compute the lower bound x of X , the new interval for x, function da (+0, −0, ↓) is called. These arguments fall in case a 1 , which yields −0 with rounding mode ↓. Indeed, when the rounding mode is ↓, the sum of −0 and +0 is −0, which is clearly the lowest result that can be obtained with the current choice of Y and Z. For the upper bound x u , the algorithm calls da u (5, 8, ↓). This falls in the case in which both operands are positive numbers (y u , z u ∈ R + ), and therefore x u = y u r u z u = 13. In conclusion, the new interval for x is X = [−0, 13].
If any other rounding mode was selected (say, r = r u = n), the new interval computed by the projection would have been X = [+0, 13].
Inverse propagation For inverse propagation, i.e., the process that infers a new interval for y (or for z) starting from the interval from x and z (x and y, resp.) we define Algorithm 2 and functions ia in Fig. 3 and ia u in Fig. 4, where ≡ indicates the syntactic substitution of expressions. Since the inverse operation of addition is subtraction, note that the values of x and z that minimize y are x and z u ; analogously, the values of x and z that maximize y are x u and z .
When the round-to-nearest rounding mode is in effect, addition presents some nice properties. Indeed, several expressions for lower and upper bounds can be easily computed without approximations, using floating-point operations. In more detail, it can be shown (see the proof of Theorem 2) that when x is subnormal ∇ n+ 2 (x) and ∇ n− 2 (x) are negligible. This allows us to define tight bounds in this case. On the contrary, when the terms ∇ n− 2 (x ) and ∇ n+ 2 (x u ) are non negligible, we need to approximate the values of expressions e and e u . This can always be done with reasonable efficiency [27], but we leave this as an imple- The next result assures us that our algorithm computes a correct inverse projection, as claimed by its postcondition.

Theorem 2 Algorithm 2 satisfies its contract.
Proof Given the constraint x = y S z with x ∈ X = [x , x u ], y ∈ Y = [y , y u ] and z ∈ Z = [z , z u ], Algorithm 2 computes a new and refined domain Y for variable y.
First, observe that the newly computed interval [y , y u ] is either intersected with the old domain Y , so that Y = [y , y u ] ∩ Y , or set to Y = ∅. Hence, Y ⊆ Y holds.
Proposition 2 and the monotonicity of allow us to find a lower bound for y by exploiting the constraint y r z u = x , and an upper bound for y by exploiting the constraint Fig. 3 Inverse projection of addition: function ia Fig. 4 Inverse projection of addition: function ia u y r u z = x u . We will now prove that the case analyses of functions ia , described in Fig. 3, and ia u , described in Fig. 4, express such bounds correctly.
Concerning the operand combinations in which ia takes the value described by the case analysis a 4 , remember that, by the IEEE 754 Standard [24], whenever the sum of two operands with opposite sign is zero, the result of that sum is +0 in all roundingdirection attributes except roundTowardNegative: in that case the result is −0. Then, since z u ↓ (−z u ) = −0, whenr =↓, y can safely be set to succ(−z u ).
As for the case in which ia takes one of the values determined by a 5 , the IEEE 754 Standard [24] asserts that +0 ↓ +0 = +0, while −0 ↓ +0 = −0: the correct lower bound for y is y = +0, in this case. As we already pointed out, for any other rounding-direction attribute +0 −0 = +0 holds, which allows us to include −0 in the new domain.
Concerning cases of ia that give the result described by the case analysis a 6 , we clearly must have y = +∞ ifr =↓; ifr =↑, it should be y + z u > f max and thus y > f max − z u and, by (28) of Proposition 5, y succ(f max ↓ z u ). Ifr = n, there are two cases: In this case, y must be greater than f max , since f max + z u < f max + ∇ n+ 2 (f max )/2 implies that f max n z u = f max < +∞. Note that in this case ∇ n+ Since odd(f max ), for x = +∞ we need y to be greater than or equal to f max + ∇ n+ allows us to apply (29) of Proposition 5, concluding y f max ↑ ∇ n+ 2 (f max )/2 ↑ z u . Equality (35) holds because either the application of ' ↑ ' is exact or the application of ' ↑ ' is exact. In fact, since z u = m·2 e ≥ ∇ n+ 2 (f max )/2 = 2 e max −p , for some 1 ≤ m < 2, there are two cases: either e = e max or e max − p ≤ e < e max .
Suppose first that e = e max : we have and thus Since if e = e max the application of ' ↑ ' is not exact, we prove that the application of ' ↑ ' is exact. Hence, if m > 1, we prove that . It is worth noting that 2 k (2 − m) can be represented by a normalized mantissa; moreover, since 1 ≤ k ≤ p − 1, e min ≤ e max − k ≤ e max , hence, and, also in this case, f max + (∇ n+ 2 (f max )/2 ↑ z u ) ∈ F. Suppose now that e max −p ≤ e < e max and let h def = e−e max +p so that 0 ≤ h ≤ p−1.
In this case we show that the application of ' ↑ ' is exact. Indeed, we have If e = e max − p and m = 1, then h = 0, m − 2 −h = 0 and thus ∇ n+ which is an element of F.
Dual arguments w.r.t. the ones used to justify cases of ia that give the result described by a 4 , a 6 and a 5 can be used to justify the cases of ia u described by a 9 , a 10 and a 7 .
We now tackle the entries of ia described by a 3 , and those of ia u described by a 8 . Exploiting x y z and x y z, by Proposition 4, we have ifr = n and even(x); > x + ∇ n− 2 (x)/2, ifr = n and odd(x).
We can now exploit Proposition 5 and obtain ifr =↓ and x = z u ; succ pred(x ) ↓ z u , ifr =↑; In fact, if x = z u , then, according to IEEE 754 [24, Section 6.3], for each non-NaN, nonzero and finite w ∈ F, −0 is the least value for y that satisfies w = y ↓ w. If x = z u , then case (29) of Proposition 5 applies and we have y x ↑ z u . Suppose now that pred(x ) = z u , then pred(x ) ↓ z u ≡ 0 and succ pred(x ) ↓ z u = f min , coherently with the fact that, for each non-NaN, nonzero and finite w ∈ F, f min is the least value for y that satisfies w = y ↑ pred(w). If pred(x ) = z u , then case (26) of Proposition 5 applies and we have y succ pred(x ) ↓ z u . A symmetric argument justifies (39).
For the remaining cases, we first show that when ∇ n+ 2 (x) = f min , The previous equality has the following main consequences: we can perform the computation in F, that is, we do not need to compute ∇ n+ 2 (x)/2 and, since [ , we can apply (30) of Proposition 5, obtaining a tight bound for y u .
Let us prove (40). Suppose ∇ n+ 2 (x u ) = f min , and assume x u = z . There are two cases: In the case where  whereas, for the case wherer u = n, we have Note that −f max is the lowest value that variable y could take, since there is no value for z ∈ Z that summed with −∞ gives a value in X. Indeed, if we take z = +f max , then we have −f max r +f max = +0 ∈ X for any r ∈ R. On the other hand, +∞ is clearly the highest value y could take, since +∞ r z = +∞ ∈ X for any value of z ∈ Z\{−∞}. In this case, our projections yield a more refined result than the competing tool FPSE [10], which computes the wider interval Y = [−∞, +∞]. One of the reasons the inverse projection for addition is not optimal is because floating point numbers present some peculiar properties that are not related in any way to those of real numbers. For interval-based consistency approaches, [29] identified a property of the representation of floating-point numbers and proposed to exploit it in filtering algorithms for addition and subtraction constraints. In [4,5] some of these authors revised and corrected the Michel and Marre filtering algorithm on intervals for addition/subtraction constraints under the round to nearest rounding mode. A generalization of such algorithm to the all rounding modes should be used to enhance the precision of the classical inverse projection of addition. Indeed, classical and maximum ULP filtering [5] for addition are orthogonal: both should be applied in order to obtain optimal results. Therefore, inverse projections for addition, as the one proposed above, have to be intersected with a filter based on the Michel and Marre property in order to obtain more precise results. By applying maximum ULP filtering [5,29], we obtain the much tighter intervals Y, Z = [−1.1 · · · 1 × 2 24 , 1.0 × 2 25 ]. These are actually optimal as −1.1 · · · 1 × 2 24 S 1.0 × 2 25 = 1.0 × 2 25 S −1.1 · · · 1 × 2 24 = 2.0. This example shows that filtering by maximum ULP can be stronger than our interval-consistency based filtering. However, the opposite phenomenon is also possible. Consider again X = [1.0, 2.0] and Z = [1.0, 5.0]. Filtering by maximum ULP projection gives Z = [−1.1 · · · 1 × 2 24 , 1.0 × 2 25 ]; in contrast, our inverse projection exploits the available information on Z to obtain Y = [−4, 1.0 · · · 01]. As we already stated, our filtering and maximum ULP filtering should both be applied in order to obtain precise results.
Exploiting the commutative property of addition, the refinement Z of Z can be defined analogously.

Division
In this section we deal with constraints of the form x = y S z with S ⊆ R.
Direct propagation For direct propagation, interval Z is partitioned into the sign- . This is needed because the sign of operand z determines the monotonicity with respect to y, and therefore the interval bounds to be used for propagation depend on it. Hence, once Z has been partitioned into sign-homogeneous intervals, we use the interval Y and W = Z − , to obtain the new interval [x − , x − u ], and Y and W = Z + , to obtain [x + , x + u ]. The appropriate bounds for interval propagation are chosen by function τ of Fig. 5. Note that the sign of z is, by construction, constant over interval W . The selected values are then taken as arguments by functions dd and dd u of Fig. 6, which return the correct bounds for the aforementioned new intervals for X. The intervals X ∩ [x − , x − u ] and X ∩ [x + , x + u ] are eventually joined using convex union, denoted by , to obtain the refining interval X . It can be proved that Algorithm 3 computes a correct and optimal direct projection, as ensured by its postconditions.   The following result assures us that Algorithm 4 computes a correct first inverse projection, as ensured by its postcondition.

Theorem 4 Algorithm 4 satisfies its contract.
Once again, in order to obtain more precise results in some cases, the first inverse projection for division has to be intersected with a filter based on an extension of the Michel and Marre property originally proposed in [29] and extended to multiplication and division in [5]. Indeed, when interval X does not contain zeroes and interval Z contains zeros and infinities, the proposed filtering by maximum ULP algorithm is able to derive more precise bounds than the ones obtained with the inverse projection we are proposing. Thus, for division (and for multiplication as well), the indirect projection and filtering by maximum ULP are mutually exclusive: one applies when the other cannot derive anything useful [5]. Inverse propagation second projection The second inverse projection for division computes a new interval for operand z. For this projection, we need to partition interval X into sign-homogeneous intervals X − def = X ∩ [−∞, −0] and X + def = X ∩ [+0, +∞] since, in this case, it is the sign of X that matters for deriving correct bounds for Z. Once X has been partitioned, we use intervals X − and Y to obtain the interval [z − , z − u ]; intervals X + and Y to obtain [z + , z + u ]. The new bounds for z are computed by functions id s of Fig. 10 and id s u of Fig. 11, after the appropriate interval extrema of Y and V = X − (or V = X + ) have been selected by function τ . The intervals Z ∩ [z − , z − u ] and Z ∩ [z + , z + u ] will be then joined with convex union to obtain Z .
In order to obtain more precise results, the result of our second inverse projection can also be intersected with the interval obtained by the maximum ULP filter proposed in [5]. Indeed, when interval X does not contain zeros and interval Y contains zeros and infinities, the proposed filtering by maximum ULP algorithm is able to derive tighter bounds than those obtained with the inverse projection presented in this work.

Experimental evaluation
The main aim of this section is to motivate the need of provably correct filtering algorithms, by highlighting the issues caused by the unsoundness of most available implementations of similar methods.

Software verification
As we reported in Section 1.3, we implemented our work in the commercial tool ECLAIR. While the initial results on a wide range of self-developed tests looked very promising, we wanted to compare them with the competing tools presented in the literature, in order to better assess the strength of our approach with respect to the state of the art. Unfortunately, most of these tools were either unavailable, or not sufficiently equipped to analyze real-world C/C++ programs. We could, however, do a comparison with the results obtained in [38]. It presents a tool called seVR-fpe, for floating-point exception detection based on symbolic execution and value-range analysis. The same task can be carried out by the constraint-based symbolic model checker we included in ECLAIR. The authors of seVR-fpe tested their tool both on a self-developed benchmark suite and on real-world programs. Upon contacting them, they were unfortunately unable to provide us with more detailed data regarding their analysis of real world programs. This prevents us from doing an in-depth comparison of the tools, since we only know the total number of bugs found, but not their exact nature and location. Data with this level of detail was instead available for (most of) their selfdeveloped benchmarks. The results obtained by running ECLAIR on them are reported in Table 1. ECLAIR could find a number of possible bugs significantly higher than seVR-fpe. As expected, due to the provable correctness of the algorithms employed in ECLAIR, no false positives were detected among the inputs it generated. This confirms the solid results obtainable by means of the algorithms presented in this paper.

SMT solvers
We compare ECLAIR with several SMT solvers that support floating-point arithmetic by executing them on a benchmark suite devised to test their floating-point theory for soundness. Since our aim is to evaluate filtering algorithms for floating-point addition, subtraction, multiplication and division, we only included tests that do not rely on other theories, such as arrays, bit-vectors and uninterpreted functions, as well as those containing quantifiers. Such features are typical of SMT, and are tackled by techniques which are out of the scope of this paper. The suite is made of a total of 151,432 tests, of which: -10,380 are randomly generated tests by Florian Schanda. (random directory from the benchmark suite 6 used in [12]). The experiments were carried out on a high-end laptop with an x86 64 CPU (6 cores @2.20GHz) and 16 GB of RAM, running Ubuntu 20.04. Each solver's version is reported in the table. MathSAT has been executed with the option -theory.fp.mode=2, which enables the ACDL-based solver for the floating-point theory.
Results are reported in Table 2. The main observation we can make is that ECLAIR is the only tool based on interval reasoning to be completely sound. On the contrary, Colibri and MathSAT are unsound on numerous tests, even though they did not explicitly report any error. This hinders their use for program verification. On the other hand, bit-blasting based tools CVC4 and Z3 do not present such issues, because their bit-vector encodings for floating-point arithmetic are derived from solid and formally verified circuit designs such as [34]. This demonstrates that the provably-correct filters presented in this paper are needed to achieve reliable implementations of interval-based constraint solving methods.
The execution times seem to be mainly determined by implementation details such as the programming language used. In fact, Colibri and ECLAIR, the slowest tools, were written respectively in ECLiPSe Prolog 9 and SWI Prolog, 10 while other tools were written in C++.

Discussion and conclusion
With the increasing use of floating-point computations in mission-and safety-critical settings, the issue of reliably verifying their correctness has risen to a point in which testing or other informal techniques are not acceptable any more. Indeed, this phenomenon has been fostered by the wide adoption of the IEEE 754 floating-point format, which has significantly simplified the use of floating-point numbers, by providing a precise, sound, and reasonably cross-platform specification of floating-point representations, operations and their semantics. The approach we propose in this paper exploits these solid foundations to enable a wide range of floating-point program verification techniques. It is based on the solution of constraint satisfaction problems by means of interval-based constraint propagation, which is enabled by the filtering algorithms we presented. These algorithms cover the whole range of possible floating-point values, including symbolic values, with respect to interval-based reasoning. Moreover, they not only support all IEEE 754 available rounding-modes, but they also allow to take care of uncertainty on the rounding-mode in use. Some important implementation aspects are also taken into account, by allowing both the use of machine floating-point arithmetic for all computations (for increased performance), and of extendedprecision arithmetic (for better precision with the round-to-nearest rounding mode). In both cases, correctness is guaranteed, so that no valid solutions can erroneously be removed from the constraint system. This is supported by the extensive correctness proofs of all algorithms and tables, which allow us to claim that neither false positives, nor false negatives may be produced. The experimental evaluation of Section 5 shows that soundness is, indeed, a widespread issue in several floating-point verification tools. Our work provides solid foundations to soundly develop such kind of tools.
Several aspects of the constraint-based verification of floating-point programs remain, however, open problems, both from a theoretical and a practical perspective. As we showed throughout the paper, the filtering algorithms we presented are not optimal, i.e., they may not yield the tightest possible intervals containing all solutions to the constraint system. They must be interleaved with the filtering algorithms of [5], and they may require multiple passes before reaching the maximum degree of variable-domain pruning they are capable of. Therefore, the next possible advance in this direction would be conceiving optimal filtering algorithms, that reduce variable domains to intervals as tight as possible with a single application. This has been achieved in [21], but only for addition.
However, filtering algorithms only represent a significant, but to some extent limited, part of the constraint solving process. Indeed, even an optimally pruned interval may contain values that are not solutions to the constraint system, due to the possible non-linearity thereof. If the framework in use supports multi-intervals, this issue is dealt with by means of labeling techniques: when a constraint-solving process reaches quiescence, i.e., the application of filtering algorithms fails to prune variable domains any further, such intervals are split into two or more sub-intervals, and the process continues on each partition separately. In this context, the main issues are where to split intervals, and in how many parts. These issues are currently addressed with heuristic labeling strategies. Indeed, significant improvements to the constraint-propagation process could be achieved by investigating better labeling strategies. To this end, possible advancements would include the identification of objective criteria for the evaluation of labeling strategies on floating point-numbers, and the conception of labeling strategies tailored to the properties of constraint systems most commonly generated by numeric programs.
In conclusion, we believe the work presented in this paper can be an extensive reference for the readers interested in realizing applications for formal reasoning on floating-point computations, as well as a solid foundation for further improvements in the state of the art.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.