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

TL;DR

This paper develops and proves the correctness of algorithms to accurately bound floating-point computations under IEEE 754 standards, aiding formal program verification involving floating-point arithmetic.

Contribution

It introduces the first formally verified filtering algorithms for bounding floating-point variables considering all IEEE 754 rounding modes, even with partial mode knowledge.

Findings

Algorithms precisely bound floating-point variables

Correctness is formally proven for all IEEE 754 rounding modes

Supports verification of floating-point programs with partial rounding mode information

Abstract

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.

Tables

Table 1. Table 1: Number of exceptions found by ECLAIR and seVR-fpe on the self-developed benchmarks of WuLZ17 .

Exception type	ECLAIR	seVR-fpe	Difference
total	135	66	69
overflow	55	26	29
underflow	30	13	17
invalid	47	8	39
divbyzero	3	3	0
false positives	0	15	-15

Table 2. Table 2: Results of the evaluation of SMT-solvers.

Solver	Version	Solved	Errors	Unsound	Time (h:m:s)
Colibri	2176	148,766	0	67	01:02:42.03
CVC4	1.8	148,833	0	0	00:09:32.11
ECLAIR	—	148,833	0	0	01:59:28.80
MathSAT ACDL	5.6.5	147,958	0	875	00:07:05.89
z3	4.8.10	148,833	0	0	00:08:38.33

Equations

$$\mathtt{result}:={\mathtt{operand}}_1\mathtt{operator}{\mathtt{operand}}_2.$$

$$m=(d_0.d_1{d}_2\dots d_{p-1}{)}_2=\sum _{i=0}^{p-1}{d}_i{2}^{-i};$$

$$\mathrm{succ}(x)\frac{\text{buildrel}\mathrm{def}}{=}\{\begin{array}{ll}+\infty ,& \text{if x=fmax;}\\ \min \{y\in \mathbb{F}\mid y>x\},& \text{if −fmax≤x<−fmin}\\ & \text{ or fmin≤x<fmax;}\\ f_{\min },& \text{if x≡0;}\\ -0,& \text{if x=−fmin;}\\ -f_{\max },& \text{if x=−∞;}\\ \text{undefined},& \text{otherwise.}\end{array}$$

$$[x{]}_{\uparrow }$$

$$[x{]}_{\downarrow }$$

$$[x{]}_0$$

$$[x{]}_{\mathrm{n}}$$

$$[x{]}_{\downarrow }\le x\le [x{]}_{\uparrow },$$

$$[x{]}_{\downarrow }\le [x{]}_0\le [x{]}_{\uparrow },$$

$$[x{]}_{\downarrow }\le [x{]}_{\mathrm{n}}\le [x{]}_{\uparrow }.$$

$$[x{]}_{\downarrow }=-[-x{]}_{\uparrow }.$$

$$\mathcal{I}_{\mathbb{F}}\mathrel{\buildrel\mathrm{def}\over{=}}\{\varnothing\}\cup\bigl{\{}\,[\ell,u]\bigm{|}\ell,u\in\mathbb{F},\ell\preccurlyeq u\,\bigr{\}}.$$

$$r_{\ell }(S,y,\text{boxcircle},z)$$

$$r_u(S,y,\text{boxcircle},z)$$

$${\hat{r}}_{\ell }(S,\text{boxcircle},b)$$

$${\hat{r}}_u(S,b)$$

$$\displaystyle\bigl{(}\bar{r}_{\ell}^{\ell}(S,b,\boxcircle,a),\bar{r}_{u}^{\ell}(S,b,\boxcircle,a)\bigr{)}$$

$$\displaystyle\bigl{(}\bar{r}_{\ell}^{r}(S,b,\boxcircle,a),\bar{r}_{u}^{r}(S,b,\boxcircle,a)\bigr{)}$$

$$y{\text{boxcircle}}_{{\mathrm{r}}_{\mathrm{\ell }}}z\preccurlyeq y{\text{boxcircle}}_{\mathrm{r}}z\preccurlyeq y{\text{boxcircle}}_{{\mathrm{r}}_{\mathrm{u}}}z.$$

$$y{\text{boxcircle}}_{{\mathrm{r}}_{\mathrm{\ell }}}z=y{\text{boxcircle}}_{{\mathrm{r}}^{\prime }}z\text{and}y{\text{boxcircle}}_{{\mathrm{r}}_{\mathrm{u}}}z=y{\text{boxcircle}}_{{\mathrm{r}}^{\prime \prime }}z.$$

$${\hat{y}}^{\prime }\preccurlyeq y\preccurlyeq {\tilde{y}}^{\prime \prime }.$$

$${\nabla }^{\downarrow }(x)$$

$${\nabla }^{\uparrow }(x)$$

$${\nabla }_2^{\mathrm{n}-}(x)$$

$${\nabla }_2^{\mathrm{n}+}(x)$$

$$\displaystyle\min_{x_{\ell}\leq x\leq x_{u}}\bigl{(}x+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)/2\bigr{)}=x_{\ell}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x_{\ell})/2,$$

$$\displaystyle\max_{x_{\ell}\leq x\leq x_{u}}\bigl{(}x+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)/2\bigr{)}=x_{u}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x_{u})/2.$$

$$(\mathrm{succ}(-f_{\max })-f_{\max })/2=-2^{e_{\max }}(2-2^{-p})=x_{\ell }+{\nabla }_2^{\mathrm{n}-}(x_{\ell })/2,$$

$$x\preccurlyeq y{\text{boxcircle}}_{\downarrow }z$$

$$x\preccurlyeq y{\text{boxcircle}}_{\uparrow }z$$

$$x\preccurlyeq y{\text{boxcircle}}_{\mathrm{n}}z$$

$$x\succcurlyeq y{\text{boxcircle}}_{\downarrow }z$$

$$x\succcurlyeq y{\text{boxcircle}}_{\uparrow }z$$

$$x\succcurlyeq y{\text{boxcircle}}_{\mathrm{n}}z$$

$$[[e{]]}_{\downarrow }$$

$$[[e{]]}_{\uparrow }$$

$$x\ge e$$

$$\text{if [[e]]↓=+∞, }x>e$$

$$x\le e$$

$$\text{if [[e]]↓=−∞, }x<e$$

$$x\ge e$$

Full text

††supplemental: Appendices

1 Introduction

Programs using floating-point numbers are notoriously difficult to reason about Monniaux08 . Many factors complicate the task:

compilers may transform the code in a way that does not preserve the semantics of floating-point computations; 2. 2.

floating-point formats are an implementation-defined aspect of most programming languages; 3. 3.

there are different, incompatible implementations of the operations for the same floating-point format; 4. 4.

mathematical libraries often come with little or no guarantee about what is actually computed; 5. 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. 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 DelmasGPSTV09 .

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:

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. 2.

Even though the used floating-point formats are implementation-defined aspects of, say, C and C**++**111This is not relevant if one analyzes machine or sufficiently low-level intermediate code. the wide adoption of the IEEE 754 standard for binary floating-point arithmetic IEEE-754-2008 has improved things considerably. 3. 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 CuytKVV02 .222The indeterminacy described in CuytKVV02 is present also in the 2008 edition of IEEE 754 IEEE-754-2008 . 4. 4.

A pragmatic, yet effective approach to support formal reasoning on commonly used implementation of mathematical functions has been recently proposed in BagnaraCGB21 . 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. 5.

A static analysis for detecting floating-point exceptions based on abstract interpretation has been presented in Mine04 . A few attempts at this task have been made using other techniques BarrVLS13 ; WuLZ17 but, as we argue in Sections 1.4 and 5, they present precision and soundness issues. 6. 6.

Most verification approaches in the literature assume the round-to-nearest rounding mode BotellaGM06 , or over-approximate by always considering worst-case rounding modes Mine04 . Analyses based on SMT solvers BrainTRW15 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., Rump13 ; RumpO07 );

• •

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 Fiedler10 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 Watte08 .

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 MichelRL01 ; Michel02 , 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 Intelx8664 and ARM Armv8 — 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.

1.1 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

[TABLE]

A further refinement is the computation of the static single assignment form (SSA) AhoLSU06 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 $\mathtt{z}\mathrel{:=}\mathtt{z}*\mathtt{z}+\mathtt{z}$ is $\mathtt{t}\mathrel{:=}\mathtt{z}*\mathtt{z};\;\mathtt{z}\mathrel{:=}\mathtt{t}+\mathtt{z}$, which in an SSA form becomes $\mathtt{t}_{1}\mathrel{:=}\mathtt{z}_{1}*\mathtt{z}_{1};\;\mathtt{z}_{2}\mathrel{:=}\mathtt{t}_{1}+\mathtt{z}_{1}$. These, in turn, can be regarded as the conjunction of the constraints $t_{1}=z_{1}\boxdot z_{1}$ and $z_{2}=t_{1}\boxplus z_{1}$, where by $\boxdot$ and $\boxplus$ 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 King76 ; ClarkeR85 . 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 BotellaGM06 ; BagnaraCGG13ICST .

1.2 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\boxdot x=x$, then any value outside the set $\{\mathrm{NaN},+0,1,+\infty\}$ 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}\boxplus 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.

1.3 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 floating-point 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.

Definition 1

(Absorption.) Let $x,y,z\in\mathbb{F}$ with $y,z\in\mathbb{R}$, let $\mathord{\boxcircle}$ be any IEEE 754 floating-point operator, and let $x=y\boxcircle z$. Then $y\boxcircle z$ gives rise to absorption if

• •

$\mathord{\boxcircle}=\mathord{\boxplus}$ and either $x=y$ and $z\neq 0$, or $x=z$ and $y\neq 0$;

• •

$\mathord{\boxcircle}=\mathord{\boxminus}$ and either $x=y$ and $z\neq 0$, or $x=-z$ and $y\neq 0$;

• •

$\mathord{\boxcircle}=\mathord{\boxdot}$ and either $x=\pm y$ and $z\neq\pm 1$, or $x=\pm z$ and $y\neq\pm 1$;

• •

$\mathord{\boxcircle}=\mathord{\boxslash}$, $x=\pm y$ and $z\neq\pm 1$.

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,333https://bugseng.com/eclair, last accessed on October 28th, 2021. 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 CousotC77 , automatic generation of test-cases, and symbolic model checking. The latter two are based on symbolic execution and constraint satisfaction problems GotliebBR98 ; GotliebBR00 , whose solution is based on multi-interval 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 BagnaraCGB21 .

To demonstrate the capabilities of the techniques presented in this paper, we applied them to the C code excerpt of Figure 1. It is part of the implementation of the Bessel functions in the GNU Scientific Library,444https://www.gnu.org/software/gsl/, last accessed on October 28th, 2021. 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 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 Figure 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 $\mathtt{x}=\texttt{-0x1.8p-1021}\approx-6.6752\times 10^{-308}$. If the function is executed with $\mathtt{x}=\texttt{-0x1p+1023}\approx-8.9885\times 10^{307}$, the multiplication of line 29 yields a negative infinity. Since $\mathtt{ax}=|\mathtt{x}|$, we know $\mathtt{x}=\texttt{0x1p+1023}$ would also cause the overflow. The same value of $\mathtt{x}$ causes an overflow in line 30 as well. The division in the same line produces a NaN if the function is executed with $\mathtt{x}=-\infty$.

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 Figure 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 $\mathtt{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.

1.4 Related Work

1.4.1 Filtering Algorithms

In Michel02 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 BotellaGM06 , the authors extended the approach of Michel02 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 BotellaGM06 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, Gallois-WongBC20 presented optimal inverse projections for addition under the round-to-nearest 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 MarreM10 . The authors are able to prove the optimality of the lower bounds computed by their algorithms. However, Gallois-WongBC20 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 MarreM10 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 BagnaraCGG13ICST ; BagnaraCGG16IJOC ). In this paper we discuss the cases in which the filtering algorithms in BagnaraCGG13ICST ; BagnaraCGG16IJOC ; MarreM10 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.

1.4.2 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 BarrettT18 . The need for verifying floating-point programs lead to the introduction of a floating-point theory BrainTRW15 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 BarrettCDHJKRT11 ; BrainSS19 , MathSAT CimattiGSS13 and Z3 MouraB08 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 BrainDGHK14 . Colibri MarreBC17 uses constraint programming techniques, with filtering algorithms such as those in BotellaGM06 ; BagnaraCGG16IJOC and those presented in this paper. However, MarreBC17 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 BrainSS19 .

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.

1.4.3 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 Mine04 and implemented in the tool Astrée.555https://www.absint.com/astree/index.htm, last accessed on October 28th, 2021. 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, Mine04 does not offer a treatment of symbolic values (e.g., infinities) as exhaustive as the one we offer in this paper.

In BarrVLS13 the authors proposed a symbolic execution system for detecting floating-point 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 WuLZ17 , 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 floating-point exception detection presented in WuLZ17 is compared with the same analysis based on our propagation algorithms. As expected, no false positives were detected among the results of our analysis.

1.5 Contribution

This paper improves the state of the art in several directions:

all rounding modes are treated and there is no assumption that the rounding mode in effect is known and unchangeable (increased generality); 2. 2.

utilization, to a large extent, of machine floating-point arithmetic in the analyzer with few rounding mode changes (increased performance); 3. 3.

accurate treatment of round half to even —the default rounding mode of IEEE 754— (increased precision); 4. 4.

explicit and complete treatment of intervals containing symbolic values (i.e., infinities and signed zeros); 5. 5.

application of floating-point constraint propagation techniques to enable detection of program anomalies such as overflows, underflows, absorption, generation of NaNs.

1.6 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.

2 Preliminaries

We will denote by $\mathbb{R}_{+}$ and $\mathbb{R}_{-}$ the sets of strictly positive and strictly negative real numbers, respectively. The set of affinely extended reals, $\mathbb{R}\cup\{-\infty,+\infty\}$, is denoted by $\overline{\mathbb{R}}$.

Definition 2

(IEEE 754 binary floating-point numbers.)

A set of IEEE 754 binary floating-point numbers IEEE-754-2008 is uniquely identified by: $p\in\mathbb{N}$, the number of significant digits (precision); ${e_{\mathrm{max}}}\in\mathbb{N}$, the maximum exponent, the minimum exponent being ${e_{\mathrm{min}}}\mathrel{\buildrel\mathrm{def}\over{=}}1-{e_{\mathrm{max}}}$. The set of binary floating-point numbers $\mathbb{F}(p,{e_{\mathrm{max}}},{e_{\mathrm{min}}})$ includes:

• •

all signed zero and non-zero numbers of the form $(-1)^{s}\cdot 2^{e}\cdot m$, where

• –

$s$ is the sign bit;

• –

the exponent $e$ is any integer such that ${e_{\mathrm{min}}}\leq e\leq{e_{\mathrm{max}}}$;

• –

the mantissa $m$, with $0\leq m<2$, is a number represented by a string of $p$ binary digits with a “binary point” after the first digit:

[TABLE]

• •

the infinities $+\infty$ and $-\infty$; the NaNs: $\mathrm{qNaN}$ (quiet NaN) and $\mathrm{sNaN}$ (signaling NaN).

Numbers such that $d_{0}=1$ are called normal. The smallest positive normal floating-point number is ${f^{\mathrm{nor}}_{\mathrm{min}}}\mathrel{\buildrel\mathrm{def}\over{=}}2^{{e_{\mathrm{min}}}}$ and the largest is ${f_{\mathrm{max}}}\mathrel{\buildrel\mathrm{def}\over{=}}2^{{e_{\mathrm{max}}}}(2-2^{1-p})$. The non-zero floating-point numbers such that $d_{0}=0$ are called subnormal: their absolute value is less than $2^{{e_{\mathrm{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_{\mathrm{min}}}\mathrel{\buildrel\mathrm{def}\over{=}}2^{{e_{\mathrm{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 $\mathop{\mathrm{even}}\nolimits(x)$ (resp., $\mathop{\mathrm{odd}}\nolimits(x)$) to signify that the least significant digit of $x$’s mantissa, $d_{p-1}$, is [math] (resp., $1$).

In the sequel we will only be concerned with IEEE 754 binary floating-point numbers and we will write simply $\mathbb{F}$ for $\mathbb{F}(p,{e_{\mathrm{max}}},{e_{\mathrm{min}}})$ when there is no risk of confusion.

Definition 3

(Floating-point symbolic order.)

Let $\mathbb{F}$ be any IEEE 754 floating-point format. The relation $\mathord{\prec}\subseteq\mathbb{F}\times\mathbb{F}$ is such that, for each $x,y\in\mathbb{F}$, $x\prec y$ if and only if both $x$ and $y$ are not NaNs and either: $x=-\infty$ and $y\neq-\infty$, or $x\neq+\infty$ and $y=+\infty$, or $x=-0$ and $y\in\{+0\}\cup\mathbb{R}_{+}$, or $x\in\mathbb{R}_{-}\cup\{-0\}$ and $y=+0$, or $x,y\in\mathbb{R}$ and $x<y$. The partial order $\mathord{\preccurlyeq}\subseteq\mathbb{F}\times\mathbb{F}$ is such that, for each $x,y\in\mathbb{F}$, $x\preccurlyeq y$ if and only if both $x$ and $y$ are not NaNs and either $x\prec y$ or $x=y$.

Note that $\mathbb{F}$ without the NaNs is linearly ordered with respect to ‘$\prec$’.

For $x\in\mathbb{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 [math]. Thus, when we write, e.g., $x<y$ we mean that $x$ is numerically less than $y$ (for example, we have $-0\prec+0$ though $-0\nless+0$, but note that $x\preccurlyeq y$ implies $x\leq y$). Numerical equivalence will be denoted by ‘$\equiv$’ so that $x\equiv 0$, $x\equiv+0$ and $x\equiv-0$ all denote $(x=+0)\lor(x=-0)$.

Definition 4

(Floating-point predecessors and successors.)

The partial function $\mathord{\mathop{\mathrm{succ}}\nolimits}\colon\mathbb{F}\rightarrowtail\mathbb{F}$ is such that, for each $x\in\mathbb{F}$,

[TABLE]

The partial function $\mathord{\mathop{\mathrm{pred}}\nolimits}\colon\mathbb{F}\rightarrowtail\mathbb{F}$ is defined by reversing the ordering, so that, for each $x\in\mathbb{F}$, $\mathop{\mathrm{pred}}\nolimits(x)=-\mathop{\mathrm{succ}}\nolimits(-x)$ whenever $\mathop{\mathrm{succ}}\nolimits(x)$ is defined.

Let $\circ\in\{\mathord{+},\mathord{-},\mathord{\cdot},\mathord{/}\}$ denote the usual arithmetic operations over the reals. Let $R\mathrel{\buildrel\mathrm{def}\over{=}}\{{\mathord{\downarrow}},{0},{\mathord{\uparrow}},{\mathrm{n}}\}$ denote the set of IEEE 754 rounding modes (rounding-direction attributes): round towards minus infinity (roundTowardNegative, ${\mathord{\downarrow}}$), round towards zero (roundTowardZero, ${0}$), round towards plus infinity (roundTowardPositive, ${\mathord{\uparrow}}$), and round to nearest (roundTiesToEven, ${\mathrm{n}}$). We will use the notation $\mathord{\boxcircle_{r}}$, where $\mathord{\boxcircle}\in\{\mathord{\boxplus},\mathord{\boxminus},\mathord{\boxdot},\mathord{\boxslash}\}$ and $r\in R$, to denote an IEEE 754 floating-point operation with rounding $r$.

The rounding functions are defined as follows. Note that they are not defined for [math]: the IEEE 754 standard, in fact, for operations whose exact result is [math], bases the choice between $+0$ and $-0$ on the operation itself and on the sign of the arguments (IEEE-754-2008, , Section 6.3).

Definition 5

(Rounding functions.)

The rounding functions defined by IEEE 754, $\mathord{[\cdot]_{{\mathord{\uparrow}}}}\colon\mathbb{R}\setminus\{0\}\to\mathbb{F}$, $\mathord{[\cdot]_{{\mathord{\downarrow}}}}\colon\mathbb{R}\setminus\{0\}\to\mathbb{F}$, $\mathord{[\cdot]_{{0}}}\colon\mathbb{R}\setminus\{0\}\to\mathbb{F}$ and $\mathord{[\cdot]_{{\mathrm{n}}}}\colon\mathbb{R}\setminus\{0\}\to\mathbb{F}$, are such that, for each $x\in\mathbb{R}\setminus\{0\}$,

[TABLE]

Note that, when the result of an operation has magnitude lower than ${f^{\mathrm{nor}}_{\mathrm{min}}}$, it is rounded to a subnormal number, by adjusting it to the form $(-1)^{s}\cdot 2^{e_{\mathrm{min}}}\cdot m$, and truncating its mantissa $m$, which now starts with at least one [math], to the first $p$ digits. This phenomenon is called gradual underflow, and while it is preferred to hard underflow, which truncates a number to [math], it still may cause precision issues due to the reduced number of significant digits of subnormal numbers.

The rounding modes ${\mathord{\downarrow}}$ and ${\mathord{\uparrow}}$ are the most “extreme”, while ${\mathrm{n}}$ and ${0}$ are always contained between them. We formalize this observation as follows:

Proposition 1

(Properties of rounding functions.)*

Let $x\in\mathbb{R}\setminus\{0\}$. Then*

[TABLE]

Moreover,

[TABLE]

In this paper, we use intervals of floating-point numbers in $\mathbb{F}$ that are not NaNs.

Definition 6

(Floating-point intervals.) Let $\mathbb{F}$ be any IEEE 754 floating-point format. The set $\mathcal{I}_{\mathbb{F}}$ of floating-point intervals with boundaries in $\mathbb{F}$ is given by

[TABLE]

By $[\ell,u]$ we denote the set $\{\,x\in\mathbb{F}\mid\ell\preccurlyeq x\preccurlyeq u\,\}$. The set $\mathcal{I}_{\mathbb{F}}$ is a bounded meet-semilattice with least element $\varnothing$, greatest element $[-\infty,+\infty]$, and the meet operation, which is induced by set-intersection, will be simply denoted by $\mathord{\cap}$.

Floating-point intervals with boundaries in $\mathbb{F}$ allow to capture the extended numbers in $\mathbb{F}$: NaNs should be tracked separately.

3 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.

3.1 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_{\ell}$) and upper (function $r_{u}$) bounds of an operation in case of direct projections.

Definition 7

(Rounding mode selectors for direct projections.)

Let $\mathbb{F}$ be any IEEE 754 floating-point format and $S\subseteq R$ be a set of rounding modes. Let also $y,z\in\mathbb{F}$ and $\mathord{\boxcircle}\in\{\mathord{\boxplus},\mathord{\boxminus},\mathord{\boxdot},\mathord{\boxslash}\}$ be such that either $\mathord{\boxcircle}\neq\mathord{\boxslash}$ or $z\neq 0$. Then

[TABLE]

The following functions select the rounding modes that will be used for the lower (functions $\bar{r}_{\ell}^{r}$ and $\bar{r}_{\ell}^{\ell}$) and upper (functions $\bar{r}_{u}^{r}$ and $\bar{r}_{u}^{\ell}$) 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: $\bar{r}_{\ell}^{r}$ and $\bar{r}_{u}^{r}$ for the right one, $\bar{r}_{\ell}^{\ell}$ and $\bar{r}_{u}^{\ell}$ for the left one.

Definition 8

(Rounding mode selectors for inverse projections.)

Let $\mathbb{F}$ be any IEEE 754 floating-point format and $S\subseteq R$ be a set of rounding modes. Let also $a,b\in\mathbb{F}$ and $\mathord{\boxcircle}\in\{\mathord{\boxplus},\mathord{\boxminus},\mathord{\boxdot},\mathord{\boxslash}\}.$ First, we define

[TABLE]

Secondly, we define the following selectors:

[TABLE]

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 $\mathbb{F}$, $S$, $y$, $z$ and ‘$\mathord{\boxcircle}$’ be as in Definition 7. Let also $r_{\ell}=r_{\ell}(S,y,\mathord{\boxcircle},z)$ and $r_{u}=r_{u}(S,y,\mathord{\boxcircle},z)$. Then, for each $r\in S$

[TABLE]

Moreover, there exist $r^{\prime},r^{\prime\prime}\in S$ such that

[TABLE]

Now, consider $x=y\boxcircle_{\mathrm{r}}z$ with $x,z\in\mathbb{F}$ and $r\in S$. Let $\bar{r}_{\ell}=\bar{r}_{\ell}^{\ell}(S,x_{u},\boxcircle,z)$ and $\bar{r}_{u}=\bar{r}_{u}^{\ell}(S,x_{\ell},\boxcircle,z)$, according to Definition 8. Moreover, let $\hat{y}^{\prime}$ be the minimum $y^{\prime}\in\mathbb{F}$ such that $x=y^{\prime}\boxcircle_{\mathrm{\bar{r}_{\ell}}}z$, and let $\tilde{y}^{\prime\prime}$ be the maximum $y^{\prime\prime}\in\mathbb{F}$ such that $x=y^{\prime\prime}\boxcircle_{\mathrm{\bar{r}_{u}}}z$. Then, the following inequalities hold:

[TABLE]

The same result holds if $x=z\boxcircle_{\mathrm{r}}y$, with $\bar{r}_{\ell}=\bar{r}_{\ell}^{r}(S,x_{u},\boxcircle,z)$ and $\bar{r}_{u}=\bar{r}_{u}^{r}(S,x_{\ell},\boxcircle,z)$.

Proof

Here we only prove the claims for direct projections (namely, (9) and (10)), leaving those concerning indirect projections, which are analogous, to Appendix B.2.

First, we observe that, for each $x,y,z\in\mathbb{F}$, we have $[y\circ z]_{{\mathrm{n}}}=[y\circ z]_{{\mathord{\downarrow}}}$, or $[y\circ z]_{{\mathrm{n}}}=[y\circ z]_{{\mathord{\uparrow}}}$ or both. Then we prove that, for each $x,y,z\in\mathbb{F}$, we have $y\boxcircle_{\mathord{\downarrow}}z\preccurlyeq y\boxcircle_{\mathrm{n}}z\preccurlyeq y\boxcircle_{\mathord{\uparrow}}z$. We distinguish between the following cases, depending on $y\circ z$:

$y\circ z=+\infty\;\lor\;y\circ z=-\infty:$

in this case we have $y\boxcircle_{\mathord{\downarrow}}z=y\boxcircle_{\mathrm{n}}z=y\boxcircle_{\mathord{\uparrow}}z$ and thus $y\boxcircle_{\mathord{\downarrow}}z\preccurlyeq y\boxcircle_{\mathrm{n}}z\preccurlyeq y\boxcircle_{\mathord{\uparrow}}z$ holds.

$y\circ z\leq-{f_{\mathrm{min}}}\;\lor\;y\circ z\geq{f_{\mathrm{min}}}:$

in this case we have, by Proposition 1, $y\boxcircle_{\mathord{\downarrow}}z=[y\circ z]_{{\mathord{\downarrow}}}\leq[y\circ z]_{{\mathrm{n}}}=y\boxcircle_{\mathrm{n}}z\leq[y\circ z]_{{\mathord{\uparrow}}}=y\boxcircle_{\mathord{\uparrow}}z$; as $y\boxcircle_{\mathord{\downarrow}}z\neq 0$, $y\boxcircle_{\mathrm{n}}z\neq 0$ and $y\boxcircle_{\mathord{\uparrow}}z\neq 0$, the numerical order is reflected into the symbolic order to give $y\boxcircle_{\mathord{\downarrow}}z\preccurlyeq y\boxcircle_{\mathrm{n}}z\preccurlyeq y\boxcircle_{\mathord{\uparrow}}z$.

$-{f_{\mathrm{min}}}<y\circ z<0:$

in this case we have $y\boxcircle_{\mathord{\downarrow}}z=-{f_{\mathrm{min}}}\leq y\boxcircle_{\mathrm{n}}z\leq y\boxcircle_{\mathord{\uparrow}}z=-0$ by Definition 5; since either $[y\circ z]_{{\mathrm{n}}}=-{f_{\mathrm{min}}}$ or $[y\circ z]_{{\mathrm{n}}}=-0$, we have$[y\circ z]_{{\mathrm{n}}}\neq+0$, thus $y\boxcircle_{\mathord{\downarrow}}z\preccurlyeq y\boxcircle_{\mathrm{n}}z\preccurlyeq y\boxcircle_{\mathord{\uparrow}}z$.

$0<y\circ z<{f_{\mathrm{min}}}:$

in this case we have $y\boxcircle_{\mathord{\downarrow}}z=+0\leq y\boxcircle_{\mathrm{n}}z\leq y\boxcircle_{\mathord{\uparrow}}z={f_{\mathrm{min}}}$ by Definition 5; again, since either $[y\circ z]_{{\mathrm{n}}}=+0$ or $[y\circ z]_{{\mathrm{n}}}={f_{\mathrm{min}}}$ we know that $[y\circ z]_{{\mathrm{n}}}\neq-0$, and thus $y\boxcircle_{\mathord{\downarrow}}z\preccurlyeq y\boxcircle_{\mathrm{n}}z\preccurlyeq y\boxcircle_{\mathord{\uparrow}}z$.

$y\circ z=0:$

in this case, for multiplication and division the result is the same for all rounding modes, i.e., $+0$ or $-0$ depending on the sign of the arguments (IEEE-754-2008, , Section 6.3); for addition or subtraction we have $y\boxcircle_{\mathord{\downarrow}}z\neq-0$ while $y\boxcircle_{\mathrm{n}}z=y\boxcircle_{\mathord{\uparrow}}z=+0$; hence, also in this case, $y\boxcircle_{\mathord{\downarrow}}z\preccurlyeq y\boxcircle_{\mathrm{n}}z\preccurlyeq y\boxcircle_{\mathord{\uparrow}}z$ holds.

Note now that, by Definition 5, if $y\circ z>0$ then $y\boxcircle_{0}z=y\boxcircle_{\mathord{\downarrow}}z$ whereas, if $y\circ z>0$, then $y\boxcircle_{0}z=y\boxcircle_{\mathord{\uparrow}}z$. Therefore we can conclude that:

• •

if $y\circ z>0$, then $y\boxcircle_{\mathord{\downarrow}}z=y\boxcircle_{0}z\preccurlyeq y\boxcircle_{\mathrm{n}}z\preccurlyeq y\boxcircle_{\mathord{\uparrow}}z$ while,

• •

if $y\circ z<0$, then $y\boxcircle_{\mathord{\downarrow}}z\preccurlyeq y\boxcircle_{\mathrm{n}}z\preccurlyeq y\boxcircle_{0}z=y\boxcircle_{\mathord{\uparrow}}z;$ moreover,

• •

if $y\circ z=0$ and $\circ\notin\{\mathord{+},\mathord{-}\}$, then $y\boxcircle_{\mathord{\downarrow}}z=y\boxcircle_{\mathrm{n}}z=y\boxcircle_{0}z=y\boxcircle_{\mathord{\uparrow}}z$ while,

• •

if $y\circ z=0$ and $\circ\in\{\mathord{+},\mathord{-}\}$, then $y\boxcircle_{\mathord{\downarrow}}z\preccurlyeq y\boxcircle_{\mathrm{n}}z=y\boxcircle_{0}z=y\boxcircle_{\mathord{\uparrow}}z.$

In order to prove inequality (9), it is now sufficient to consider all possible sets $S\subseteq R$ and use the relations above.

For claim (10), observe that $r_{\ell}(S,y,\mathord{\boxcircle},z)\in S$ for any combination of rounding modes in $S$ except for one case: that is when $y\circ z>0$, and ${0}\in S$ but ${\mathord{\downarrow}}\notin S$. In this case, however, by Definition 5, $y\boxcircle_{\mathord{\downarrow}}z=y\boxcircle_{\mathrm{n}}z$. Similarly, $r_{u}(S,y,\mathord{\boxcircle},z)\in S$ except for the case when $y\circ z\leq 0$, ${0}\in S$ but ${\mathord{\uparrow}}\notin S$. First, assume that $y\circ z<0$: in this case, by Definition 5, $y\boxcircle_{\mathord{\uparrow}}z=y\boxcircle_{\mathrm{n}}z$. For the remaining case, that is $y\circ z=0$, we observe that for multiplication and division the result is the same for all rounding modes (IEEE-754-2008, , Section 6.3), while for addition or subtraction we have $y\boxcircle_{0}z=y\boxcircle_{\mathrm{n}}z=y\boxcircle_{\mathord{\uparrow}}z=+0$.

Thanks to Proposition 2 we need not be concerned with sets of rounding modes, as any such set $S\subseteq 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 $x=y\boxcircle_{S}z$, meaning “$x$ is obtained as the result of $y\boxcircle_{r}z$ for some $r\in S$.” Of course, $x=y\boxcircle_{S}z$ implies $x\succcurlyeq y\boxcircle_{S}z$ and $x\preccurlyeq y\boxcircle_{S}z$, which, by Proposition 2, imply $x\succcurlyeq y\boxcircle_{r_{\ell}}z$ and $x\preccurlyeq y\boxcircle_{r_{u}}z$, where $r_{\ell}\mathrel{\buildrel\mathrm{def}\over{=}}r_{\ell}(S,y,\mathord{\boxcircle},z)$ and $r_{u}\mathrel{\buildrel\mathrm{def}\over{=}}r_{u}(S,y,\mathord{\boxcircle},z)$. The results obtained by projection functions that only consider $r_{\ell}$ and $r_{u}$ are consequently valid for any $r\in S$.

3.2 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 $\mathord{\mathop{\nabla^{\uparrow}}\nolimits}\colon\mathbb{F}\rightarrowtail\overline{\mathbb{R}}$, $\mathord{\mathop{\nabla^{\downarrow}}\nolimits}\colon\mathbb{F}\rightarrowtail\overline{\mathbb{R}}$, $\mathord{\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits}\colon\mathbb{F}\rightarrowtail\overline{\mathbb{R}}$ and $\mathord{\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits}\colon\mathbb{F}\rightarrowtail\overline{\mathbb{R}}$ are defined as follows, for each $x\in\mathbb{F}$ that is not a NaN:

[TABLE]

An interesting observation is that the values of the functions introduced in Definition 9 are always representable in $\mathbb{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, $\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits$ and $\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits$ have been defined as twice the approximation error bounds: the absolute value of the bounds themselves, being ${f_{\mathrm{min}}}/2$, is not representable in $\mathbb{F}$ for each $x\in\mathbb{F}$ such that $|x|\leq{f^{\mathrm{nor}}_{\mathrm{min}}}$.

When the round-to-nearest rounding mode is in effect, Proposition 3 relates the bounds of a floating-point interval $[x_{\ell},x_{u}]$ with those of the corresponding interval of $\overline{\mathbb{R}}$ it represents.

Proposition 3

Let $x_{\ell},x_{u}\in\mathbb{F}\cap\mathbb{R}$. Then

[TABLE]

Proof (sketch)

To prove (15), we separately consider the two cases defined by (13).

If $x_{\ell}=-{f_{\mathrm{max}}}$, we prove that $x_{\ell}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x_{\ell})/2=-2^{{e_{\mathrm{max}}}}(2-2^{-p})$, while for all $x_{\ell}<x\leq x_{u}$ we have $x+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)/2=(x+\mathop{\mathrm{pred}}\nolimits(x))/2$. By monotonicity of ‘$\mathop{\mathrm{pred}}\nolimits$’, the minimum value of $(x+\mathop{\mathrm{pred}}\nolimits(x))/2$ occurs when $x=\mathop{\mathrm{succ}}\nolimits(-{f_{\mathrm{max}}})$, and

[TABLE]

which proves (15) in this case.

If, instead, $x_{\ell}>-{f_{\mathrm{max}}}$, applying monotonicity of ‘$\mathop{\mathrm{pred}}\nolimits$’ suffices.

The full proof is in Appendix B.2, together with the one of (16), which is symmetric.

3.3 Real Approximations of Floating-Point Constraints

In this section we show how inequalities of the form $x\succcurlyeq y\boxcircle_{r}z$ and $x\preccurlyeq y\boxcircle_{r}z$, with $r\in\{{\mathord{\downarrow}},{\mathord{\uparrow}},{\mathrm{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\in\mathbb{F}$ such that all the involved expressions do not evaluate to NaN, for each floating-point operation $\mathord{\boxcircle}\in\{\mathord{\boxplus},\mathord{\boxminus},\mathord{\boxdot},\mathord{\boxslash}\}$ and the corresponding extended real operation $\circ\in\{\mathord{+},\mathord{-},\mathord{\cdot},\mathord{/}\}$, where the entailed inequalities are to be interpreted over $\overline{\mathbb{R}}$:

[TABLE]

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.

3.4 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_{\mathbb{F}}$ of abstract syntax trees with leafs labelled by constants in $\mathbb{F}$ and internal nodes labeled with a symbol in $\{\mathord{+},\mathord{-},\mathord{\cdot},\mathord{/}\}$ 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: $\mathord{\llbracket\cdot\rrbracket_{{\mathord{\downarrow}}}}\colon E_{\mathbb{F}}\to\mathbb{F}$ and $\mathord{\llbracket\cdot\rrbracket_{{\mathord{\uparrow}}}}\colon E_{\mathbb{F}}\to\mathbb{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., KornerupLLM09 ).

Definition 10

(Evaluation functions.)

The two partial functions $\mathord{\llbracket\cdot\rrbracket_{{\mathord{\downarrow}}}}\colon E_{\mathbb{F}}\rightarrowtail\mathbb{F}$ and $\mathord{\llbracket\cdot\rrbracket_{{\mathord{\uparrow}}}}\colon E_{\mathbb{F}}\rightarrowtail\mathbb{F}$ are such that, for each $e\in\mathbb{F}$ that evaluates on $\overline{\mathbb{R}}$ to a nonzero value,

[TABLE]

Proposition 5

Let $x\in\mathbb{F}$ be a non-NaN floating point number and $e\in E_{\mathbb{F}}$ an expression that evaluates on $\overline{\mathbb{R}}$ to a nonzero value. The following implications hold:

[TABLE]

The implications of Proposition 5 can be derived from Definition 10 and Proposition 1. Their proof is postponed to Appendix B.2.

4 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.

4.1 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 $\mathord{\boxcircle}\in\{\mathord{\boxplus},\mathord{\boxminus},\mathord{\boxdot},\mathord{\boxslash}\}$ and $S\subseteq R$. Consider a constraint $x=y\boxcircle_{S}z$ with $x\in X=[x_{\ell},x_{u}]$, $y\in Y=[y_{\ell},y_{u}]$ and $z\in Z=[z_{\ell},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^{\prime}=[x^{\prime}_{\ell},x^{\prime}_{u}]\subseteq X$, such that

[TABLE]

Property (31) is known as the direct propagation correctness property.

Of course it is always possible to take $X^{\prime}=X$, but the objective of constraint propagation is to compute a “small”, possibly the smallest $X^{\prime}$ enjoying (31), compatibly with the available information. The smallest $X^{\prime}$ that satisfies (31) is called optimal and is such that

[TABLE]

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, $x=y\boxcircle_{S}z$, it means computing a possibly refined interval for $y$, $Y^{\prime}=[y^{\prime}_{\ell},y^{\prime}_{u}]\subseteq Y$, such that

[TABLE]

Property (33) is known as the inverse propagation correctness property. Again, taking $Y^{\prime}=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

[TABLE]

Property (34) is called the inverse propagation optimality property. Satisfying this last property can be very difficult.

4.2 The Boolean Domain for NaN

From now on, we will consider floating-point intervals with boundaries in $\mathbb{F}$. They allow for capturing the extended numbers in $\mathbb{F}$ only: NaNs (quiet NaNs and signaling NaNs) should be tracked separately. To this purpose, a Boolean domain $\mathcal{N}\mathrel{\buildrel\mathrm{def}\over{=}}\{\top,\bot\}$, where $\top$ stands for “may be NaN” and $\bot$ means “cannot be NaN”, can be used and coupled with the arithmetic filtering algorithms.

Let be $x=y\mathbin{\boxcircle}z$ an arithmetic constraint over floating-point numbers, and $(X,\mathrm{NaN}_{x})$, $(Y,\mathrm{NaN}_{y})$ and $(Z,\mathrm{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^{\prime},\mathrm{NaN}_{x}^{\prime})$, $(Y^{\prime},\mathrm{NaN}_{y}^{\prime})$ and $(Z^{\prime},\mathrm{NaN}_{z}^{\prime})$ 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 (IEEE-754-2008, , 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., $+\infty\boxplus-\infty$ or $+\infty\boxminus+\infty$), and also $\pm\infty\boxdot 0=\mathrm{NaN}$, $0\boxdot\pm\infty=\mathrm{NaN}$, $0\boxslash 0=\mathrm{NaN}$, and $\pm\infty\boxslash\pm\infty=\mathrm{NaN}$.

Thus, direct projections are such that if $\mathrm{NaN}_{y}=\top$ or $\mathrm{NaN}_{z}=\top$, then also $\mathrm{NaN}_{x}^{\prime}=\top$; indirect projections yield $\mathrm{NaN}_{y}^{\prime}=\mathrm{NaN}_{z}^{\prime}=\top$ if $\mathrm{NaN}_{x}=\top$. Moreover, e.g., if $\boxcircle=\boxplus$, then the direct projection yields $\mathrm{NaN}_{x}^{\prime}=\top$ also if $\pm\infty\in Y$ and $\mp\infty\in Z$, and the indirect one allows for $\pm\infty$ in $Y^{\prime}$ and $\mp\infty$ in $Z^{\prime}$ only if $\mathrm{NaN}_{x}=\top$, and so on for the other operators.

4.3 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 ($\boxplus$) and division ($\boxslash$) 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.

4.3.1 Addition

Here we deal with constraints of the form $x=y\boxplus_{S}z$ with $S\subseteq R$. Let $X=[x_{\ell},x_{u}]$, $Y=[y_{\ell},y_{u}]$ and $Z=[z_{\ell},z_{u}]$.

Thanks to Proposition 2, any set of rounding modes $S\subseteq 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$.

Direct Propagation.

For direct propagation, i.e., the process that infers a new interval for $x$ starting from the interval for $y$ and $z$, we propose Algorithm 1 and functions $\mathrm{da}_{\ell}$ and $\mathrm{da}_{u}$, as defined in Figure 2. Functions $\mathrm{da}_{\ell}$ and $\mathrm{da}_{u}$ yield new bounds for interval $X$. In particular, function $\mathrm{da}_{\ell}$ gives the new lower bound, while function $\mathrm{da}_{u}$ provides the new upper bound of the interval. Functions $\mathrm{da}_{\ell}$ and $\mathrm{da}_{u}$ handle all rounding modes and, in order to be as precise as possible, they distinguish between several cases, depending on the values of the bounds of intervals $Y$ and $Z$. These cases are infinities ($-\infty$ and $+\infty$), zeroes ($-0$ and $+0$), negative values ($\mathbb{R}_{-}$) and positive values ($\mathbb{R}_{+}$).

It can be proved that Algorithm 1 computes a correct and optimal direct projection, as stated by its postconditions.

Theorem 4.1

Algorithm 1* satisfies its contract.*

Proof

Given the constraint $x=y\boxplus_{S}z$ with $x\in X=[x_{\ell},x_{u}]$, $y\in Y=[y_{\ell},y_{u}]$ and $z\in Z=[z_{\ell},z_{u}]$, Algorithm 1 sets $X^{\prime}=[x^{\prime}_{\ell},x^{\prime}_{u}]\cap X$; hence, we have $X^{\prime}\subseteq X$. Moreover, by Proposition 2, for each $y\in Y$, $z\in Z$ and $r\in S$, we have $y\boxplus_{r_{\ell}}z\preccurlyeq y\boxplus_{r}z\preccurlyeq y\boxplus_{r_{u}}z$, and because $a\preccurlyeq b$ implies $a\leq b$ for any $a,b\in\mathbb{F}$ according to Definition 3, we know that $y\boxplus_{r_{\ell}}z\leq y\boxplus_{r}z\leq y\boxplus_{r_{u}}z$. Thus, by monotonicity of $\boxplus$, we have $y_{\ell}\boxplus_{r_{\ell}}z_{\ell}\leq y\boxplus_{r_{\ell}}z\leq y\boxplus_{r}z\leq y\boxplus_{r_{u}}z\leq y_{u}\boxplus_{r_{u}}z_{u}$. Therefore, we can focus on finding a lower bound for $y_{\ell}\boxplus_{r_{\ell}}z_{\ell}$ and an upper bound for $y_{u}\boxplus_{r_{u}}z_{u}$.

Such bounds are given by the functions $\mathrm{da}_{\ell}$ and $\mathrm{da}_{u}$ of Figure 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 IEEE-754-2008 ; just two cases need further explanation. Concerning the entry of $\mathrm{da}_{\ell}$ in which $y_{\ell}=-\infty$ and $z_{\ell}=+\infty$, note that $z_{\ell}=+\infty$ implies $z_{u}=+\infty$. Then for any $y>y_{\ell}=-\infty$, $y\boxplus+\infty=+\infty$. On the other hand, by the IEEE 754 Standard IEEE-754-2008 , $-\infty\boxplus+\infty$ is an invalid operation. For the symmetric case, i.e., the entry of $\mathrm{da}_{u}$ in which $y_{u}=-\infty$ and $z_{u}=+\infty$, we can reason dually.

We are now left to prove that $\forall X^{\prime\prime}\subset X^{\prime}\mathrel{:}\exists r\in S,y\in Y,z\in Z\mathrel{:}y\boxplus_{r}z\not\in X^{\prime\prime}$. Let us focus on the lower bound $x^{\prime}_{\ell}$, proving that there always exists a $r\in S$ such that $y_{\ell}\boxplus_{r}z_{\ell}=x^{\prime}_{\ell}$.

First, consider the cases in which $y_{\ell}\not\in(\mathbb{R}_{-}\cup\mathbb{R}_{+})$ or $z_{\ell}\not\in(\mathbb{R}_{-}\cup\mathbb{R}_{+})$. In these cases, a case analysis proves that $\mathrm{da}_{\ell}(y_{\ell},z_{\ell},r_{\ell})$ is equal to $y_{\ell}\boxplus_{r_{\ell}}z_{\ell}$. Indeed, if either of the operands (say $y_{\ell}$) is $-\infty$ and the other one (say $z_{\ell}$) is not $+\infty$, then according to the IEEE 754 Standard we have $y_{\ell}\boxplus_{r}z_{\ell}=-\infty$ for any $r\in R$. Symmetrically, $y_{\ell}\boxplus_{r}z_{\ell}=+\infty$ if one operand is $+\infty$ and the other one is not $-\infty$. If, w.l.o.g., $y_{\ell}=+\infty$ and $z_{\ell}=-\infty$, the set $X^{\prime}$ is non-empty only if $z_{u}\neq-\infty$, and $y_{\ell}\boxplus_{r}z_{u}=+\infty$ for any $r\in R$.

For the cases in which $y_{\ell}\in(\mathbb{R}_{-}\cup\mathbb{R}_{+})$ and $z_{\ell}\in(\mathbb{R}_{-}\cup\mathbb{R}_{+})$ we have $x^{\prime}_{\ell}=y_{\ell}\boxplus_{r_{\ell}}z_{\ell}$, by definition of $\mathrm{da}_{\ell}$ of Figure 2. Remember that, by Proposition 2, there exists $r\in S$ such that $y_{\ell}\boxplus_{\mathrm{r_{\ell}}}z_{\ell}=y_{\ell}\boxplus_{\mathrm{r}}z_{\ell}$. Since $y_{\ell}\in Y$ and $z_{\ell}\in Z$, we can conclude that for any $X^{\prime\prime}\subseteq X^{\prime}$, $x^{\prime}_{\ell}\not\in X^{\prime\prime}$ implies $y_{\ell}\boxplus_{\mathrm{r}}z_{\ell}\not\in X^{\prime\prime}$.

An analogous argument allows us to conclude that there exists an $r\in S$ for which the following holds: for any $X^{\prime\prime}\subseteq X^{\prime}$, $x^{\prime}_{u}\not\in X^{\prime\prime}$ implies $y_{u}\boxplus_{\mathrm{r}}z_{u}\not\in X^{\prime\prime}$.

The following example will better illustrate how the tables in Figure 2 should be used to compute functions $\mathrm{da}_{\ell}$ and $\mathrm{da}_{u}$. All examples in this section refer to the IEEE 754 binary single precision format.

Example 1

Assume $Y=[+0,5]$, $Z=[-0,8]$, and that the selected rounding mode is $r_{\ell}=r_{u}={\mathord{\downarrow}}$. In order to compute the lower bound $x_{\ell}^{\prime}$ of $X^{\prime}$, the new interval for $x$, function $\mathrm{da}_{\ell}(+0,-0,{\mathord{\downarrow}})$ is called. These arguments fall in case $a_{1}$, which yields $-0$ with rounding mode ${\mathord{\downarrow}}$. Indeed, when the rounding mode is ${\mathord{\downarrow}}$, 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}^{\prime}$, the algorithm calls $\mathrm{da}_{u}(5,8,{\mathord{\downarrow}})$. This falls in the case in which both operands are positive numbers ($y_{u},z_{u}\in\mathbb{R}_{+}$), and therefore $x_{u}=y_{u}\boxplus_{r_{u}}z_{u}=13$. In conclusion, the new interval for $x$ is $X^{\prime}=[-0,13]$.

If any other rounding mode was selected (say, $r_{\ell}=r_{u}={\mathrm{n}}$), the new interval computed by the projection would have been $X^{\prime\prime}=[+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 $\mathrm{ia}_{\ell}$ in Figure 3 and $\mathrm{ia}_{u}$ in Figure 4, where $\equiv$ 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_{\ell}$ and $z_{u}$; analogously, the values of $x$ and $z$ that maximize $y$ are $x_{u}$ and $z_{\ell}$.

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 4.2) that when $x$ is subnormal $\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)$ and $\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)$ are negligible. This allows us to define tight bounds in this case. On the contrary, when the terms $\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x_{\ell})$ and $\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x_{u})$ are non negligible, we need to approximate the values of expressions $e_{\ell}$ and $e_{u}$. This can always be done with reasonable efficiency KornerupLLM09 , but we leave this as an implementation choice, thus accounting for the case when the computation is exact ($\llbracket e_{\ell}\rrbracket_{{\mathord{\uparrow}}}=[e_{\ell}]_{{\mathord{\uparrow}}}$ and $\llbracket e_{u}\rrbracket_{{\mathord{\downarrow}}}=[e_{u}]_{{\mathord{\downarrow}}}$) as well as when it is not ($\llbracket e_{\ell}\rrbracket_{{\mathord{\uparrow}}}>[e_{\ell}]_{{\mathord{\uparrow}}}$ and $\llbracket e_{u}\rrbracket_{{\mathord{\downarrow}}}<[e_{u}]_{{\mathord{\downarrow}}}$).

The next result assures us that our algorithm computes a correct inverse projection, as claimed by its postcondition.

Theorem 4.2

Algorithm 2* satisfies its contract.*

Proof

Given the constraint $x=y\boxplus_{S}z$ with $x\in X=[x_{\ell},x_{u}]$, $y\in Y=[y_{\ell},y_{u}]$ and $z\in Z=[z_{\ell},z_{u}]$, Algorithm 2 computes a new and refined domain $Y^{\prime}$ for variable $y$.

First, observe that the newly computed interval $[y^{\prime}_{\ell},y^{\prime}_{u}]$ is either intersected with the old domain $Y$, so that $Y^{\prime}=[y^{\prime}_{\ell},y^{\prime}_{u}]\cap Y$, or set to $Y^{\prime}=\varnothing$. Hence, $Y^{\prime}\subseteq Y$ holds.

Proposition 2 and the monotonicity of $\boxplus$ allow us to find a lower bound for $y$ by exploiting the constraint $y\boxplus_{\bar{r}_{\ell}}z_{u}=x_{\ell}$, and an upper bound for $y$ by exploiting the constraint $y\boxplus_{\bar{r}_{u}}z_{\ell}=x_{u}$. We will now prove that the case analyses of functions $\mathrm{ia}_{\ell}$, described in Figure 3, and $\mathrm{ia}_{u}$, described in Figure 4, express such bounds correctly.

Concerning the operand combinations in which $\mathrm{ia}_{\ell}$ takes the value described by the case analysis $a_{4}$, remember that, by the IEEE 754 Standard IEEE-754-2008 , whenever the sum of two operands with opposite sign is zero, the result of that sum is $+0$ in all rounding-direction attributes except roundTowardNegative: in that case the result is $-0$. Then, since $z_{u}\boxplus_{{\mathord{\downarrow}}}(-z_{u})=-0$, when $\bar{r}_{\ell}={\mathord{\downarrow}}$, $y_{\ell}$ can safely be set to $\mathop{\mathrm{succ}}\nolimits(-z_{u})$.

As for the case in which $\mathrm{ia}_{\ell}$ takes one of the values determined by $a_{5}$, the IEEE 754 Standard IEEE-754-2008 asserts that $+0\boxplus_{{\mathord{\downarrow}}}+0=+0$, while $-0\boxplus_{{\mathord{\downarrow}}}+0=-0$: the correct lower bound for $y$ is $y^{\prime}_{\ell}=+0$, in this case. As we already pointed out, for any other rounding-direction attribute $+0\boxplus-0=+0$ holds, which allows us to include $-0$ in the new domain.

Concerning cases of $\mathrm{ia}_{\ell}$ that give the result described by the case analysis $a_{6}$, we clearly must have $y=+\infty$ if $\bar{r}_{\ell}={\mathord{\downarrow}}$; if $\bar{r}_{\ell}={\mathord{\uparrow}}$, it should be $y+z_{u}>{f_{\mathrm{max}}}$ and thus $y>{f_{\mathrm{max}}}-z_{u}$ and, by (28) of Proposition 5, $y\succcurlyeq\mathop{\mathrm{succ}}\nolimits({f_{\mathrm{max}}}\boxminus_{\mathord{\downarrow}}z_{u})$. If $\bar{r}_{\ell}={\mathrm{n}}$, there are two cases:

$z_{u}<\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2.$

In this case, $y$ must be greater than ${f_{\mathrm{max}}}$, since ${f_{\mathrm{max}}}+z_{u}<{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2$ implies that ${f_{\mathrm{max}}}\boxplus_{{\mathrm{n}}}z_{u}={f_{\mathrm{max}}}<+\infty$. Note that in this case $\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2\boxminus_{\mathord{\uparrow}}z_{u}\geq{f_{\mathrm{min}}}$, hence {f_{\mathrm{max}}}\boxplus_{\mathord{\uparrow}}\bigl{(}\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2\boxminus_{\mathord{\uparrow}}z_{u}\bigr{)}=+\infty.

$z_{u}\geq\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2.$

Since $\mathop{\mathrm{odd}}\nolimits({f_{\mathrm{max}}})$, for $x_{\ell}=+\infty$ we need $y$ to be greater than or equal to ${f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2-z_{u}$. Note that $y\geq{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2-z_{u}$ together with

[TABLE]

allows us to apply (29) of Proposition 5, concluding y\succcurlyeq{f_{\mathrm{max}}}\boxplus_{\mathord{\uparrow}}\bigl{(}\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2\boxminus_{\mathord{\uparrow}}z_{u}\bigr{)}.

Equality (35) holds because either the application of ‘$\mathord{\boxminus_{\mathord{\uparrow}}}$’ is exact or the application of ‘$\mathord{\boxplus_{\mathord{\uparrow}}}$’ is exact. In fact, since $z_{u}=m\cdot 2^{e}\geq\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2=2^{{e_{\mathrm{max}}}-p}$, for some $1\leq m<2$, there are two cases: either $e={e_{\mathrm{max}}}$ or ${e_{\mathrm{max}}}-p\leq e<{e_{\mathrm{max}}}$.

Suppose first that $e={e_{\mathrm{max}}}$: we have

[TABLE]

and thus

[TABLE]

Since if $e={e_{\mathrm{max}}}$ the application of ‘$\mathord{\boxminus_{\mathord{\uparrow}}}$’ is not exact, we prove that the application of ‘$\mathord{\boxplus_{\mathord{\uparrow}}}$’ is exact. Hence, if $m>1$, we prove that

[TABLE]

where k\mathrel{\buildrel\mathrm{def}\over{=}}-\bigl{\lfloor}\log_{2}(2-m)\bigr{\rfloor}. It is worth noting that $2^{k}(2-m)$ can be represented by a normalized mantissa; moreover, since $1\leq k\leq p-1$, ${e_{\mathrm{min}}}\leq{e_{\mathrm{max}}}-k\leq{e_{\mathrm{max}}}$, hence, ${f_{\mathrm{max}}}+(\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2\boxminus_{\mathord{\uparrow}}z_{u})\in\mathbb{F}$. If, instead, $m=1$,

[TABLE]

and, also in this case, ${f_{\mathrm{max}}}+(\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2\boxminus_{\mathord{\uparrow}}z_{u})\in\mathbb{F}$.

Suppose now that ${e_{\mathrm{max}}}-p\leq e<{e_{\mathrm{max}}}$ and let $h\mathrel{\buildrel\mathrm{def}\over{=}}e-{e_{\mathrm{max}}}+p$ so that $0\leq h\leq p-1$. In this case we show that the application of ‘$\mathord{\boxminus_{\mathord{\uparrow}}}$’ is exact. Indeed, we have

[TABLE]

If $e={e_{\mathrm{max}}}-p$ and $m=1$, then $h=0$, $m-2^{-h}=0$ and thus $\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2-z_{u}=0$. Otherwise, let k\mathrel{\buildrel\mathrm{def}\over{=}}-\bigl{\lfloor}\log_{2}(m-2^{-h})\bigr{\rfloor}. We have

[TABLE]

which is an element of $\mathbb{F}$.

Dual arguments w.r.t. the ones used to justify cases of $\mathrm{ia}_{\ell}$ that give the result described by $a_{4}$, $a_{6}$ and $a_{5}$ can be used to justify the cases of $\mathrm{ia}_{u}$ described by $a_{9}$, $a_{10}$ and $a_{7}$.

We now tackle the entries of $\mathrm{ia}_{\ell}$ described by $a_{3}$, and those of $\mathrm{ia}_{u}$ described by $a_{8}$. Exploiting $x\preccurlyeq y\boxplus z$ and $x\succcurlyeq y\boxplus z$, by Proposition 4, we have

[TABLE]

The same case analysis gives us

[TABLE]

We can now exploit the fact that $x\in[x_{\ell},x_{u}]$ and $z\in[z_{\ell},z_{u}]$ with $x_{\ell},x_{u},z_{\ell},z_{u}\in\mathbb{F}$ to obtain, using Proposition 3 and the monotonicity of ‘$\mathop{\mathrm{pred}}\nolimits$’ and ‘$\mathop{\mathrm{succ}}\nolimits$’:

[TABLE]

We can now exploit Proposition 5 and obtain

[TABLE]

In fact, if $x_{\ell}=z_{u}$, then, according to IEEE 754 (IEEE-754-2008, , Section 6.3), for each non-NaN, nonzero and finite $w\in\mathbb{F}$, $-0$ is the least value for $y$ that satisfies $w=y\boxplus_{\mathord{\downarrow}}w$. If $x_{\ell}\neq z_{u}$, then case (29) of Proposition 5 applies and we have $y\succcurlyeq x_{\ell}\boxminus_{\mathord{\uparrow}}z_{u}$. Suppose now that $\mathop{\mathrm{pred}}\nolimits(x_{\ell})=z_{u}$, then $\mathop{\mathrm{pred}}\nolimits(x_{\ell})\boxminus_{\mathord{\downarrow}}z_{u}\equiv 0$ and \mathop{\mathrm{succ}}\nolimits\bigl{(}\mathop{\mathrm{pred}}\nolimits(x_{\ell})\boxminus_{\mathord{\downarrow}}z_{u}\bigr{)}={f_{\mathrm{min}}}, coherently with the fact that, for each non-NaN, nonzero and finite $w\in\mathbb{F}$, ${f_{\mathrm{min}}}$ is the least value for $y$ that satisfies $w=y\boxplus_{\mathord{\uparrow}}\mathop{\mathrm{pred}}\nolimits(w)$. If $\mathop{\mathrm{pred}}\nolimits(x_{\ell})\neq z_{u}$, then case (26) of Proposition 5 applies and we have y\succcurlyeq\mathop{\mathrm{succ}}\nolimits\bigl{(}\mathop{\mathrm{pred}}\nolimits(x_{\ell})\boxminus_{\mathord{\downarrow}}z_{u}\bigr{)}. A symmetric argument justifies (39).

For the remaining cases, we first show that when $\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)={f_{\mathrm{min}}}$,

[TABLE]

The previous equality has the following main consequences: we can perform the computation in $\mathbb{F}$, that is, we do not need to compute $\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)/2$ and, since $[x_{u}-z_{\ell}]_{{\mathord{\downarrow}}}=\llbracket x_{u}-z_{\ell}\rrbracket_{{\mathord{\downarrow}}}$, we can apply (30) of Proposition 5, obtaining a tight bound for $y_{u}^{\prime}$.

Let us prove (40). Suppose $\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x_{u})={f_{\mathrm{min}}}$, and assume $x_{u}\neq z_{\ell}$. There are two cases:

$[x_{u}-z_{\ell}]_{{\mathord{\downarrow}}}=x_{u}-z_{\ell}:$

then we have $y\leq[x_{u}-z_{\ell}]_{{\mathord{\downarrow}}}=x_{u}-z_{\ell}$ since the addition of $\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x_{u})/2={f_{\mathrm{min}}}/2$ is insufficient to reach $\mathop{\mathrm{succ}}\nolimits(x_{u}-z_{\ell})$, whose distance from $x_{u}-z_{\ell}$ is at least ${f_{\mathrm{min}}}$.

$[x_{u}-z_{\ell}]_{{\mathord{\downarrow}}}<x_{u}-z_{\ell}<[x_{u}-z_{\ell}]_{{\mathord{\uparrow}}}:$

since by Definition 2 every finite floating-point number is an integral multiple of ${f_{\mathrm{min}}}$, so are $x_{u}-z_{\ell}$ and $[x_{u}-z_{\ell}]_{{\mathord{\uparrow}}}$. Therefore, again, $y\leq[x_{u}-z_{\ell}]_{{\mathord{\downarrow}}}$, since the addition of $\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x_{u})/2={f_{\mathrm{min}}}/2$ $x_{u}-z_{\ell}$ is insufficient to reach $[x_{u}-z_{\ell}]_{{\mathord{\uparrow}}}$, whose distance from $x_{u}-z_{\ell}$ is at least ${f_{\mathrm{min}}}$.

In the case where $x_{u}=z_{\ell}$ we have \bigl{[}x_{u}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x_{u})/2-z_{\ell}\bigr{]}_{{\mathord{\downarrow}}}=[0+{f_{\mathrm{min}}}/2]_{{\mathord{\downarrow}}}=+0, hence (40) holds. As we have already pointed out, this allows us to apply (30) of Proposition 5 to the case $\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x_{u})={f_{\mathrm{min}}}$, obtaining the bound $y\preccurlyeq[x_{u}-z_{\ell}]_{{\mathord{\downarrow}}}$.

Similar arguments can be applied to $\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x_{\ell})$ whenever $\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x_{\ell})=-{f_{\mathrm{min}}}$ to prove that \bigl{[}x_{\ell}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x_{\ell})/2-z_{u}\bigr{]}_{{\mathord{\uparrow}}}=[x_{\ell}-z_{u}]_{{\mathord{\uparrow}}}. Then, by (29) of Proposition 5, we obtain $y\succcurlyeq[x_{\ell}-z_{u}]_{{\mathord{\uparrow}}}$.

When the terms $\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x_{\ell})$ and $\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x_{u})$ are non-negligible, we need to approximate the values of the expressions $e_{\ell}\mathrel{\buildrel\mathrm{def}\over{=}}x_{\ell}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x_{\ell})/2-z_{u}$ and $e_{u}\mathrel{\buildrel\mathrm{def}\over{=}}x_{u}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x_{u})/2-z_{\ell}$. Hence, we have the cases $\llbracket e_{\ell}\rrbracket_{{\mathord{\uparrow}}}=[e_{\ell}]_{{\mathord{\uparrow}}}$ and $\llbracket e_{u}\rrbracket_{{\mathord{\downarrow}}}=[e_{u}]_{{\mathord{\downarrow}}}$ as well as $\llbracket e_{\ell}\rrbracket_{{\mathord{\uparrow}}}>[e_{\ell}]_{{\mathord{\uparrow}}}$ and $\llbracket e_{u}\rrbracket_{{\mathord{\downarrow}}}<[e_{u}]_{{\mathord{\downarrow}}}$. Thus, when $\llbracket e_{u}\rrbracket_{{\mathord{\downarrow}}}<[e_{u}]_{{\mathord{\downarrow}}}$ by (37) and (27) of Proposition 5 we obtain $y\preccurlyeq\llbracket e_{u}\rrbracket_{{\mathord{\uparrow}}}$, while, when $\llbracket e_{\ell}\rrbracket_{{\mathord{\downarrow}}}>[e_{\ell}]_{{\mathord{\downarrow}}}$ by (37) and (25) of Proposition 5 we obtain $y\succcurlyeq\llbracket e_{\ell}\rrbracket_{{\mathord{\downarrow}}}$. Finally, when $\mathop{\mathrm{odd}}\nolimits(x_{u})$, by (37) and (28) of Proposition 5, we obtain y\preccurlyeq\mathop{\mathrm{pred}}\nolimits\bigl{(}\llbracket e_{u}\rrbracket_{{\mathord{\uparrow}}}\bigr{)}. Dually, when $\mathop{\mathrm{odd}}\nolimits(x_{\ell})$ by (36) and (26) of Proposition 5, we obtain y\succcurlyeq\mathop{\mathrm{succ}}\nolimits\bigl{(}\llbracket e_{\ell}\rrbracket_{{\mathord{\downarrow}}}\bigr{)}.

Thus, for the case $\bar{r}_{\ell}={\mathrm{n}}$ we have

[TABLE]

whereas, for the case where $\bar{r}_{u}={\mathrm{n}}$, we have

[TABLE]

Example 2

Let $X=[+0,+\infty]$ and $Z=[-\infty,+\infty]$. Regardless of the rounding mode, the calls to functions $\mathrm{ia}_{\ell}(+0,+\infty,\bar{r}_{\ell})$ and $\mathrm{ia}_{u}(+\infty,-\infty,\bar{r}_{u})$ yield $Y^{\prime}=[-{f_{\mathrm{max}}},+\infty]$. Note that $-{f_{\mathrm{max}}}$ is the lowest value that variable $y$ could take, since there is no value for $z\in Z$ that summed with $-\infty$ gives a value in $X$. Indeed, if we take $z=+{f_{\mathrm{max}}}$, then we have $-{f_{\mathrm{max}}}\boxplus_{r}+{f_{\mathrm{max}}}=+0\in X$ for any $r\in R$. On the other hand, $+\infty$ is clearly the highest value $y$ could take, since $+\infty\boxplus_{r}z=+\infty\in X$ for any value of $z\in Z\setminus\{-\infty\}$. In this case, our projections yield a more refined result than the competing tool FPSE BotellaGM06 , which computes the wider interval $Y^{\prime}=[-\infty,+\infty]$.

Example 3

Consider also $X=[1.0,2.0]$ and $Z=[-{1.0}\times{2^{30}},{1.0}\times{2^{30}}]$ and $S=\{{\mathrm{n}}\}.$ With our inverse projection we obtain $Y=[-{1.1\cdots 1}\times{2^{29}},{1.0}\times{2^{30}}]$ which is correct but not optimal. For example, pick $y={1.0}\times{2^{30}}$: for $z=-{1.0}\times{2^{30}}$ we have $y\boxplus_{S}z=0$ and $y\boxplus_{S}z^{+}=64$. By monotonicity of $\boxplus_{S}$, for no $z\in[-{1.0}\times{2^{30}},{1.0}\times{2^{30}}]$ we can have $y\boxplus_{S}z\in[1.0,2.0]$.

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, MarreM10 identified a property of the representation of floating-point numbers and proposed to exploit it in filtering algorithms for addition and subtraction constraints. In BagnaraCGG13ICST ; BagnaraCGG16IJOC 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 BagnaraCGG16IJOC 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.

Example 4

Assume, again, $X=[1.0,2.0]$ and $Z=[-{1.0}\times{2^{30}},{1.0}\times{2^{30}}]$ and $S=\{{\mathrm{n}}\}$. By applying maximum ULP filtering MarreM10 ; BagnaraCGG16IJOC , we obtain the much tighter intervals $Y,Z=[-{1.1\cdots 1}\times{2^{24}},{1.0}\times{2^{25}}]$. These are actually optimal as $-{1.1\cdots 1}\times{2^{24}}\boxplus_{S}{1.0}\times{2^{25}}={1.0}\times{2^{25}}\boxplus_{S}-{1.1\cdots 1}\times{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\cdots 1}\times{2^{24}},{1.0}\times{2^{25}}]$; in contrast, our inverse projection exploits the available information on $Z$ to obtain $Y=[-4,1.0\cdots 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^{\prime}$ of $Z$ can be defined analogously.

4.3.2 Division

In this section we deal with constraints of the form $x=y\boxslash_{S}z$ with $S\subseteq R$.

Direct Propagation.

For direct propagation, interval $Z$ is partitioned into the sign-homogeneous intervals $Z_{-}\mathrel{\buildrel\mathrm{def}\over{=}}Z\cap[-\infty,-0]$ and $Z_{+}\mathrel{\buildrel\mathrm{def}\over{=}}Z\cap[+0,+\infty]$. 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^{-}_{\ell},x^{-}_{u}]$, and $Y$ and $W=Z_{+}$, to obtain $[x^{+}_{\ell},x^{+}_{u}]$. The appropriate bounds for interval propagation are chosen by function $\tau$ of Figure 5. Note that the sign of $z$ is, by construction, constant over interval $W$. The selected values are then taken as arguments by functions $\mathrm{dd}_{\ell}$ and $\mathrm{dd}_{u}$ of Figure 6, which return the correct bounds for the aforementioned new intervals for $X$. The intervals $X\cap[x^{-}_{\ell},x^{-}_{u}]$ and $X\cap[x^{+}_{\ell},x^{+}_{u}]$ are eventually joined using convex union, denoted by $\biguplus$, to obtain the refining interval $X^{\prime}$.

It can be proved that Algorithm 3 computes a correct and optimal direct projection, as ensured by its postconditions.

Theorem 4.3

Algorithm 3* satisfies its contract.*

Example 5

Consider $Y=[-0,42]$, $Z=[-3,6]$ and any value of $S$. First, $Z$ is split into $Z_{-}=[-3,-0]$ and $Z_{+}=[+0,6]$. For the negative interval, the third case of $\tau(-0,42,-3,-0)$ applies, yielding $(y_{L},y_{U},w_{L},w_{U})=(42,-0,-0,-0)$. Then, the projection functions are invoked, and we have $\mathrm{dd}_{\ell}(42,-0,r_{\ell})=-\infty$ and $\mathrm{dd}_{u}(-0,-0,r_{u})=+0$, i.e., $[x^{-}_{\ell},x^{-}_{u}]=[-\infty,+0]$. For the positive part, we have $\tau(-0,42,+0,6)=(-0,42,+0,+0)$ (sixth case). From the projections we obtain $\mathrm{dd}_{\ell}(-0,+0,r_{\ell})=-0$ and $\mathrm{dd}_{u}(42,+0,r_{u})=+\infty$, and $[x^{+}_{\ell},x^{+}_{u}]=[-0,+\infty]$. Finally, $X^{\prime}=[x^{-}_{\ell},x^{-}_{u}]\biguplus[x^{+}_{\ell},x^{+}_{u}]=[-\infty,+\infty]$.

Inverse Propagation (First Projection).

The inverse projections of division must be handled separately for each operand. The projection on $y$ is the first inverse projection. This case requires, as explained for Algorithm 3, to split $Z$ into the sign-homogeneous intervals $Z_{-}\mathrel{\buildrel\mathrm{def}\over{=}}Z\cap[-\infty,-0]$ and $Z_{+}\mathrel{\buildrel\mathrm{def}\over{=}}Z\cap[+0,+\infty]$. Then, in order to select the extrema that determine the appropriate lower and upper bound for $y$, function $\sigma$ of Figure 7 is applied.

Example 6

Suppose $X=[-42,+0]$, $Z=[-{1.0}\times{2^{100}},-0]$ and $S=\{{\mathrm{n}}\}$. In this case, $Z_{-}=Z$, and $Z_{+}=\varnothing$. We obtain $\sigma(-{1.0}\times{2^{100}},-0,-42,+0)=(-{1.0}\times{2^{100}},-{1.0}\times{2^{100}},+0,-42)$ from the sixth case of $\sigma$. Then, $\mathrm{id}^{f}_{\ell}(+0,-{1.0}\times{2^{100}},{\mathrm{n}})=({f_{\mathrm{min}}}\boxdot_{\mathord{\uparrow}}(-{1.0}\times{2^{100}}))/2=-{1.0}\times{2^{-50}}$, because the lowest value of $y_{\ell}$ is obtained when a division by $-{1.0}\times{2^{100}}$ underflows. Moreover, $\mathrm{id}^{f}_{u}(-42,-{1.0}\times{2^{100}},{\mathrm{n}})={1.0101}\times{2^{105}}$. Therefore, the projected interval is $Y^{\prime}=[-{1.0}\times{2^{-50}},{1.0101}\times{2^{105}}]$.

The following result assures us that Algorithm 4 computes a correct first inverse projection, as ensured by its postcondition.

Theorem 4.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 MarreM10 and extended to multiplication and division in BagnaraCGG16IJOC . 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 BagnaraCGG16IJOC .

Example 7

Consider the IEEE 754 single-precision constraint $x=y\boxslash_{S}z$ with initial intervals $X=[-{1.0}\times{2^{-110}},-{1.0}\times{2^{-121}}]$ and $Y=Z=[-\infty,+\infty]$. When $S=\{{\mathrm{n}}\}$, filtering by maximum ULP results in the possible refinement $Y^{\prime}=[-{1.1\cdots 1}\times{2^{17}},{1.1\cdots 1}\times{2^{17}}]$, while Algorithm 4 would return the less precise $Y^{\prime}=[-{f_{\mathrm{max}}},{f_{\mathrm{max}}}]$, with any rounding mode.

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_{-}\mathrel{\buildrel\mathrm{def}\over{=}}X\cap[-\infty,-0]$ and $X_{+}\mathrel{\buildrel\mathrm{def}\over{=}}X\cap[+0,+\infty]$ 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^{-}_{\ell},z^{-}_{u}]$; intervals $X_{+}$ and $Y$ to obtain $[z^{+}_{\ell},z^{+}_{u}]$. The new bounds for $z$ are computed by functions $\mathrm{id}^{s}_{\ell}$ of Figure 10 and $\mathrm{id}^{s}_{u}$ of Figure 11, after the appropriate interval extrema of $Y$ and $V=X_{-}$ (or $V=X_{+}$) have been selected by function $\tau$. The intervals $Z\cap[z^{-}_{\ell},z^{-}_{u}]$ and $Z\cap[z^{+}_{\ell},z^{+}_{u}]$ will be then joined with convex union to obtain $Z^{\prime}$.

Our algorithm computes a correct second inverse projection.

Theorem 4.5

Algorithm 5* satisfies its contract.*

Example 8

Consider $X=[6,+\infty]$, $Y=[+0,42]$ and $S=\{{\mathrm{n}}\}$. In this case, we only have $X_{+}=X$, and $X_{-}=\varnothing$. With this input, $\tau(+0,42,6,+\infty)=(+0,42,+\infty,6)$ (case 5). Therefore, we obtain $\mathrm{id}^{s}_{\ell}(+0,+\infty,{\mathrm{n}})=+0$, because any number in $Y$ except $+0$ yields $+\infty$ when divided by $+0$. If we compute intermediate values exactly, $\mathrm{id}^{s}_{u}(42,6)=7$ and the refined interval is $Z^{\prime}=[+0,7]$. If not, then $z^{\prime}_{u}={1.110\cdots 01}\times{2^{2}}=\mathop{\mathrm{succ}}\nolimits(7)$.

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 BagnaraCGG16IJOC . 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.

Example 9

Consider the IEEE 754 single-precision division constraint $x=y\boxslash_{S}z$ with initial intervals $x\in[{1.0\cdots 010}\times{2^{110}},{1.0}\times{2^{121}}]$ and $Y=Z=[-\infty,+\infty]$. When $S=\{{\mathrm{n}}\}$, filtering by maximum ULP results in the possible refinement $Z^{\prime}=[-{1.0}\times{2^{18}},{1.0}\times{2^{18}}]$, while Algorithm 5 would compute $Z^{\prime}=[-{f_{\mathrm{max}}},{f_{\mathrm{max}}}]$, regardless of the rounding mode.

5 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.

5.1 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 WuLZ17 . 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 self-developed 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.

5.2 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 suite666https://github.com/florianschanda/smtlib_schanda, last accessed on October 28th, 2021. used in BrainSS19 .)

• •

11,544 are randomly generated tests by Christoph M. Wintersteiger. (QF_FP/wintersteiger directory from the SMT-LIB benchmark repository.777https://smtlib.cs.uiowa.edu/benchmarks.shtml, last accessed on October 28th, 2021.)

• •

126,909 tests were generated from the IBM Test Suite for IEEE 754R Compliance888https://www.research.ibm.com/haifa/projects/verification/fpgen/test_suite_download.shtml, last accessed on October 28th, 2021. created with FPgen AharoniAFKN03 .

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 MullerP00 . 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 Prolog999https://eclipseclp.org/, last accessed on October 28th, 2021. and SWI Prolog,101010https://www.swi-prolog.org/, last accessed on October 28th, 2021. while other tools were written in C**++**.

6 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 extended-precision 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 BagnaraCGG16IJOC , 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 Gallois-WongBC20 , 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.

Acknowledgements.

The authors express their gratitude to Roberto Amadini, Isacco Cattabiani and Laura Savino for their careful reading of early versions of this paper.

Appendix A Filtering algorithms: Subtraction and Multiplication

A.1 Subtraction

Here we deal with constraints of the form $x=y\boxminus_{S}z$.

Assume $X=[x_{\ell},x_{u}]$, $Y=[y_{\ell},y_{u}]$ and $Z=[z_{\ell},z_{u}]$.

Again, thanks to Proposition  2 we need not be concerned with sets of rounding modes, as any such set $S\subseteq R$ can always be mapped to a pair of “worst-case rounding modes” which, in addition are never round-to-zero.

Direct Propagation.

For direct propagation, we use Algorithm 6 and functions $\mathrm{ds}_{\ell}$ and $\mathrm{ds}_{u}$, as defined in Figure 12.

Theorem A.1

Algorithm 6* satisfies its contract.*

Inverse Propagation.

For inverse propagation, we have to deal with two different cases depending on which variable we are computing: the first inverse projection on $y$ or the second inverse projection on $z$.

The first inverse projection of subtraction is somehow similar to the direct projection of addition. In this case we define Algorithm 7 and functions $\mathrm{is}^{f}_{\ell}$ and $\mathrm{is}^{f}_{u}$, as defined in Figure 13 and 14 respectively.

Theorem A.2

Algorithm 7* satisfies its contract.*

The second inverse projection of subtraction is quite similar to the case of direct projection of subtraction. Here we define Algorithm 8 and functions $\mathrm{is}^{s}_{\ell}$ and $\mathrm{is}^{s}_{u}$, as defined in Figures 15 and 16 respectively.

Theorem A.3

Algorithm 8* is correct.*

Since subtraction is very closely related to addition, the proofs of Theorems A.2 and A.3 can be obtained by reasoning in the same way as for the projections of addition. Moreover, it is worth noting that in order to obtain more precise results, inverse projections for subtraction need to be intersected with maximum ULP filtering BagnaraCGG16IJOC , as in the case of addition.

A.2 Multiplication

Here we deal with constraints of the form $x=y\boxdot_{S}z$. As usual, assume $X=[x_{\ell},x_{u}]$, $Y=[y_{\ell},y_{u}]$ and $Z=[z_{\ell},z_{u}]$.

Direct Propagation.

For direct propagation, a case analysis is performed in order to select the interval extrema $y_{L}$ and $z_{L}$ (resp., $y_{U}$ and $z_{U}$) to be used to compute the new lower (resp., upper) bound for $x$.

Firstly, whenever $\mathop{\mathrm{sgn}}\nolimits(y_{\ell})\neq\mathop{\mathrm{sgn}}\nolimits(y_{u})$ and $\mathop{\mathrm{sgn}}\nolimits(z_{\ell})\neq\mathop{\mathrm{sgn}}\nolimits(z_{u})$, there is no unique choice for $y_{L}$ and $z_{L}$ (resp., $y_{U}$ and $z_{U}$); therefore we need to compute the two candidate lower (and upper) bounds for $x$ and then choose the minumum (the maximum, resp).

The choice is instead unique in all cases where the signs of one among $y$ and $z$, or both of them, are constant over the respective intervals. Function $\sigma$ of Figure 7 determines the extrema of $y$ and $z$ useful to compute the new lower (resp., upper) bound for $y$ when the sign of $z$ is constant. When the sign of $y$ is constant, the appropriate choice for the extrema of $y$ and $z$ can be determined by swapping the role of $y$ and $z$ in function $\sigma$.

Once the extrema $(y_{L},y_{U},z_{L},z_{U})$ have been selected, functions $\mathrm{dm}_{\ell}$ and $\mathrm{dm}_{u}$ of Figure 17 are used to find new bounds for $x$. It is worth noting that it is not necessary to compute new values of $r_{\ell}$ and $r_{u}$ for the application of functions $\mathrm{dm}_{\ell}$ and $\mathrm{dm}_{u}$ at line 6 of Algorithm 9. This is true because, by Definition 7, the choice of $r_{\ell}$ (of $r_{u}$, resp.) is driven by the sign of $y_{L}\boxdot z_{L}$ (of $y_{U}\boxdot z_{U}$, resp.) only. Since, in this case, the sign of $y_{L}\boxdot z_{L}$ (of $y_{U}\boxdot z_{U}$, resp.) as defined at line 2 and the sign of $y_{L}\boxdot z_{L}$ (of $y_{U}\boxdot z_{U}$, resp.) as defined at line 5 are the same, we do not need to compute $r_{\ell}$ and $r_{u}$ another time.

Theorem A.4

Algorithm 9* satisfies its contract.*

Inverse Propagation.

For inverse propagation, Algorithm 10 partitions interval $Z$ into the sign-homogeneous intervals $Z_{-}\mathrel{\buildrel\mathrm{def}\over{=}}Z\cap[-\infty,-0]$ and $Z_{+}\mathrel{\buildrel\mathrm{def}\over{=}}Z\cap[+0,+\infty]$. This is done because the sign of $Z$ must be taken into account in order to derive correct bounds for $Y$. Hence, once $Z$ has been partitioned into sign-homogeneous intervals, we use intervals $X$ and $Z_{-}$ to obtain interval $[y^{-}_{\ell},y^{-}_{u}]$, and $X$ and $Z_{+}$ to obtain $[y^{+}_{\ell},y^{+}_{u}]$. To do so, the algorithm determines the appropriate extrema of intervals $X$ and $W=Z_{-}$ or $W=Z_{+}$ to be used for constraint propagation. To this aim, function $\tau$ of Figure 5 is employed; note that the sign of $W$ is, by construction, constant over the interval. The chosen extrema are then passed as parameters to functions $\mathrm{im}_{\ell}$ of Figure 18 and $\mathrm{im}_{u}$ of Figure 19, that compute the new, refined bounds for $y$, by using the inverse operation of multiplication, i.e., division. The so obtained intervals $Y\cap[y^{-}_{\ell},y^{-}_{u}]$ and $Y\cap[y^{+}_{\ell},y^{+}_{u}]$ will be then joined with convex union, denoted by $\biguplus$, to obtain $Y^{\prime}$.

Theorem A.5

Algorithm 10* satisfies its contract.*

Of course, the refinement $Z^{\prime}$ of $Z$ can be defined analogously.

Appendix B Proofs of Results

B.1 Proofs of Results in Section 2

Proof (of Proposition 1)

In order to prove (5), we first prove that $[x]_{{\mathord{\downarrow}}}\leq x$. To this aim, consider the following cases on $x\in\mathbb{R}\setminus\{0\}$:

$-{f_{\mathrm{max}}}\leq x<0\;\lor\;{f_{\mathrm{min}}}\leq x:$

by (2) we have $[x]_{{\mathord{\downarrow}}}=\max\{\,z\in\mathbb{F}\mid z\leq x\,\}$, hence $[x]_{{\mathord{\downarrow}}}\leq x$;

$0<x<{f_{\mathrm{min}}}:$

by (2) we have $[x]_{{\mathord{\downarrow}}}=-0\leq x$;

$x<-{f_{\mathrm{max}}}:$

by (2) we have $[x]_{{\mathord{\downarrow}}}=-\infty\leq x$.

We now prove that $x\leq[x]_{{\mathord{\uparrow}}}$. Consider the following cases on $x\in\mathbb{R}\setminus\{0\}$:

$x>{f_{\mathrm{max}}}:$

by (1) we have $[x]_{{\mathord{\uparrow}}}=+\infty$ and thus $x\leq[x]_{{\mathord{\uparrow}}}$ holds;

$x\leq-{f_{\mathrm{min}}}\;\lor\;0<x\leq{f_{\mathrm{max}}}:$

by (1) we have $[x]_{{\mathord{\uparrow}}}=\min\{\,z\in\mathbb{F}\mid z\geq x\,\}$, hence $x\leq[x]_{{\mathord{\uparrow}}}$ holds;

$-{f_{\mathrm{min}}}<x<0:$

by (1) we have $[x]_{{\mathord{\uparrow}}}=-0$ hence $x\leq[x]_{{\mathord{\uparrow}}}$ holds.

In order to prove (6), consider the following cases on $x\in\mathbb{R}\setminus\{0\}$:

$x>0:$

by (3) we have $[x]_{{0}}=[x]_{{\mathord{\downarrow}}}\leq[x]_{{\mathord{\uparrow}}}$;

$x<0:$

by (3) we have $[x]_{{\mathord{\downarrow}}}\leq[x]_{{\mathord{\uparrow}}}=[x]_{{0}}$.

In order to prove (7), consider the following cases on $x\in\mathbb{R}\setminus\{0\}$:

$-{f_{\mathrm{max}}}\leq x\leq{f_{\mathrm{max}}}:$

we have the following cases

\bigl{|}[x]_{{\mathord{\downarrow}}}-x\bigr{|}<\bigl{|}[x]_{{\mathord{\uparrow}}}-x\bigr{|}\vee\bigl{(}\bigl{|}[x]_{{\mathord{\downarrow}}}-x\bigr{|}=\bigl{|}[x]_{{\mathord{\uparrow}}}-x\bigr{|}\bigr{)}\wedge\mathop{\mathrm{even}}\nolimits\bigl{(}[x]_{{\mathord{\downarrow}}}\bigr{)}:

by (4), we have $[x]_{{\mathrm{n}}}=[x]_{{\mathord{\downarrow}}}\leq[x]_{{\mathord{\uparrow}}}$;

\bigl{|}[x]_{{\mathord{\downarrow}}}-x\bigr{|}>\bigl{|}[x]_{{\mathord{\uparrow}}}-x\bigr{|}\vee\bigl{(}\bigl{|}[x]_{{\mathord{\downarrow}}}-x\bigr{|}=\bigl{|}[x]_{{\mathord{\uparrow}}}-x\bigr{|}\wedge\neg\mathop{\mathrm{even}}\nolimits\bigl{(}[x]_{{\mathord{\downarrow}}}\bigr{)}\bigr{)}:

by (4) we have $[x]_{{\mathord{\downarrow}}}\leq[x]_{{\mathord{\uparrow}}}=[x]_{{\mathrm{n}}}$.

$-{f_{\mathrm{max}}}>x:$

we have the following cases

$-2^{{e_{\mathrm{max}}}}(2-2^{-p})<x<-{f_{\mathrm{max}}}$:

by (4) we have $[x]_{{\mathord{\downarrow}}}\leq[x]_{{\mathord{\uparrow}}}=[x]_{{\mathrm{n}}}$.

$x\leq-2^{{e_{\mathrm{max}}}}(2-2^{-p})$:

by (4) we have $[x]_{{\mathrm{n}}}=[x]_{{\mathord{\downarrow}}}\leq[x]_{{\mathord{\uparrow}}}$;

${f_{\mathrm{max}}}<x:$

we have the following cases

$2^{{e_{\mathrm{max}}}}(2-2^{-p})>x>{f_{\mathrm{max}}}$:

by (4) we have $[x]_{{\mathrm{n}}}=[x]_{{\mathord{\downarrow}}}\leq[x]_{{\mathord{\uparrow}}}$;

$x\geq 2^{{e_{\mathrm{max}}}}(2-2^{-p})$:

by (4) we have $[x]_{{\mathord{\downarrow}}}\leq[x]_{{\mathord{\uparrow}}}=[x]_{{\mathrm{n}}}$.

In order to prove (8), let us compute $-[-x]_{{\mathord{\uparrow}}}$. There are the following cases:

$-x>{f_{\mathrm{max}}}:$

this implies that $x<-{f_{\mathrm{max}}}$ and, by (1), $[-x]_{{\mathord{\uparrow}}}=+\infty$; hence, by (2), $-[-x]_{{\mathord{\uparrow}}}=-\infty=[x]_{{\mathord{\downarrow}}}$;

$-x\leq-{f_{\mathrm{min}}}\;\lor\;0<-x\leq{f_{\mathrm{max}}}:$

this implies that $x\geq{f_{\mathrm{min}}}\;\lor\;-{f_{\mathrm{max}}}\geq x>0$ and, by (1), we have $[-x]_{{\mathord{\uparrow}}}=\min\{\,z\in\mathbb{F}\mid z\geq-x\,\}$; therefore, by (2), $-[-x]_{{\mathord{\uparrow}}}=-\min\{\,z\in\mathbb{F}\mid z\geq-x\,\}=\max\{\,z\in\mathbb{F}\mid z\leq x\,\}=[x]_{{\mathord{\downarrow}}}$,

$-{f_{\mathrm{min}}}<-x<0:$

this implies that $0<x<{f_{\mathrm{min}}}$ and, by (1), $[-x]_{{\mathord{\uparrow}}}=-0$; hence, by (2), $-[-x]_{{\mathord{\uparrow}}}=+0=[x]_{{\mathord{\downarrow}}}$.

B.2 Proofs of Results in Section 3

Proof (Rest of the proof of Proposition 2)

We prove the second part of Proposition 2, regarding rounding mode selectors for inverse propagators. Before doing so, we need to prove the following result. Let $\mathord{\boxcircle}\in\{\mathord{\boxplus},\mathord{\boxminus},\mathord{\boxdot},\mathord{\boxslash}\},$ and let $\mathrm{r}$ and $\mathrm{s}$ be two IEEE 754 rounding modes, such that for any $a,b\in\mathbb{F}$,

[TABLE]

Moreover, let $x,z\in\mathbb{F}$, and let $\bar{y}_{\mathrm{s}}$ be the minimum $y_{\mathrm{s}}\in\mathbb{F}$ such that $x=y_{\mathrm{s}}\boxcircle_{\mathrm{s}}z$. Then, for any $y_{\mathrm{r}}\in\mathbb{F}$ such that $x=y_{\mathrm{r}}\boxcircle_{\mathrm{r}}z$ we have

[TABLE]

This leads us to write

[TABLE]

which, due to the isotonicity of all IEEE 754 rounding modes, implies

[TABLE]

Finally, if operator ‘$\circ$’ is isotone we have

[TABLE]

which implies that $\bar{y}_{\mathrm{s}}$ is the minimum $y\in\mathbb{F}$ such that $x=y\boxcircle_{r}z$ or $x=y\boxcircle_{s}z$. On the other hand, if ‘$\circ$’ is antitone we have

[TABLE]

and $\bar{y}_{\mathrm{s}}$ is the maximum $y\in\mathbb{F}$ such that $x=y\boxcircle_{r}z$ or $x=y\boxcircle_{s}z$. An analogous result can be proved regarding the upper bound for $y$ in case the operator is isotone, and regarding the lower bound for $y$ in case it is antitone.

The above claim allows us to prove the following. Assume first that $\boxcircle$ is isotone with respect to $y$ in $x=y\boxcircle z$. Let $\hat{y}_{{\mathord{\uparrow}}}$ be the minimum $y_{{\mathord{\uparrow}}}\in\mathbb{F}$ such that $x=y_{{\mathord{\uparrow}}}\boxcircle_{\mathord{\uparrow}}z=[y_{\uparrow}\boxcircle z]_{{\mathord{\uparrow}}}$, let $\hat{y}_{\mathrm{n}}$ be the minimum $y_{\mathrm{n}}\in\mathbb{F}$ such that $x=y_{{\mathrm{n}}}\boxcircle_{\mathrm{n}}z=[y_{\mathrm{n}}\boxcircle z]_{{\mathrm{n}}}$ and, finally, let $\hat{y}_{\mathord{\downarrow}}$ be the minimum $y_{\mathord{\downarrow}}\in\mathbb{F}$ such that $x=y_{{\mathord{\downarrow}}}\boxcircle_{\mathord{\downarrow}}z=[y_{\mathord{\downarrow}}\boxcircle z]_{{\mathord{\downarrow}}}$. We will prove that

[TABLE]

Since we assumed that $\boxcircle$ is isotone with respect to $y$ in $x=y\boxcircle z$, the rounding mode that gives the minimal $y$ solution of $x=[y\boxcircle z]_{\mathrm{r}}$ is the one that yields a bigger (w.r.t. $\preccurlyeq$ order) floating point number, as we proved before. We must now separately treat the following cases:

$y\boxcircle z\not=0:$

By (7), we have $[y\boxcircle z]_{{\mathord{\downarrow}}}\leq[y\boxcircle z]_{{\mathrm{n}}}\leq[y\boxcircle z]_{{\mathord{\uparrow}}}$. Since in this case $y\boxcircle z\not=0$, we have that $[y\boxcircle z]_{{\mathord{\downarrow}}}\preccurlyeq[y\boxcircle z]_{{\mathrm{n}}}\preccurlyeq[y\boxcircle z]_{{\mathord{\uparrow}}}$. This implies $\hat{y}_{\mathord{\uparrow}}\preccurlyeq\hat{y}_{\mathrm{n}}\preccurlyeq\hat{y}_{\mathord{\downarrow}}$.

$y\boxcircle z=0:$

In this case, $[y\boxcircle z]_{{\mathord{\downarrow}}}\preccurlyeq[y\boxcircle z]_{{\mathrm{n}}}=[y\boxcircle z]_{{\mathord{\uparrow}}}$. This implies $\hat{y}_{\mathord{\uparrow}}\preccurlyeq\hat{y}_{\mathrm{n}}\preccurlyeq\hat{y}_{\mathord{\downarrow}}$.

Moreover, let $\tilde{y}_{{\mathord{\uparrow}}}$ be the maximum $y_{{\mathord{\uparrow}}}\in\mathbb{F}$ such that $x=y_{{\mathord{\uparrow}}}\boxcircle_{\mathord{\uparrow}}z=[y_{\uparrow}\boxcircle z]_{{\mathord{\uparrow}}}$, let $\tilde{y}_{\mathrm{n}}$ be the maximum $y_{\mathrm{n}}\in\mathbb{F}$ such that $x=y_{{\mathrm{n}}}\boxcircle_{\mathrm{n}}z=[y_{\mathrm{n}}\boxcircle z]_{{\mathrm{n}}}$ and, finally, let $\tilde{y}_{\mathord{\downarrow}}$ be the maximum $y_{\mathord{\downarrow}}\in\mathbb{F}$ such that $x=y_{{\mathord{\downarrow}}}\boxcircle_{\mathord{\downarrow}}z=[y_{\mathord{\downarrow}}\boxcircle z]_{{\mathord{\downarrow}}}$. We will prove the fact that

[TABLE]

Since we assumed that $\boxcircle$ is isotone with respect to $y$ in $x=y\boxcircle z$, the rounding mode that gives a maximum $y$ solution of $x=y\boxcircle z_{\mathrm{r}}$ is the one that gives a smaller (w.r.t. $\preccurlyeq$ order) floating point number. We must now deal with the following cases:

$y\boxcircle z\not=0:$

By (7), we have $[y\boxcircle z]_{{\mathord{\downarrow}}}\leq[y\boxcircle z]_{{\mathrm{n}}}\leq[y\boxcircle z]_{{\mathord{\uparrow}}}$. Since in this case $y\boxcircle z\not=0$, we have $[y\boxcircle z]_{{\mathord{\downarrow}}}\preccurlyeq[y\boxcircle z]_{{\mathrm{n}}}\preccurlyeq[y\boxcircle z]_{{\mathord{\uparrow}}}$. This implies $\tilde{y}_{\mathord{\uparrow}}\preccurlyeq\tilde{y}_{\mathrm{n}}\preccurlyeq\tilde{y}_{\mathord{\downarrow}}$.

$y\boxcircle z=0:$

In this case $[y\boxcircle z]_{{\mathord{\downarrow}}}\preccurlyeq[y\boxcircle z]_{{\mathrm{n}}}=[y\boxcircle z]_{{\mathord{\uparrow}}}$. This implies $\tilde{y}_{\mathord{\uparrow}}\preccurlyeq\tilde{y}_{\mathrm{n}}\preccurlyeq\tilde{y}_{\mathord{\downarrow}}$.

The inequalities $\hat{y}_{\mathord{\uparrow}}\preccurlyeq\hat{y}_{\mathrm{n}}\preccurlyeq\hat{y}_{\mathord{\downarrow}}$ and $\tilde{y}_{\mathord{\uparrow}}\preccurlyeq\tilde{y}_{\mathrm{n}}\preccurlyeq\tilde{y}_{\mathord{\downarrow}}$ allow us to claim that the rounding mode selectors $\hat{r}_{\ell}(S,\boxcircle,b)$ and $\hat{r}_{u}(S,b)$ are correct when $\boxcircle$ is isotone with respect to $y$. In a similar way it is possible to prove that, in case $\boxcircle$ is antitone with respect to argument $y$, the above-mentioned rounding mode selectors can be exchanged: $\hat{r}_{u}(S,b)$ can be used to obtain the lower bound for $y$, while $\hat{r}_{\ell}(S,\boxcircle,b)$ can be used to obtain the upper bound.

Note that, in general, the roundTowardZero rounding mode is equivalent to roundTowardPositive if the result of the rounded operation is negative, and to roundTowardNegative if it is positive. The only case in which this is not true is when the result is $+0$ and the operation is a sum or a subtraction: this value can come from the rounding toward negative infinity of a strictly positive exact result, or the sum of $+0$ and $-0$, which behaves like roundTowardPositive, yielding $+0$. This case must be treated separately, and it is significant only in $\hat{r}_{\ell}(S,\boxcircle,b)$, which is used when seeking for the lowest possible value of the variable to be refined that yields $+0$.

Definition 8 also contains selectors that can choose between rounding mode selectors $\hat{r}_{\ell}(S,b)$ and $\hat{r}_{u}(S,b)$ by distinguishing whether the operator is isotone or antitone with respect to the operand $y$ to be derived by propagation; they take the result of the operation $b$ and the known operand $a$ into account. In particular, $\bar{r}_{\ell}^{\ell}(S,b,\boxcircle,a)$, $\bar{r}_{u}^{\ell}(S,b,\boxcircle,a)$ choose the appropriate selector for the leftmost operand, and $\bar{r}_{\ell}^{r}(S,b,\boxcircle,a)$, $\bar{r}_{u}^{r}(S,b,\boxcircle,a)$ are valid for the rightmost one.

Proof (of Proposition 3)

We first prove (15). By Definition 9, we have the following cases:

$x_{\ell}=-{f_{\mathrm{max}}}$:

Then,

[TABLE]

On the other hand, consider any $x$ such that $x_{\ell}<x\leq x_{u}$. Since $x\in\mathbb{F}$, this implies that $\mathop{\mathrm{succ}}\nolimits(-{f_{\mathrm{max}}})\leq x\leq x_{u}$. In this case

[TABLE]

Since ‘$\mathop{\mathrm{pred}}\nolimits$’ is monotone, the minimum can be found when $x=\mathop{\mathrm{succ}}\nolimits(-{f_{\mathrm{max}}})$. In this case, we have that

[TABLE]

Hence we can conclude that \min_{x_{\ell}\leq x\leq x_{u}}\bigl{(}x+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)/2\bigr{)}=x_{\ell}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x_{\ell})/2.

$x_{\ell}>-{f_{\mathrm{max}}}$:

In this case

[TABLE]

Since ‘$\mathop{\mathrm{pred}}\nolimits$’ is monotone, \min_{x_{\ell}\leq x\leq x_{u}}\bigl{(}x+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x_{\ell})/2\bigr{)}=x_{\ell}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x_{\ell})/2.

We now prove (16). By Definition 9, we have the following cases:

$x_{u}={f_{\mathrm{max}}}$:

Then,

[TABLE]

Now, consider any $x$ such that $x_{\ell}\leq x<x_{u}$. Since $x\in\mathbb{F}$, this implies that $x_{\ell}\leq x\leq\mathop{\mathrm{pred}}\nolimits({f_{\mathrm{max}}})$. In this case

[TABLE]

Since ‘$\mathop{\mathrm{succ}}\nolimits$’ is monotone, the maximum can be found when $x=\mathop{\mathrm{pred}}\nolimits({f_{\mathrm{max}}})$. In this case, we have that

[TABLE]

Hence we can conclude that \max_{x_{\ell}\leq x\leq x_{u}}\bigl{(}x+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)/2\bigr{)}=x_{u}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x_{u})/2.

$x_{u}<{f_{\mathrm{max}}}$:

In this case

[TABLE]

Since ‘$\mathop{\mathrm{succ}}\nolimits$’ is monotone, \max_{x_{\ell}\leq x\leq x_{u}}\bigl{(}x+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)/2\bigr{)}=x_{u}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x_{u})/2.

We now introduce and prove Proposition 6, which contains properties of the rounding error functions that are only needed in the proof of Proposition 4.

Proposition 6

For each $r\in\mathbb{R}\setminus\{0\}$ we have

[TABLE]

where the two inequalities of (45) are strict if $[r]_{{\mathrm{n}}}$ is odd.

Proof

Suppose $r\in\mathbb{R}$ was rounded down to $x\in\mathbb{F}$. Then the error that was committed, $r-x$, is a nonnegative extended real that is strictly bounded from above by $\mathop{\nabla^{\downarrow}}\nolimits(x)=\mathop{\mathrm{succ}}\nolimits(x)-x$, that is, $0\leq r-x<\mathop{\mathrm{succ}}\nolimits(x)-x$, for otherwise we would have $r\geq\mathop{\mathrm{succ}}\nolimits(x)$ or $r<x$ and, in both cases $r$ would not have been rounded down to $x$. Note that $\mathop{\nabla^{\downarrow}}\nolimits({f_{\mathrm{max}}})=+\infty$, coherently with the fact that the error is unbounded from above in this case.

Dually, if $r\in\mathbb{R}$ was rounded up to $x\in\mathbb{F}$ the error that was committed, $r-x$, is a nonpositive extended real that is strictly bounded from below by $\mathop{\nabla^{\uparrow}}\nolimits(x)=\mathop{\mathrm{pred}}\nolimits(x)-x$, that is, $\mathop{\mathrm{pred}}\nolimits(x)-x<r-x\leq 0$ since, clearly, $\mathop{\mathrm{pred}}\nolimits(x)<r\leq x$. Note that $\mathop{\nabla^{\uparrow}}\nolimits(-{f_{\mathrm{max}}})=-\infty$, coherently with the fact that the error is unbounded from below in this case.

Suppose now that $r\in\mathbb{R}$ was rounded-to-nearest to $x\in\mathbb{F}$. Then the error that was committed, $r-x$, is such that $\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)/2\leq r-x\leq\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)/2$, where the two inequalities are strict if $x$ is odd.

In fact, if $x\notin\{-\infty,-{f_{\mathrm{max}}}\}$, then \mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)/2=\bigl{(}\mathop{\mathrm{pred}}\nolimits(x)-x\bigr{)}/2\leq r-x, for otherwise $r$ would be closer to $\mathop{\mathrm{pred}}\nolimits(x)$. If $x=-\infty$, then $\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)/2=+\infty$ and $r-x=+\infty$, so $\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)/2\leq r-x$ holds. If $x=-{f_{\mathrm{max}}}$, then

[TABLE]

and thus, considering that $-{f_{\mathrm{max}}}$ is odd, $\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)/2<r-x$ is equivalent to

[TABLE]

which must hold, for otherwise $r$ would have been rounded to $-\infty$ (IEEE-754-2008, , Section 4.3.1).

Suppose now $x\notin\{+\infty,{f_{\mathrm{max}}}\}$: then \mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)/2=\bigl{(}\mathop{\mathrm{succ}}\nolimits(x)-x\bigr{)}/2\geq r-x, for otherwise $r$ would be closer to $\mathop{\mathrm{succ}}\nolimits(x)$. If $x=+\infty$, then $\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)/2=-\infty$ and $r-x=-\infty$, and thus $\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)/2\geq r-x$ holds. If $x={f_{\mathrm{max}}}$, then

[TABLE]

and thus, considering that ${f_{\mathrm{max}}}$ is odd, $\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)/2>r-x$ is equivalent to

[TABLE]

which must hold, for otherwise $r$ would have been rounded to $+\infty$.

Proof (of Proposition 4)

In order to prove (17), first observe that $x\preccurlyeq y\boxcircle_{\mathord{\downarrow}}z$ implies that $x\leq y\boxcircle_{\mathord{\downarrow}}z$. Assume first that $y\boxcircle_{\mathord{\downarrow}}z\in\mathbb{R}_{+}\cup\mathbb{R}_{-}$. In this case, $y\boxcircle_{\mathord{\downarrow}}z=[y\circ z]_{{\mathord{\downarrow}}}$. By inequality (5) of Proposition 1, $y\boxcircle_{\mathord{\downarrow}}z=[y\circ z]_{{\mathord{\downarrow}}}\leq y\circ z$. Therefore, $x\leq y\boxcircle_{\mathord{\downarrow}}z=[y\circ z]_{{\mathord{\downarrow}}}\leq y\circ z$. Then, assume that $y\boxcircle_{\mathord{\downarrow}}z=+\infty$. In this case, since the rounding towards minus infinity never rounds to $+\infty$, it follows that $y\boxcircle_{\mathord{\downarrow}}z=y\circ z$. Hence, $x\leq y\circ z=+\infty$, holds. Assume now that $y\boxcircle_{\mathord{\downarrow}}z=-\infty$. In this case it must be that $x=-\infty$ then $x\leq y\circ z$, holds. Finally, assume that $y\boxcircle_{\mathord{\downarrow}}z=+0$ or $y\boxcircle_{\mathord{\downarrow}}z=-0$. In any case $x\preccurlyeq+0$ that implies $x\leq 0$. On the other hand, we have two cases, $y\circ z\neq 0$ or $y\circ z=0$. For the first case, by Definition 5, $0\leq y\circ z<{f_{\mathrm{min}}}$, then $x\leq y\circ z$, holds. For the second case, since $x\leq 0$ then $x\leq y\circ z$.

In order to prove (18), as before observe that $x\preccurlyeq y\boxcircle_{\mathord{\uparrow}}z$ implies that $x\leq y\boxcircle_{\mathord{\uparrow}}z$. Note that $x+\mathop{\nabla^{\uparrow}}\nolimits(x)=\mathop{\mathrm{pred}}\nolimits(x)$. So we are left to prove $\mathop{\mathrm{pred}}\nolimits(x)<y\circ z$. Assume now that $0<y\circ z\leq{f_{\mathrm{max}}}$ or $x\leq-{f_{\mathrm{min}}}$. Moreover, note that it cannot be the case that $\mathop{\mathrm{pred}}\nolimits(x)\geq y\circ z$, otherwise, by Definition 5, $y\boxcircle_{\mathord{\uparrow}}z\leq\mathop{\mathrm{pred}}\nolimits(x)$ and, therefore, $x\leq y\boxcircle_{\mathord{\uparrow}}z$ would not hold. Then, in this case, we can conclude $\mathop{\mathrm{pred}}\nolimits(x)<y\circ z$. Now, assume that $-{f_{\mathrm{min}}}<y\circ z<0$. In this case $y\boxcircle_{\mathord{\uparrow}}z=-0$. Hence, $x\leq 0$. By Definition 4, $\mathop{\mathrm{pred}}\nolimits(x)\leq-{f_{\mathrm{min}}}$. Hence, $\mathop{\mathrm{pred}}\nolimits(x)\leq y\circ z$, holds. Next, assume $y\circ z>{f_{\mathrm{max}}}$. In this case $y\boxcircle_{\mathord{\uparrow}}z=\infty$. Hence, $x\leq\infty$. By Definition 4, $\mathop{\mathrm{pred}}\nolimits(x)\leq{f_{\mathrm{max}}}$. Hence $\mathop{\mathrm{pred}}\nolimits(x)<y\circ z$, holds. Next assume $y\circ z=0$. In this case $y\boxcircle_{\mathord{\uparrow}}z=+0$ or $y\boxcircle_{\mathord{\uparrow}}z=-0$. Hence, $x\leq 0$. By Definition 4, $\mathop{\mathrm{pred}}\nolimits(x)\leq-{f_{\mathrm{min}}}$. Hence $\mathop{\mathrm{pred}}\nolimits(x)<y\circ z$, holds. Finally assume $y\circ z=\infty$. In this case $y\boxcircle_{\mathord{\uparrow}}z=\infty$. Hence $x\preccurlyeq\infty$ and therefore $x\leq\infty$. By Definition 4, $\mathop{\mathrm{pred}}\nolimits(x)\leq{f_{\mathrm{max}}}$. Hence $\mathop{\mathrm{pred}}\nolimits(x)<y\circ z$, holds.

In order to prove (19), as the previous two cases, note that $x\preccurlyeq y\boxcircle_{\mathrm{n}}z$ implies that $x\leq y\boxcircle_{\mathrm{n}}z$. First observe that for $x\neq-\infty$, $x+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)/2<x$. Indeed, assume first that $x\neq-{f_{\mathrm{max}}}$, then, by Definition 9, $\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)=x-\mathop{\mathrm{succ}}\nolimits(x)$. Hence $x+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)/2=x+(x-\mathop{\mathrm{succ}}\nolimits(x))/2=(3x-\mathop{\mathrm{succ}}\nolimits(x))/2$. Since $x<\mathop{\mathrm{succ}}\nolimits(x)$, we can conclude that $x+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)/2<x$. Assume now that $x=-{f_{\mathrm{max}}}$. By Definition 9, $\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)=\mathop{\mathrm{pred}}\nolimits(x)-x$. Hence $x+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)/2=x+(\mathop{\mathrm{pred}}\nolimits(x)-x)/2=(x+\mathop{\mathrm{pred}}\nolimits(x))/2$. Since $x>\mathop{\mathrm{pred}}\nolimits(x)$, we can conclude that $x+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)/2<x$.

Now, by Definition 5, we have to consider the following cases for $x\boxcircle_{\mathrm{n}}y\in\mathbb{R}_{+}\cup\mathbb{R}_{-}$:

$y\boxcircle_{\mathrm{n}}z=[y\circ z]_{{\mathord{\downarrow}}}$. In this case, by inequality (5) of Proposition 1, $x+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)/2<x\leq y\boxcircle_{\mathrm{n}}z=[y\circ z]_{{\mathord{\downarrow}}}\leq y\circ z$. Therefore, $x+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)/2<y\circ z$, holds.

$y\boxcircle_{\mathrm{n}}z=[y\circ z]_{{\mathord{\uparrow}}}$. Assume first that $x<y\boxcircle_{\mathrm{n}}z$. In this case, by Definition 5, since $x\in\mathbb{F}$ and $x<y\boxcircle_{\mathrm{n}}z$, it must be the case that $x<y\circ z$. Then, we can conclude that $x+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)/2<x<y\circ z$. Therefore, $x+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)/2<y\circ z$, holds. Assume now that $x=y\boxcircle_{\mathrm{n}}z$ and $\mathop{\mathrm{even}}\nolimits(x)$. In this case, by Proposition 6, we have that \mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits\bigl{(}[y\circ z]_{{\mathrm{n}}}\bigr{)}/2\leq(y\circ z)-[y\circ z]_{{\mathrm{n}}}. Since, in this case $x=y\boxcircle_{\mathrm{n}}z$, we obtain $\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)/2\leq(y\circ z)-x$. Hence, $x+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)/2\leq y\circ z$. If $\mathop{\mathrm{odd}}\nolimits(x)$, by Proposition 6, we have that \mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits\bigl{(}[y\circ z]_{{\mathrm{n}}}\bigr{)}/2<(y\circ z)-[y\circ z]_{{\mathrm{n}}}. Hence, $x+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)/2<y\circ z$.

Consider now the case that $y\boxcircle_{\mathrm{n}}z=+0$ or $y\boxcircle_{\mathrm{n}}z=-0$. If $y\circ z\neq 0$, then $y\boxcircle_{\mathrm{n}}z=[y\circ z]_{{\mathord{\downarrow}}}$ or $y\boxcircle_{\mathrm{n}}z=[y\circ z]_{{\mathord{\uparrow}}}$. In this case we can reason as above. Assume then that $y\circ z=0$. Since $x\preccurlyeq+0$ or $x\preccurlyeq-0$ implies that $x\leq 0$. Therefore, we can conclude that $x+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)<x\leq 0$ holds. Assume now that $y\boxcircle_{\mathrm{n}}z=+\infty$. If $y\circ z\neq\infty$ then $y\boxcircle_{\mathrm{n}}z=[y\circ z]_{{\mathord{\uparrow}}}$. In this case we can reason as above. On the other hand if $y\circ z=+\infty$ then $x+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x)\leq\infty$ holds.

In order to prove (20), remember that $x\succcurlyeq y\boxcircle_{\mathord{\downarrow}}z$ implies that $x\geq y\boxcircle_{\mathord{\downarrow}}z$. Note that $x+\mathop{\nabla^{\downarrow}}\nolimits(x)=\mathop{\mathrm{succ}}\nolimits(x)$. So we are left to prove $\mathop{\mathrm{succ}}\nolimits(x)>y\circ z$. Assume now that $-{f_{\mathrm{max}}}<y\circ z<0$ or ${f_{\mathrm{min}}}<y\circ z\leq{f_{\mathrm{max}}}$. Note that it cannot be the case that $\mathop{\mathrm{succ}}\nolimits(x)\leq y\circ z$, otherwise, by Definition 5, $y\boxcircle_{\mathord{\downarrow}}z\geq\mathop{\mathrm{succ}}\nolimits(x)$ and $x\geq y\boxcircle_{\mathord{\downarrow}}z$ would not hold. Then, in this case, we can conclude that $\mathop{\mathrm{succ}}\nolimits(x)>y\circ z$. Next, assume that $0<y\circ z<{f_{\mathrm{min}}}$. In this case $y\boxcircle_{\mathord{\downarrow}}z=+0$. Hence, $x\geq 0$. By Definition 4, $\mathop{\mathrm{succ}}\nolimits(x)\geq{f_{\mathrm{min}}}$. Hence $\mathop{\mathrm{succ}}\nolimits(x)\geq y\circ z$, holds. Next, assume $y\circ z<-{f_{\mathrm{max}}}$. In this case $y\boxcircle_{\mathord{\downarrow}}z=-\infty$. Hence $x\geq-\infty$. By Definition 4, $\mathop{\mathrm{succ}}\nolimits(x)\geq-{f_{\mathrm{max}}}$. Hence $\mathop{\mathrm{succ}}\nolimits(x)>y\circ z$, holds. Next assume $y\circ z=0$. In this case $y\boxcircle_{\mathord{\downarrow}}z=+0$ or $y\boxcircle_{\mathord{\downarrow}}z=-0$. In any case, $x\geq 0$. By Definition 4, $\mathop{\mathrm{succ}}\nolimits(x)\geq{f_{\mathrm{min}}}$. Hence $\mathop{\mathrm{succ}}\nolimits(x)>y\circ z$, holds. Finally assume $y\circ z=-\infty$. In this case $y\boxcircle_{\mathord{\downarrow}}z=-\infty$. Hence, since $x\succcurlyeq-\infty$, $x\geq-\infty$. By Definition 4, $\mathop{\mathrm{succ}}\nolimits(x)\geq-{f_{\mathrm{max}}}$. Hence, $\mathop{\mathrm{succ}}\nolimits(x)>y\circ z$, holds.

In order to prove (21), as before, observe that $x\succcurlyeq y\boxcircle_{\mathord{\uparrow}}z$ implies that $x\geq y\boxcircle_{\mathord{\uparrow}}z$. Assume first that $y\boxcircle_{\mathord{\uparrow}}z\in\mathbb{R}_{+}\cup\mathbb{R}_{-}$. In this case, $y\boxcircle_{\mathord{\uparrow}}z=[y\circ z]_{{\mathord{\uparrow}}}$. By (5) from Proposition 1, $y\boxcircle_{\mathord{\uparrow}}z=[y\circ z]_{{\mathord{\uparrow}}}\geq y\circ z$. Then, assume that $y\boxcircle_{\mathord{\uparrow}}z=-\infty$. In this case, since the rounding towards plus infinity never rounds to $-\infty$, it follows that $y\boxcircle_{\mathord{\uparrow}}z=y\circ z$. Hence, $x\geq y\circ z=-\infty$, holds. Assume now that $y\boxcircle_{\mathord{\uparrow}}z=+\infty$. In this case, $x=+\infty$ then $x\geq y\circ z$, holds. Finally, assume that $y\boxcircle_{\mathord{\uparrow}}z=+0$ or $y\boxcircle_{\mathord{\uparrow}}z=-0$. In any case $x\succcurlyeq-0$ that implies $x\geq 0$. On the other hand, we have two cases, $y\circ z\neq 0$ or $y\circ z=0$. For the first case, by Definition 5, $-{f_{\mathrm{min}}}<y\circ z<0$, then $x\geq y\circ z$, holds. For the second case, since $x\geq 0$ then $x\geq y\circ z$.

In order to prove (22), note that $x\succcurlyeq y\boxcircle_{\mathrm{n}}z$ implies that $x\geq y\boxcircle_{\mathrm{n}}z$. First observe that for $x\neq+\infty$, $x+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)/2>x$. Indeed, assume first that $x\neq{f_{\mathrm{max}}}$, then, by Definition 9, $\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)=x-\mathop{\mathrm{pred}}\nolimits(x)$. Hence $x+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)/2=x+(x-\mathop{\mathrm{pred}}\nolimits(x))/2=(3x-\mathop{\mathrm{pred}}\nolimits(x))/2$. Since $x>\mathop{\mathrm{pred}}\nolimits(x)$, we can conclude that $x+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)/2>x$. Assume now that $x={f_{\mathrm{max}}}$. By Definition 9, $\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)=\mathop{\mathrm{succ}}\nolimits(x)-x$. Hence $x+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)/2=x+(\mathop{\mathrm{succ}}\nolimits(x)-x)/2=(x+\mathop{\mathrm{succ}}\nolimits(x))/2$. Since $x<\mathop{\mathrm{succ}}\nolimits(x)$, we can conclude that $x+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)/2>x$.

By Definition 5, we have to consider the following cases for $x\boxcircle_{\mathrm{n}}y\in\mathbb{R}_{+}\cup\mathbb{R}_{-}$:

$y\boxcircle_{\mathrm{n}}z=[y\circ z]_{{\mathord{\uparrow}}}$. In this case, by inequality (5) of Proposition 1, $x+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)/2>x\geq y\boxcircle_{\mathrm{n}}z=[y\circ z]_{{\mathord{\uparrow}}}\geq y\circ z$. Therefore, $x+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)/2>y\circ z$, holds.

$y\boxcircle_{\mathrm{n}}z=[y\circ z]_{{\mathord{\downarrow}}}$. Assume first that $x>y\boxcircle_{\mathrm{n}}z$. In this case, by Definition 5, since $x\in\mathbb{F}$ and $x>y\boxcircle_{\mathrm{n}}z$, it must be the case that $x>y\circ z$. Hence, as in the previous case, by inequality (5) Proposition 1, $x+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)/2>x>y\circ z$. Therefore, $x+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)/2>y\circ z$, holds. Assume now that $x=y\boxcircle_{\mathrm{n}}z$ and $\mathop{\mathrm{even}}\nolimits(x)$. In this case, by Proposition 6, we have that \mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits\bigl{(}[y\circ z]_{{\mathrm{n}}}\bigr{)}/2\geq(y\circ z)-[y\circ z]_{{\mathrm{n}}}. Since, in this case $x=y\boxcircle_{\mathrm{n}}z$, we obtain $\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)/2\geq(y\circ z)-x$. Hence, $x+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)/2\geq y\circ z$. If $\mathop{\mathrm{odd}}\nolimits(x)$, by Proposition 6, we have that \mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits\bigl{(}[y\circ z]_{{\mathrm{n}}}\bigr{)}/2>y\circ z-[y\circ z]_{{\mathrm{n}}}. Hence, $x+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)/2>y\circ z$.

Consider now the case that $y\boxcircle_{\mathrm{n}}z=+0$ or $y\boxcircle_{\mathrm{n}}z=-0$. If $y\circ z\neq 0$, then $y\boxcircle_{\mathrm{n}}z=[y\circ z]_{{\mathord{\downarrow}}}$ or $y\boxcircle_{\mathrm{n}}z=[y\circ z]_{{\mathord{\uparrow}}}$. In this case we can reason as above. Assume now that $y\circ z=0$. Since $x\succcurlyeq+0$ or $x\succcurlyeq-0$ implies that $x\geq 0$, we can conclude that $x+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)/2>x\geq 0$ holds. Assume now that $y\boxcircle_{\mathrm{n}}z=-\infty$. If $y\circ z\neq-\infty$ then $y\boxcircle_{\mathrm{n}}z=[y\circ z]_{{\mathord{\downarrow}}}$. In this case we can reason as above. On the other hand if $y\circ z=-\infty$ then $x+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x)/2\geq-\infty$ holds.

Proof (of Proposition 5)

We first prove (25). By inequality (5) from Proposition 1, $e\geq[e]_{{\mathord{\downarrow}}}$. Hence, $x\geq[e]_{{\mathord{\downarrow}}}$. Since by hypothesis, $e\in E_{\mathbb{F}}$ is an expression that evaluates on $\overline{\mathbb{R}}$ to a nonzero value, we have three cases:

$[e]_{{\mathord{\downarrow}}}\neq 0$ and $x\neq 0$:

In this case $x\geq[e]_{{\mathord{\downarrow}}}$ implies $x\succcurlyeq[e]_{{\mathord{\downarrow}}}$.

$[e]_{{\mathord{\downarrow}}}=+0$:

In this case, $0<e<{f_{\mathrm{min}}}$. Then, it must be the case that $x>0$. Therefore $x\succcurlyeq[e]_{{\mathord{\downarrow}}}$ holds.

$x=0$:

In this case $x$ must be strictly greater than $e$ since $e\in E_{\mathbb{F}}$ evaluates to a nonzero value. Therefore, $e<0$. Hence, by Definition 5, $[e]_{{\mathord{\downarrow}}}\leq-{f_{\mathrm{min}}}$. Then $x\succcurlyeq[e]_{{\mathord{\downarrow}}}$ holds.

In all cases, we have that $x\succcurlyeq[e]_{{\mathord{\downarrow}}}$. By Definition 10, we conclude that $x\succcurlyeq\llbracket e\rrbracket_{{\mathord{\downarrow}}}$.

We now prove (26). By inequality (5) from Proposition 1, as in the previous case, $e\geq[e]_{{\mathord{\downarrow}}}$. Hence, $x>[e]_{{\mathord{\downarrow}}}$. Since by hypothesis $e\in E_{\mathbb{F}}$ is an expression that evaluates on $\overline{\mathbb{R}}$ to a nonzero value, we have three cases:

$[e]_{{\mathord{\downarrow}}}\neq 0$ and $x\neq 0$:

In this case $x>[e]_{{\mathord{\downarrow}}}$ implies $x\succ[e]_{{\mathord{\downarrow}}}$.

$[e]_{{\mathord{\downarrow}}}=+0$:

In this case, $0<e<{f_{\mathrm{min}}}$. Hence, $x>0$. Therefore $x\succ[e]_{{\mathord{\downarrow}}}$ holds.

$x=0$:

In this case $x$ must be strictly greater than $e$ since $e\in E_{\mathbb{F}}$ evaluates to a nonzero value. Therefore, $e<0$. Hence, by Definition 5, $[e]_{{\mathord{\downarrow}}}\leq-{f_{\mathrm{min}}}$. Then $x\succ[e]_{{\mathord{\downarrow}}}$ holds.

In all cases, we have that $x\succ[e]_{{\mathord{\downarrow}}}$. By Definition 10, we conclude that $x\succ\llbracket e\rrbracket_{{\mathord{\downarrow}}}$. Then, by Definition 4, we have the following cases on $\llbracket e\rrbracket_{{\mathord{\downarrow}}}$:

$\llbracket e\rrbracket_{{\mathord{\downarrow}}}={f_{\mathrm{max}}}$:

In this case $\mathop{\mathrm{succ}}\nolimits(\llbracket e\rrbracket_{{\mathord{\downarrow}}})=+\infty$. Since $x\succ\llbracket e\rrbracket_{{\mathord{\downarrow}}}$, this implies that $x=+\infty$. Then $x\succcurlyeq\mathop{\mathrm{succ}}\nolimits(\llbracket e\rrbracket_{{\mathord{\downarrow}}})$, holds.

$-{f_{\mathrm{max}}}\leq\llbracket e\rrbracket_{{\mathord{\downarrow}}}<-{f_{\mathrm{min}}}$ or ${f_{\mathrm{min}}}\leq\llbracket e\rrbracket_{{\mathord{\downarrow}}}<{f_{\mathrm{max}}}$:

In this case $\mathop{\mathrm{succ}}\nolimits(\llbracket e\rrbracket_{{\mathord{\downarrow}}})=\min\{\,y\in\mathbb{F}\mid y>\llbracket e\rrbracket_{{\mathord{\downarrow}}}\,\}$. Since $x>\llbracket e\rrbracket_{{\mathord{\downarrow}}}$, $x\in\{\,y\in\mathbb{F}\mid y>\llbracket e\rrbracket_{{\mathord{\downarrow}}}\,\}$. Hence, $x\succcurlyeq\mathop{\mathrm{succ}}\nolimits(\llbracket e\rrbracket_{{\mathord{\downarrow}}})$, holds.

$\llbracket e\rrbracket_{{\mathord{\downarrow}}}=+0$ or $\llbracket e\rrbracket_{{\mathord{\downarrow}}}=-0$:

In this case $\mathop{\mathrm{succ}}\nolimits(\llbracket e\rrbracket_{{\mathord{\downarrow}}})={f_{\mathrm{min}}}$. Since $x>\llbracket e\rrbracket_{{\mathord{\downarrow}}}$, this implies that $x\geq{f_{\mathrm{min}}}$. Hence, $x\succcurlyeq\mathop{\mathrm{succ}}\nolimits(\llbracket e\rrbracket_{{\mathord{\downarrow}}})$, holds.

$\llbracket e\rrbracket_{{\mathord{\downarrow}}}=-{f_{\mathrm{min}}}$:

In this case $\mathop{\mathrm{succ}}\nolimits(\llbracket e\rrbracket_{{\mathord{\downarrow}}})=-0$. Since $x>\llbracket e\rrbracket_{{\mathord{\downarrow}}}=-{f_{\mathrm{min}}}$, $x\succcurlyeq-0$. Hence, $x\succcurlyeq\mathop{\mathrm{succ}}\nolimits(\llbracket e\rrbracket_{{\mathord{\downarrow}}})$, holds.

$\llbracket e\rrbracket_{{\mathord{\downarrow}}}=-\infty$:

In this case $\mathop{\mathrm{succ}}\nolimits(\llbracket e\rrbracket_{{\mathord{\downarrow}}})=-{f_{\mathrm{max}}}$. Since $x>\llbracket e\rrbracket_{{\mathord{\downarrow}}}=-\infty$, $x\succcurlyeq-{f_{\mathrm{max}}}$. Hence, $x\succcurlyeq\mathop{\mathrm{succ}}\nolimits(\llbracket e\rrbracket_{{\mathord{\downarrow}}})$, holds.

We now prove (27). By inequality (5) from Proposition 1, $e\leq[e]_{{\mathord{\uparrow}}}$. Hence, like before, $x\leq[e]_{{\mathord{\uparrow}}}$. Since by hypothesis $e\in E_{\mathbb{F}}$ is an expression that evaluates on $\overline{\mathbb{R}}$ to a nonzero value, we have three cases:

$[e]_{{\mathord{\uparrow}}}\neq 0$ and $x\neq 0$:

In this case $x\leq[e]_{{\mathord{\uparrow}}}$ implies $x\preccurlyeq[e]_{{\mathord{\uparrow}}}$.

$[e]_{{\mathord{\uparrow}}}=-0$:

In this case, $-{f_{\mathrm{min}}}<e<0$. Hence, $x<0$. Therefore $x\preccurlyeq[e]_{{\mathord{\uparrow}}}$ holds.

$x=0$:

In this case it must be the case that $x$ is strictly smaller than $e$, since $e\in E_{\mathbb{F}}$ evaluates to a nonzero value. Therefore, $e>0$. Hence, by Definition 5, $[e]_{{\mathord{\uparrow}}}\geq{f_{\mathrm{min}}}$. Then $x\preccurlyeq[e]_{{\mathord{\uparrow}}}$ holds.

In any case, $x\preccurlyeq[e]_{{\mathord{\uparrow}}}$ holds. By Definition 10, we conclude that $x\preccurlyeq\llbracket e\rrbracket_{{\mathord{\uparrow}}}$.

Next we prove (28). By, again, inequality (5) from Proposition 1, $e\leq[e]_{{\mathord{\uparrow}}}$. Hence, $x<[e]_{{\mathord{\uparrow}}}$. Since by hypothesis $e\in E_{\mathbb{F}}$ is an expression that evaluates on $\overline{\mathbb{R}}$ to a nonzero value, we have three cases:

$[e]_{{\mathord{\uparrow}}}\neq 0$ and $x\neq 0$:

In this case $x<[e]_{{\mathord{\uparrow}}}$ implies $x\prec[e]_{{\mathord{\uparrow}}}$.

$[e]_{{\mathord{\uparrow}}}=-0$:

In this case, $-{f_{\mathrm{min}}}<e<0$. Hence, $x<0$. Therefore $x\prec[e]_{{\mathord{\uparrow}}}$ holds.

$x=0$:

In this case it must be the case that $x$ is strictly smaller than $e$, since $e\in E_{\mathbb{F}}$ evaluates to a nonzero value. Therefore, $e>0$. Hence, by Definition 5, $[e]_{{\mathord{\uparrow}}}\geq{f_{\mathrm{min}}}$. Then $x\prec[e]_{{\mathord{\uparrow}}}$ holds.

In any case, $x\prec[e]_{{\mathord{\uparrow}}}$ holds. By Definition 10, we conclude that $x\prec\llbracket e\rrbracket_{{\mathord{\uparrow}}}$. By Definition 4, we have the following cases on $\llbracket e\rrbracket_{{\mathord{\uparrow}}}$:

$\llbracket e\rrbracket_{{\mathord{\uparrow}}}=-{f_{\mathrm{max}}}$:

In this case $\mathop{\mathrm{pred}}\nolimits(\llbracket e\rrbracket_{{\mathord{\uparrow}}})=-\infty$. Since $x\prec\llbracket e\rrbracket_{{\mathord{\uparrow}}}$, this implies that $x=-\infty$. Then $x\preccurlyeq\mathop{\mathrm{pred}}\nolimits(\llbracket e\rrbracket_{{\mathord{\uparrow}}})$, holds.

${f_{\mathrm{min}}}<\llbracket e\rrbracket_{{\mathord{\uparrow}}}\leq{f_{\mathrm{max}}}$ or $-{f_{\mathrm{max}}}<\llbracket e\rrbracket_{{\mathord{\uparrow}}}\leq-{f_{\mathrm{min}}}$:

In this case $\mathop{\mathrm{pred}}\nolimits(\llbracket e\rrbracket_{{\mathord{\uparrow}}})=\max\{\,y\in\mathbb{F}\mid y<\llbracket e\rrbracket_{{\mathord{\uparrow}}}\,\}$. Since $x<\llbracket e\rrbracket_{{\mathord{\uparrow}}}$, $x\in\{\,y\in\mathbb{F}\mid y<\llbracket e\rrbracket_{{\mathord{\uparrow}}}\,\}$. Hence, $x\preccurlyeq\mathop{\mathrm{pred}}\nolimits(\llbracket e\rrbracket_{{\mathord{\uparrow}}})$, holds.

$\llbracket e\rrbracket_{{\mathord{\uparrow}}}=+0$ or $\llbracket e\rrbracket_{{\mathord{\uparrow}}}=-0$:

In this case $\mathop{\mathrm{pred}}\nolimits(\llbracket e\rrbracket_{{\mathord{\uparrow}}})=-{f_{\mathrm{min}}}$. Since $x<\llbracket e\rrbracket_{{\mathord{\uparrow}}}$, this implies that $x\leq-{f_{\mathrm{min}}}$. Hence, $x\preccurlyeq\mathop{\mathrm{pred}}\nolimits(\llbracket e\rrbracket_{{\mathord{\uparrow}}})$, holds.

$\llbracket e\rrbracket_{{\mathord{\uparrow}}}={f_{\mathrm{min}}}$ :

In this case $\mathop{\mathrm{pred}}\nolimits(\llbracket e\rrbracket_{{\mathord{\uparrow}}})=+0$. Since $x<\llbracket e\rrbracket_{{\mathord{\uparrow}}}={f_{\mathrm{min}}}$, $x\preccurlyeq+0$. Hence, $x\preccurlyeq\mathop{\mathrm{pred}}\nolimits(\llbracket e\rrbracket_{{\mathord{\uparrow}}})$, holds.

$\llbracket e\rrbracket_{{\mathord{\downarrow}}}=+\infty$ :

In this case $\mathop{\mathrm{pred}}\nolimits(\llbracket e\rrbracket_{{\mathord{\uparrow}}})={f_{\mathrm{max}}}$. Since $x<\llbracket e\rrbracket_{{\mathord{\uparrow}}}=\infty$, $x\preccurlyeq{f_{\mathrm{max}}}$. Hence, $x\preccurlyeq\mathop{\mathrm{pred}}\nolimits(\llbracket e\rrbracket_{{\mathord{\uparrow}}})$, holds.

In order to prove (29) we first want to prove that $x\succcurlyeq[e]_{{\mathord{\uparrow}}}$. To this aim consider the following cases for $e$:

$e>{f_{\mathrm{max}}}$:

In this case $[e]_{{\mathord{\uparrow}}}=+\infty$. On the hand, $x\geq e>{f_{\mathrm{max}}}$. Since $x\in\mathbb{F}$ implies that $x=+\infty$. Hence $x\succcurlyeq[e]_{{\mathord{\uparrow}}}$.

$e\leq-{f_{\mathrm{min}}}$ or $0<e\leq{f_{\mathrm{max}}}$:

In this case $[e]_{{\mathord{\uparrow}}}=\min\{\,z\in\mathbb{F}\mid z\geq e\,\}$. Since $x\geq e$, $x\in\{\,z\in\mathbb{F}\mid z\geq e\,\}$. Hence, $x\geq[e]_{{\mathord{\uparrow}}}$, holds and also $x\succcurlyeq[e]_{{\mathord{\uparrow}}}$.

$-{f_{\mathrm{min}}}<e<0$:

In this case $[e]_{{\mathord{\uparrow}}}=-0$. Since $x\geq e$ and $x\in\mathbb{F}$, $x\succcurlyeq-0$, holds.

$e=-\infty$:

In this case $[e]_{{\mathord{\uparrow}}}=-\infty$ and $x\succcurlyeq-\infty$ holds.

Since by hypothesis $[e]_{{\mathord{\uparrow}}}=\llbracket e\rrbracket_{{\mathord{\uparrow}}}$, we can conclude that $x\succcurlyeq\llbracket e\rrbracket_{{\mathord{\uparrow}}}$ holds.

In order to prove (30) we first want to prove that $x\preccurlyeq[e]_{{\mathord{\downarrow}}}$. To this aim consider the following cases for $e$:

$e<-{f_{\mathrm{max}}}$:

In this case $[e]_{{\mathord{\downarrow}}}=-\infty$. On the hand, $x\leq e<-{f_{\mathrm{max}}}$. Since $x\in\mathbb{F}$ implies that $x=-\infty$. Hence $x\preccurlyeq[e]_{{\mathord{\downarrow}}}$.

$e\geq{f_{\mathrm{min}}}$ or $-{f_{\mathrm{max}}}\leq e<0$:

In this case $[e]_{{\mathord{\downarrow}}}=\max\{\,z\in\mathbb{F}\mid z\leq e\,\}$. Since $x\leq e$, $x\in\{\,z\in\mathbb{F}\mid z\leq e\,\}$. Hence, $x\leq[e]_{{\mathord{\downarrow}}}$, holds and also $x\preccurlyeq[e]_{{\mathord{\downarrow}}}$.

$0<e<{f_{\mathrm{min}}}$:

In this case $[e]_{{\mathord{\downarrow}}}=+0$. Since $x\leq e$ and $x\in\mathbb{F}$, $x\preccurlyeq+0$, holds.

$e=+\infty$:

In this case $[e]_{{\mathord{\downarrow}}}=\infty$ and $x\preccurlyeq+\infty$ holds.

Since by hypothesis $[e]_{{\mathord{\downarrow}}}=\llbracket e\rrbracket_{{\mathord{\downarrow}}}$, we can conclude that $x\preccurlyeq\llbracket e\rrbracket_{{\mathord{\downarrow}}}$ holds.

B.3 Proofs of Results in Section 4.3

Proof (of Theorem 4.3)

Given the constraint $x=y\boxslash_{S}z$ with $x\in X=[x_{\ell},x_{u}]$, $y\in Y=[y_{\ell},y_{u}]$ and $z\in Z=[z_{\ell},z_{u}]$, Algorithm 3 computes a new refining interval $X^{\prime}$ for variable $x$. Note that $X^{\prime}=[x^{\prime}_{\ell},x^{\prime}_{u}]\cap X$, which assures us that $X^{\prime}\subseteq X$.

As for the proof of Theorem A.5, it is easy to verify that $y_{L}$ and $w_{L}$ (resp., $y_{U}$ and $w_{U}$) computed using function $\tau$ of Figure 5, are the boundaries of $Y$ and $W$ upon which $x$ touches its minimum (resp., maximum). Moreover, remember that by Proposition 2, following the same reasoning of the proofs of the previous theorems, we can focus on finding a lower bound for $y_{L}\boxslash_{r_{\ell}}w_{L}$ and an upper bound for $y_{U}\boxslash_{r_{u}}w_{U}$.

We will now comment only on the most critical entries of function $\mathrm{dd}_{\ell}$ of Figure 6: let us briefly discuss the cases in which $y_{L}=-\infty$ and $w_{L}=\pm\infty$.

$w_{L}=-\infty.$

In this case, by function $\tau$ of Figure 5 (see the first three cases), we have $y_{L}=y_{u}=-\infty$, while either $w_{L}=w_{\ell}$ or $w_{L}=w_{u}$. Since by the IEEE 754 Standard IEEE-754-2008 dividing $\pm\infty$ by $\pm\infty$ is an invalid operation, we are left to consider the case $w_{L}=w_{\ell}$. In this case, recall that by the IEEE 754 Standard IEEE-754-2008 , dividing $-\infty$ by a finite negative number yields $+\infty$. Hence, we can conclude $x_{\ell}=+\infty$.

$w_{L}=+\infty.$

By function $\tau$ of Figure 5 (see the fourth and last case), we have $y_{L}=y_{\ell}=-\infty$, while $w_{L}=w_{\ell}=+\infty$. Hence, $x_{\ell}=-0$, since dividing a negative finite number by $+\infty$ gives $-0$.

A similar reasoning applies for the cases $y_{L}=+\infty$, $w_{L}=\pm\infty$. Dually, the only critical entries of function $\mathrm{dd}_{u}$ of Figure 6 are those in which $y_{U}=\pm\infty$ and $w_{U}=\pm\infty$ and can be handled analogously.

We are left to prove that $\forall X^{\prime\prime}\subset X,\exists r\in S,y\in Y,z\in Z\mathrel{:}y\boxslash_{r}z\not\in X^{\prime\prime}$. Let us focus on the lower bound $x^{+}_{\ell}$ proving that, if $[x^{+}_{\ell},x^{+}_{\ell}]\neq\varnothing$, then there exist $r\in S,y\in Y,z\in Z$ such that $y\boxslash_{r}z=x^{+}_{\ell}$. Consider the particular values $y_{L}$, $z_{\ell}=w_{L}$ and $r_{\ell}$ that correspond to $x^{+}_{\ell}$ in Algorithm 3, i.e. $y_{L}$ and $w_{L}$ and $r_{\ell}$ are such that $\mathrm{dd}_{\ell}(y_{L},w_{L},r_{\ell})=x^{+}_{\ell}$. By Algorithm 3, such $y_{L}$ and $w_{L}$ must exist. First consider the cases in which $y_{L}\not\in(\mathbb{R}_{-}\cup\mathbb{R}_{+})$ or $w_{L}\not\in(\mathbb{R}_{-}\cup\mathbb{R}_{+})$. A brute-force verification was successfully conducted, in this cases, to prove that $y_{L}\boxslash_{r_{\ell}}w_{L}=x^{+}_{\ell}$. For the cases in which $y_{L}\in(\mathbb{R}_{-}\cup\mathbb{R}_{+})$ and $w_{L}\in(\mathbb{R}_{-}\cup\mathbb{R}_{+})$ we have, by definition of $\mathrm{dd}_{\ell}$ of Figure 6, that $x^{+}_{\ell}=y_{L}\boxslash_{r_{\ell}}w_{L}$. Remember that, by Proposition 2, there exists $r\in S$ such that $y_{L}\boxslash_{\mathrm{r_{\ell}}}w_{L}=y_{L}\boxslash_{\mathrm{r}}w_{L}$. Since $y_{L}\in Y$ and $w_{L}\in Z$, we can conclude that $x^{+}_{\ell}\not\in X^{\prime\prime}$ implies that $y_{L}\boxslash_{\mathrm{r}}w_{L}\not\in X^{\prime\prime}$, for any $X^{\prime\prime}\subseteq X^{\prime}$. An analogous reasoning applies to $x^{-}_{\ell}$, to $x^{+}_{u}$ and $x^{-}_{u}$. This allows us to prove the optimality claim.

Proof (of Theorem 4.4)

Given the constraint $x=y\boxslash_{S}z$ with $x\in X=[x_{\ell},x_{u}]$, $y\in Y=[y_{\ell},y_{u}]$ and $z\in Z=[z_{\ell},z_{u}]$, Algorithm 4 computes a new, refining interval $Y^{\prime}$ for variable $y$. It returns either $Y^{\prime}\mathrel{\mathord{:}\mathord{=}}(Y\cap[y^{-}_{\ell},y^{-}_{u}])\biguplus(Y\cap[y^{+}_{\ell},y^{+}_{u}])$ or $Y^{\prime}=\varnothing$: hence, in both cases, we are sure that $Y^{\prime}\subseteq Y$.

By Proposition 2, we can focus on finding a lower bound for $y\in Y$ by exploiting the constraint $y\boxslash_{\bar{r}_{\ell}}z=x$ and an upper bound for $y$ by exploiting the constraint $y\boxslash_{\bar{r}_{u}}z=x$.

In order to compute correct bounds for $y$, Algorithm 5 first splits the interval of $z$ into the sign-homogeneous intervals $Z_{-}$ and $Z_{+}$, since knowing the sign of $z$ is crucial to determine correct bounds for $y$. Hence, for $W=Z_{-}$ (and, analogously, for $W=Z_{+}$), it calls function $\sigma$ of Figure 7 to determine the appropriate extrema of intervals $X$ and $W$ to be used to compute the new lower and upper bounds for $y$. As we did in the proof of Theorem A.4, it is easy to verify that $x_{L}$ and $w_{L}$ (resp., $x_{U}$ and $w_{U}$), computed using function $\sigma$ of Figure 7, are the boundaries of $X$ and $W$ upon which $y$ touches its minimum (resp., maximum). Functions $\mathrm{id}^{f}_{\ell}$ of Figure 8 and $\mathrm{id}^{f}_{u}$ of Figure 9 are then used to find the new bounds for $y$. The so obtained intervals for $y$ will be eventually joined using convex union to obtain the refining interval for $y$.

We will now prove the non-trivial parts of the definitions of functions $\mathrm{id}^{f}_{\ell}$ and $\mathrm{id}^{f}_{u}$. Concerning the case analysis of $\mathrm{id}^{f}_{\ell}$ (Fig 8) marked as $a_{4}$, the result changes depending on the selected rounding mode:

$\bar{r}_{\ell}={\mathord{\uparrow}}:$

we clearly must have $y=+\infty$, according to the IEEE 754 Standard IEEE-754-2008 ;

$\bar{r}_{\ell}={\mathord{\downarrow}}:$

it must be $y/w_{L}<-{f_{\mathrm{max}}}$ and thus, since $w_{L}$ is negative, $y>-{f_{\mathrm{max}}}\cdot w_{L}$ and, by (26) of Proposition 5, $y\succcurlyeq\mathop{\mathrm{succ}}\nolimits(-{f_{\mathrm{max}}}\boxdot_{\mathord{\downarrow}}w_{L})$.

$\bar{r}_{\ell}={\mathrm{n}}:$

since $\mathop{\mathrm{odd}}\nolimits({f_{\mathrm{max}}})$, for $w_{L}=-\infty$ we need $y$ to be greater than or equal to \bigl{(}-{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(-{f_{\mathrm{max}}})/2\bigr{)}\cdot w_{L}. If \llbracket\bigl{(}-{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(-{f_{\mathrm{max}}})/2\bigr{)}\cdot w_{L}\rrbracket_{{\mathord{\uparrow}}}=[\bigl{(}-{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(-{f_{\mathrm{max}}})/2\bigr{)}\cdot w_{L}]_{{\mathord{\uparrow}}}, by (29) of Proposition 5, we can conclude y\succcurlyeq\llbracket\bigl{(}-{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(-{f_{\mathrm{max}}})/2\bigr{)}\cdot w_{L}\rrbracket_{{\mathord{\uparrow}}}. On the other hand, if \llbracket\bigl{(}-{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(-{f_{\mathrm{max}}})/2\bigr{)}\cdot w_{L}\rrbracket_{{\mathord{\uparrow}}}\neq[\bigl{(}-{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(-{f_{\mathrm{max}}})/2\bigr{)}\cdot w_{L}]_{{\mathord{\uparrow}}}, then we can only apply (25) of Proposition 5, obtaining y\succcurlyeq\llbracket\bigl{(}-{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(-{f_{\mathrm{max}}})/2\bigr{)}\cdot w_{L}\rrbracket_{{\mathord{\downarrow}}}.

The case analysis of $\mathrm{id}^{f}_{\ell}$ (Fig 8) marked as $a_{5}$ can be explained as follows:

$\bar{r}_{\ell}={\mathord{\downarrow}}:$

we must have $y=+\infty$, according to the IEEE 754 Standard IEEE-754-2008 ;

$\bar{r}_{\ell}={\mathord{\uparrow}}:$

inequality $y/w_{L}>{f_{\mathrm{max}}}$ must hold and thus, since $w_{L}$ is positive, $y>{f_{\mathrm{max}}}\cdot w_{L}$ and, by (26) of Proposition 5, $y\succcurlyeq\mathop{\mathrm{succ}}\nolimits({f_{\mathrm{max}}}\boxdot_{\mathord{\downarrow}}w_{L})$.

$\bar{r}_{\ell}={\mathrm{n}}:$

since $\mathop{\mathrm{odd}}\nolimits({f_{\mathrm{max}}})$, for $x_{L}=+\infty$ we need $y$ to be greater than or equal to \bigl{(}{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2\bigl{)}\cdot w_{L}. If \llbracket\bigl{(}{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2\bigl{)}\cdot w_{L}\rrbracket_{{\mathord{\uparrow}}}=[\bigl{(}{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2\bigl{)}\cdot w_{L}]_{{\mathord{\uparrow}}}, by (29) of Proposition 5, we can conclude y\succcurlyeq\llbracket\bigl{(}{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2\bigl{)}\cdot w_{L}\rrbracket_{{\mathord{\uparrow}}}. On the other hand, if \llbracket\bigl{(}{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2\bigl{)}\cdot w_{L}\rrbracket_{{\mathord{\uparrow}}}\neq[\bigl{(}{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2\bigl{)}\cdot w_{L}]_{{\mathord{\uparrow}}} then, we can only apply (25) of Proposition 5, obtaining y\succcurlyeq\llbracket\bigl{(}{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2\bigl{)}\cdot w_{L}\rrbracket_{{\mathord{\downarrow}}}.

The explanation for the case analysis of $\mathrm{id}^{f}_{\ell}$ (Fig 8) marked as $a_{6}$ is the following:

$\bar{r}_{\ell}={\mathord{\uparrow}}:$

the lowest value of $y$ that yields $x_{L}=+0$ with $w_{L}\in\mathbb{R}_{-}$ is clearly $y=-0$;

$\bar{r}_{\ell}={\mathord{\downarrow}}:$

inequality $y/w_{L}<{f_{\mathrm{min}}}$ should hold and thus, since $w_{L}$ is negative, $y>{f_{\mathrm{min}}}\cdot w_{L}$ and, by (26) of Proposition 5, $y\succcurlyeq\mathop{\mathrm{succ}}\nolimits({f_{\mathrm{min}}}\boxdot_{\mathord{\downarrow}}w_{L})$.

$\bar{r}_{\ell}={\mathrm{n}}:$

since $\mathop{\mathrm{odd}}\nolimits({f_{\mathrm{min}}})$, for $x_{L}=+0$ we need $y$ to be greater than or equal to $({f_{\mathrm{min}}}\cdot w_{L})/2$. Since in this case $\llbracket({f_{\mathrm{min}}}\cdot w_{L})/2\rrbracket_{{\mathord{\uparrow}}}=[({f_{\mathrm{min}}}\cdot w_{L})/2]_{{\mathord{\uparrow}}}=({f_{\mathrm{min}}}\boxdot_{\mathord{\uparrow}}w_{L})/2$, by (29) of Proposition 5, we can conclude $y\succcurlyeq({f_{\mathrm{min}}}\boxdot_{\mathord{\uparrow}}w_{L})/2$.

Concerning the case analysis of $\mathrm{id}^{f}_{\ell}$ (Fig 8) marked as $a_{7}$, we must distinguish between the following cases:

$\bar{r}_{\ell}={\mathord{\downarrow}}:$

considering $x_{L}=-0$ and $w_{L}\in\mathbb{R}_{+}$, we clearly must have $y=-0$;

$\bar{r}_{\ell}={\mathord{\uparrow}}:$

it should be $y/w_{L}>-{f_{\mathrm{min}}}$ and thus, since $w_{L}$ is positive, $y>-{f_{\mathrm{min}}}\cdot w_{L}$ and, by (26) of Proposition 5, $y\succcurlyeq\mathop{\mathrm{succ}}\nolimits(-{f_{\mathrm{min}}}\boxdot_{\mathord{\downarrow}}w_{L})$.

$\bar{r}_{\ell}={\mathrm{n}}:$

since $\mathop{\mathrm{odd}}\nolimits({f_{\mathrm{min}}})$, for $x_{L}=-0$ we need $y$ be to greater than or equal to \bigl{(}-{f_{\mathrm{min}}}\cdot w_{L}\bigr{)}/2. Since in this case \llbracket\bigl{(}-{f_{\mathrm{min}}}\cdot w_{L}\bigr{)}/2\rrbracket_{{\mathord{\uparrow}}}=[\bigl{(}-{f_{\mathrm{min}}}\cdot w_{L}\bigr{)}/2]_{{\mathord{\uparrow}}}=(-{f_{\mathrm{min}}}\boxdot_{\mathord{\uparrow}}w_{L})/2, by (29) of Proposition 5, we can conclude $y\succcurlyeq(-{f_{\mathrm{min}}}\boxdot_{\mathord{\uparrow}}w_{L})/2$.

Similar arguments can be used to prove the case analyses of $\mathrm{id}^{f}_{u}$ of Fig 9 marked as $a_{9}$, $a_{10}$, $a_{11}$ and $a_{12}$.

We will now analyze the case analyses of $\mathrm{id}^{f}_{\ell}$ of Fig 8 marked as $a^{-}_{3}$ and $a^{+}_{3}$, and the ones of $\mathrm{id}^{f}_{u}$ of Fig 9 marked as $a^{-}_{8}$ and $a^{+}_{8}$. We can assume, of course, $X=[x_{\ell},x_{u}]$, $Y=[y_{\ell},y_{u}]$ and $Z=[w_{\ell},w_{u}]$, where $x_{\ell},x_{u},w_{\ell},w_{u}\in\mathbb{F}\cap\mathbb{R}$, $x_{\ell}\leq x_{u}$, $w_{\ell}\leq w_{u}$ and $\mathop{\mathrm{sgn}}\nolimits(w_{\ell})=\mathop{\mathrm{sgn}}\nolimits(w_{u})$. Exploiting $x\preccurlyeq y\boxslash z$ and $x\succcurlyeq y\boxslash z$, by Proposition 4, we have

[TABLE]

Since the case $z=0$ is handled separately by $\mathrm{id}^{f}_{\ell}$ of Fig 8 and by $\mathrm{id}^{f}_{u}$ of Fig 9, we can assume $z\neq 0$. Thanks to the split of $Z$ into a positive and a negative part, the sign of $z$ is determinate. In the following, we will prove the case analyses marked as $a_{3}^{+}$ and $a_{8}^{+}$, hence assuming $z>0$. From the previous case analysis we can derive

[TABLE]

Note that the members of the product are independent. Therefore, we can find the minimum of the product by minimizing each member of the product. Since we are analyzing the case in which $W=Z_{+}$, let $(x_{L},x_{U},w_{L},w_{U})$ as defined in function $\sigma$ of Figure 7, replacing the role of $y$ with $z$ and the role of $z$ with $x$. Hence, by Proposition 3 and the monotonicity of ‘$\mathop{\mathrm{pred}}\nolimits$’ and ‘$\mathop{\mathrm{succ}}\nolimits$’ we obtain

[TABLE]

We can now exploit Proposition 5 and obtain:

[TABLE]

Indeed, if $\bar{r}_{\ell}={\mathord{\uparrow}}$ and $x_{L}\neq 0$, then part (29) of Proposition 5 applies and we have $y\succcurlyeq x_{L}\boxdot_{\mathord{\uparrow}}w_{L}$. On the other hand, if $x_{L}=0$, since by hypothesis $z>0$ implies $w_{L}>0$, according to IEEE 754 (IEEE-754-2008, , Section 6.3), we have $x_{L}\boxdot_{\mathord{\uparrow}}w_{L}=\mathop{\mathrm{sgn}}\nolimits(x_{L})\cdot 0$ and, indeed, for each non-NaN, nonzero and finite $w\in\mathbb{F}\cap[+0,+\infty]$, $\mathop{\mathrm{sgn}}\nolimits(x_{L})\cdot 0$ is the least value for $y$ that satisfies $\mathop{\mathrm{sgn}}\nolimits(x_{L})\cdot 0=y\boxslash_{\mathord{\downarrow}}w$.

Analogously, if $\bar{r}_{\ell}={\mathord{\uparrow}}$ and $x_{L}\neq{f_{\mathrm{min}}}$, then Proposition 5 applies and we have \mathop{\mathrm{succ}}\nolimits\bigl{(}\mathop{\mathrm{pred}}\nolimits(x_{L})\boxdot_{\mathord{\downarrow}}w_{L}\bigr{)}. On the other hand, if $x_{L}={f_{\mathrm{min}}}$, in this case, \mathop{\mathrm{succ}}\nolimits\bigl{(}\mathop{\mathrm{pred}}\nolimits(x_{L})\boxdot_{\mathord{\downarrow}}w_{L}\bigr{)}={f_{\mathrm{min}}} which is consistent with the fact that, for each non-NaN, nonzero and finite $w\in\mathbb{F}\cap[+0,+\infty]$, ${f_{\mathrm{min}}}$ is the lowest value for $y$ that satisfies ${f_{\mathrm{min}}}=y\boxslash_{\mathord{\uparrow}}w$.

A symmetric argument justifies (53).

As before, we need to approximate the values of the expressions e^{+}_{\ell}=\bigl{(}x_{L}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(x_{L})/2\bigr{)}\cdot w_{L} and e^{+}_{u}=\bigl{(}x_{U}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(x_{U})/2\bigr{)}\cdot w_{U}. We leave this as an implementation choice, thus taking into account the case $\llbracket e^{+}_{\ell}\rrbracket_{{\mathord{\uparrow}}}=[e^{+}_{\ell}]_{{\mathord{\uparrow}}}$ and $\llbracket e^{+}_{u}\rrbracket_{{\mathord{\downarrow}}}=[e^{+}_{u}]_{{\mathord{\downarrow}}}$ as well as $\llbracket e^{+}_{\ell}\rrbracket_{{\mathord{\uparrow}}}>[e^{+}_{\ell}]_{{\mathord{\uparrow}}}$ and $\llbracket e^{+}_{u}\rrbracket_{{\mathord{\downarrow}}}<[e^{+}_{u}]_{{\mathord{\downarrow}}}$. Therefore, when $\llbracket e^{+}_{u}\rrbracket_{{\mathord{\downarrow}}}<[e^{+}_{u}]_{{\mathord{\downarrow}}}$ by (51) and (27) of Proposition 5 we obtain $y\preccurlyeq\llbracket e^{+}_{u}\rrbracket_{{\mathord{\uparrow}}}$, while, when $\llbracket e^{+}_{\ell}\rrbracket_{{\mathord{\downarrow}}}>[e^{+}_{\ell}]_{{\mathord{\downarrow}}}$ by (51) and (25) of Proposition 5 we obtain $y\succcurlyeq\llbracket e^{+}_{\ell}\rrbracket_{{\mathord{\downarrow}}}$.

Thus, for the case in which $\bar{r}_{\ell}={\mathrm{n}}$, since $e^{+}_{u}\neq 0$ and $e^{+}_{\ell}\neq 0$, by Proposition 5, we have

[TABLE]

An analogous reasoning, but with $z<0$, allows us to obtain the case analyses marked as $a_{3}^{-}$ and $a_{8}^{-}$.

Proof (of Theorem 4.5)

Given the constraint $x=y\boxslash_{S}z$ with $x\in X=[x_{\ell},x_{u}]$, $y\in Y=[y_{\ell},y_{u}]$ and $z\in Z=[z_{\ell},z_{u}]$, Algorithm 5 finds a new, refined interval $Z^{\prime}$ for variable $z$.

Since it assigns either $Z^{\prime}\mathrel{\mathord{:}\mathord{=}}(Z\cap[z^{-}_{\ell},z^{-}_{u}])\biguplus(Z\cap[z^{+}_{\ell},z^{+}_{u}])$ or $Z^{\prime}=\varnothing$, in both cases we are sure that $Z^{\prime}\subseteq Z$. By Proposition 2, as in the previous proofs, we can focus on finding a lower bound for $z\in Z$ by exploiting the constraint $y\boxslash_{\bar{r}_{\ell}}z=x$ and an upper bound for $z$ by exploiting the constraint $y\boxslash_{\bar{r}_{u}}z=x$.

We first need to split interval $X$ into the sign-homogeneous intervals $X_{-}$ and $X_{+}$, because knowing the sign of $x$ is crucial for determining correct bounds for $z$. Hence, for $V=X_{-}$ (and, analogously, for $V=X_{+}$) function $\tau$ of Figure 5 determines the appropriate interval extrema of $Y$ and $V$ to be used to compute the new lower and upper bounds for $z$. As in the previous proofs (see, for example, proof of Theorem A.5), it is easy to verify that $y_{L}$ and $v_{L}$ (resp., $y_{U}$ and $v_{U}$) computed using function $\tau$ of Figure 5 are the boundaries of $Y$ and $V$ upon which $z$ touches its minimum (resp., maximum). Functions $\mathrm{id}^{s}_{\ell}$ of Figure 10 and $\mathrm{id}^{s}_{u}$ of Figure 11 are then used to find the new bounds for $z$. The so obtained intervals for $z$ will be then joined with convex union in order to obtain the refining interval for $z$.

We will prove the most important parts of the definitions of $\mathrm{id}^{s}_{\ell}$ (Figure 10) and $\mathrm{id}^{s}_{u}$ (Figure 11) only, starting with the case analysis marked as $a_{4}$. Depending on the rounding mode in effect, the following arguments are given:

$\bar{r}_{\ell}={\mathord{\downarrow}}:$

in this case, the only possible way to obtain $-0$ as the result of the division is having $z=+\infty$ (with $y\in\mathbb{R}_{-}$);

$\bar{r}_{\ell}={\mathord{\uparrow}}:$

it should be $y_{L}/z>-{f_{\mathrm{min}}}$ and thus, since $y_{L}$ and $x_{L}$ are negative, we can conclude that $z$ is positive. Thus, $y_{L}>-{f_{\mathrm{min}}}\cdot z$ implies $y_{L}/-{f_{\mathrm{min}}}<z$, and by (26) of Proposition 5, $z\succcurlyeq\mathop{\mathrm{succ}}\nolimits(z_{L}\boxslash_{\mathord{\downarrow}}-{f_{\mathrm{min}}})$.

$\bar{r}_{\ell}={\mathrm{n}}:$

since $\mathop{\mathrm{odd}}\nolimits(-{f_{\mathrm{min}}})$, for $v_{L}=-0$ we need $y_{L}/z\geq(-{f_{\mathrm{min}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(-{f_{\mathrm{min}}})/2)=(-{f_{\mathrm{min}}}+{f_{\mathrm{min}}}/2)=-{f_{\mathrm{min}}}/2$. As before, since $y_{L}$ and $v_{L}$ are negative, we can conclude that $z$ is positive: hence $y_{L}\geq(-{f_{\mathrm{min}}}/2)\cdot z$. Therefore, $z\geq y_{L}/(-{f_{\mathrm{min}}}/2)=z\geq(y_{L}/-{f_{\mathrm{min}}})\cdot 2$. Since in this case $\llbracket(y_{L}/-{f_{\mathrm{min}}})\cdot 2\rrbracket_{{\mathord{\uparrow}}}=[(y_{L}/-{f_{\mathrm{min}}})\cdot 2]_{{\mathord{\uparrow}}}=(y_{L}\boxslash_{\mathord{\uparrow}}-{f_{\mathrm{min}}})\cdot 2$, by (29) of Proposition 5, we can conclude $y\succcurlyeq(y_{L}\boxslash_{\mathord{\uparrow}}-{f_{\mathrm{min}}})\cdot 2$.

As for the case analysis of $\mathrm{id}^{s}_{\ell}$ (Figure 10) marked as $a_{5}$, we must distinguish between the following cases:

$\bar{r}_{\ell}={\mathord{\uparrow}}:$

we must have $z=+\infty$ in order to obtain $x=+0$;

$\bar{r}_{\ell}={\mathord{\downarrow}}:$

inequality $y_{L}/z<{f_{\mathrm{min}}}$ must hold and thus, since positive $y_{L}$ and $v_{L}$ imply a positive $z$, $z>y_{L}/{f_{\mathrm{min}}}$ and, by (26) of Proposition 5, $z\succcurlyeq\mathop{\mathrm{succ}}\nolimits(y_{L}\boxslash_{\mathord{\downarrow}}{f_{\mathrm{min}}})$.

$\bar{r}_{\ell}={\mathrm{n}}:$

since $\mathop{\mathrm{odd}}\nolimits({f_{\mathrm{min}}})$, for $v_{L}=+0$ we need $y_{L}/z\leq{f_{\mathrm{min}}}/2$. As $z$ is positive in this case, $(y_{L}/{f_{\mathrm{min}}})\cdot 2\leq z$. Since $\llbracket(y_{L}/{f_{\mathrm{min}}})\cdot 2\rrbracket_{{\mathord{\uparrow}}}=[(y_{L}/{f_{\mathrm{min}}})\cdot 2]_{{\mathord{\uparrow}}}=(y_{L}\boxslash_{\mathord{\uparrow}}{f_{\mathrm{min}}})\cdot 2$, by (29) of Proposition 5, we can conclude $y\succcurlyeq(y_{L}\boxslash_{\mathord{\uparrow}}{f_{\mathrm{min}}})\cdot 2$.

Concerning the case analysis of $\mathrm{id}^{s}_{\ell}$ (Fig 10) marked as $a_{6}$, we must distinguish between the following cases:

$\bar{r}_{\ell}={\mathord{\downarrow}}:$

the lowest value of $z$ that gives $x=+\infty$ with $y\in\mathbb{R}_{-}$ is $z=-0$;

$\bar{r}_{\ell}={\mathord{\uparrow}}:$

inequality $y_{L}/z>{f_{\mathrm{max}}}$ must hold; since $y_{L}$ is negative and $v_{L}$ is positive, $z$ must be negative, and therefore $y_{L}<{f_{\mathrm{max}}}\cdot z$. Hence, $y_{L}/{f_{\mathrm{max}}}<z$. By (26) of Proposition 5, we obtain $z\succcurlyeq\mathop{\mathrm{succ}}\nolimits(y_{L}\boxslash_{\mathord{\downarrow}}{f_{\mathrm{max}}})$.

$\bar{r}_{\ell}={\mathrm{n}}:$

since $\mathop{\mathrm{odd}}\nolimits({f_{\mathrm{max}}})$, for $v_{L}=+\infty$ we need $y_{L}/z\geq({f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2)$. As before, since $w_{L}$ is negative and $v_{L}$ is positive, we can conclude that $z$ is negative, and, therefore, $y_{L}\leq({f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2)\cdot z$ holds. As a consequence, $y_{L}/({f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2)\leq z$. If $\llbracket y_{L}/({f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2)\rrbracket_{{\mathord{\uparrow}}}=[y_{L}/({f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2)]_{{\mathord{\uparrow}}}$, by (29) of Proposition 5, we can conclude $z\succcurlyeq\llbracket y_{L}/({f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2)\rrbracket_{{\mathord{\uparrow}}}$. On the other hand, if $\llbracket y_{L}/({f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2)\rrbracket_{{\mathord{\uparrow}}}\neq[y_{L}/({f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2)]_{{\mathord{\uparrow}}}$ then, we can only apply (25) of Proposition 5, obtaining $z\succcurlyeq\llbracket y_{L}/({f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2)\rrbracket_{{\mathord{\downarrow}}}$.

Regarding the case analysis of $\mathrm{id}^{s}_{\ell}$ (Fig 10) marked as $a_{7}$, we have the following cases:

$\bar{r}_{\ell}={\mathord{\uparrow}}:$

the lowest value of $z$ that yields $x=-\infty$ with $y\in\mathbb{R}_{+}$ is $z=-0$;

$\bar{r}_{\ell}={\mathord{\downarrow}}:$

inequality $y_{L}/z<-{f_{\mathrm{max}}}$ must hold and thus, since a positive $y_{L}$ and a negative $v_{L}$ imply that the sign of $z$ is negative, $y_{L}>-{f_{\mathrm{max}}}\cdot z$. Hence, $y_{L}/-{f_{\mathrm{max}}}<z$. By (26) of Proposition 5, $z\succcurlyeq\mathop{\mathrm{succ}}\nolimits(y_{L}\boxslash_{\mathord{\downarrow}}-{f_{\mathrm{max}}})$.

$\bar{r}_{\ell}={\mathrm{n}}:$

since $\mathop{\mathrm{odd}}\nolimits(-{f_{\mathrm{max}}})$, for $v_{L}=-\infty$ we need $y_{L}/z\leq-{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(-{f_{\mathrm{max}}})/2$. Since $z$ in this case is negative, we obtain the inequality $z\geq y_{L}/(-{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(-{f_{\mathrm{max}}})/2)$. If \llbracket y_{L}/\bigl{(}-{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(-{f_{\mathrm{max}}})/2\bigr{)}\rrbracket_{{\mathord{\uparrow}}}=[y_{L}/\bigl{(}-{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(-{f_{\mathrm{max}}})/2\bigr{)}]_{{\mathord{\uparrow}}}, by (29) of Proposition 5, we can conclude y\succcurlyeq\llbracket y_{L}/\bigl{(}-{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(-{f_{\mathrm{max}}})/2\bigr{)}\rrbracket_{{\mathord{\uparrow}}}. On the other hand, if \llbracket y_{L}/\bigl{(}-{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(-{f_{\mathrm{max}}})/2\bigr{)}\rrbracket_{{\mathord{\uparrow}}}\neq[y_{L}/\bigl{(}-{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(-{f_{\mathrm{max}}})/2\bigr{)}]_{{\mathord{\uparrow}}}, then we can only apply (25) of Proposition 5, obtaining y\succcurlyeq\llbracket y_{L}/\bigl{(}-{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(-{f_{\mathrm{max}}})/2\bigr{)}\rrbracket_{{\mathord{\downarrow}}}.

Similar arguments can be used to prove the case analyses of function $\mathrm{id}^{s}_{u}$ of Figure 11 marked as $a_{9}$, $a_{10}$, $a_{11}$ and $a_{12}$.

We will now analyze the case analyses of $\mathrm{id}^{s}_{\ell}$ of Figure 10 marked as $a^{-}_{3}$ and $a^{+}_{3}$, and the ones of $\mathrm{id}^{s}_{u}$ of Figure 9 marked as $a^{-}_{8}$ and $a^{+}_{8}$. In this proof, we can assume $y_{L},v_{L}\in\mathbb{R}_{-}\cup\mathbb{R}_{+}$, $y_{U},v_{U}\in\mathbb{R}_{-}\cup\mathbb{R}_{+}$ and $\mathop{\mathrm{sgn}}\nolimits(v_{L})=\mathop{\mathrm{sgn}}\nolimits(v_{U})$. First, note that the argument that leads to (48) and (49) starting from $x\preccurlyeq y\boxslash z$ and $x\succcurlyeq y\boxslash z$ is in common with the proof of Theorem 4.4.

Provided that interval $X$ is split into intervals $X_{+}$ and $X_{-}$, it is worth discussing the reasons why it is not necessary to partition also $Y$ directly in Algorithm 5. Assume $Y=[-a,b]$ with $a,b>0$ and consider the partition of $Y$ into two sign-homogeneus intervals $Y\cap[-\infty,-0]$ and $Y\cap[+0,+\infty]$, as usual. Note that the values $-0\in Y\cap[-\infty,-0]=[-a,-0]$ and the values $+0\in Y\cap[+0,+\infty]=[+0,b]$ can never be the boundaries of $Y$ upon which $z$ touches its minimum (resp., maximum). This is because $y$ will be the numerator of fractions (see expressions (56) and (57)). Moreover, by the definition of functions $\mathrm{id}^{s}_{\ell}$ of Fig 10 and $\mathrm{id}^{s}_{u}$ of Fig 11, it easy to verify that the partition of $Y$ would not prevent the interval computed for $y$ from being equal to the empty set. That is, if $\mathrm{id}^{s}_{\ell}(y_{L},v_{L},\bar{r}_{\ell})=\text{unsat.}$ or $\mathrm{id}^{s}_{u}(y_{U},v_{U},\bar{r}_{u})=\text{unsat.}$, then partitioning also $Y$ into sign-homogeneus intervals and then applying the procedure of Algorithm 5 to the two distinct intervals results again into an empty refining interval for $z$.

Hence, to improve efficiency, Algorithm 5 does not split interval $Y$ into sign-homogeneous intervals. However, in this proof it is necessary to partition $Y$ into intervals $Y_{-}$ and $Y_{+}$ in order to determine the correct formulas for lower and upper bounds for $z$. In the following, for the sake of simplicity, we will analyze the special case $X_{+}$ and $Y=Y_{+}$, so that $Y$ does not need to be split because it is already a sign-homogeneous interval. The remaining cases in which $Y$ is sign-homogeneous as well as those in which it is not can be derived analogously. To sum up, in this case we assume $x\geq 0$ and $y\geq 0$, and therefore $z>0$.

Now, we need to prove the cases marked as $a_{3}^{+}$ and $a_{8}^{+}$. The case analysis of (46) and (47) yields (48) and (49). Remember that the case $x=\pm 0$ is handled separately by functions $\mathrm{id}^{s}_{\ell}$ of Figure 10 and $\mathrm{id}^{s}_{u}$ of Figure 11, hence assuming $x>0$, we obtain

[TABLE]

Since the members of the divisions are independent, we can find the minimum of said divisions by minimizing each one of their members. Let $(y_{L},y_{U},v_{L},v_{U})$ be as returned by function $\tau$ of Figure 5. By Proposition 3 and the monotonicity of ‘$\mathop{\mathrm{pred}}\nolimits$’ and ‘$\mathop{\mathrm{succ}}\nolimits$’ we obtain

[TABLE]

We can now exploit Proposition 5 and obtain:

[TABLE]

Since $y_{L}\neq 0$, then $y_{L}/\mathop{\mathrm{succ}}\nolimits(v_{L})\neq 0$. Hence, Proposition 5 applies and we have $z\succcurlyeq y_{L}\boxslash_{\mathord{\uparrow}}v_{L}$ if $\bar{r}_{\ell}={\mathord{\uparrow}}$ and z\succcurlyeq\mathop{\mathrm{succ}}\nolimits\bigl{(}y_{L}/\mathop{\mathrm{succ}}\nolimits(v_{L})\bigr{)} if $\bar{r}_{\ell}={\mathord{\downarrow}}$ and $v_{L}\neq-{f_{\mathrm{min}}}$. Analogously, since $y_{U}\neq 0$, then $y_{U}/\mathop{\mathrm{pred}}\nolimits(v_{L})\neq 0$. Hence, by Proposition 5 we obtain (61).

Note that, since division by zero is not defined on real numbers, we had to separately address the case $\bar{r}_{u}={\mathord{\uparrow}}$ and $x={f_{\mathrm{min}}}$ in (56), and the case $\bar{r}_{\ell}={\mathord{\downarrow}}$ and $x=-{f_{\mathrm{min}}}$ in (57). Division by zero is, however, defined on IEEE 754 floating-point numbers. Indeed, if we evaluate the second case of (60) with $v_{L}=-{f_{\mathrm{min}}}$, we obtain $\mathop{\mathrm{succ}}\nolimits(y_{L}\boxslash_{\mathord{\downarrow}}\mathop{\mathrm{succ}}\nolimits(-{f_{\mathrm{min}}}))=-{f_{\mathrm{max}}}$, which happens to be the correct value for $z^{\prime}_{\ell}$, provided $y_{L}>0$. The same happens for (61). Therefore, there is no need for a separate treatment when variable $x$ takes the values $\pm{f_{\mathrm{min}}}$.

As before, we need to approximate the values of the expressions e^{+}_{u}\mathrel{\buildrel\mathrm{def}\over{=}}y_{U}/\bigl{(}v_{U}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(v_{U})/2\bigr{)} and e^{+}_{\ell}\mathrel{\buildrel\mathrm{def}\over{=}}y_{L}/\bigl{(}v_{L}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits(v_{L})/2\bigr{)}. Thus, when $\llbracket e^{+}_{u}\rrbracket_{{\mathord{\downarrow}}}<[e^{+}_{u}]_{{\mathord{\downarrow}}}$ by (51) and (27) of Proposition 5 we obtain $y\preccurlyeq\llbracket e^{+}_{u}\rrbracket_{{\mathord{\uparrow}}}$, while, when $\llbracket e^{+}_{\ell}\rrbracket_{{\mathord{\downarrow}}}>[e^{+}_{\ell}]_{{\mathord{\downarrow}}}$ by (51) and (25) of Proposition 5 we obtain $y\succcurlyeq\llbracket e^{+}_{\ell}\rrbracket_{{\mathord{\downarrow}}}$. Thus, for the case where $\bar{r}_{\ell}={\mathrm{n}}$, since $e^{+}_{u}\neq 0$ and $e^{+}_{\ell}\neq 0$, by Proposition 5, we have

[TABLE]

An analogous reasoning allows us to prove the case analyses marked as $a_{3}^{-}$ and $a_{8}^{-}$.

Proof (of Theorem A.4)

Given the constraint $x=y\boxdot_{S}z$ with $x\in X=[x_{\ell},x_{u}]$, $y\in Y=[y_{\ell},y_{u}]$ and $z\in Z=[z_{\ell},z_{u}]$, then $X^{\prime}=[x^{\prime}_{\ell},x^{\prime}_{u}]\cap X$. Hence, we are sure that $X^{\prime}\subseteq X$.

It should be immediate to verify that function $\sigma$ of Figure 7, related to the case $\mathop{\mathrm{sgn}}\nolimits(y_{\ell})=\mathop{\mathrm{sgn}}\nolimits(y_{u})$, chooses the appropriate interval extrema $y_{L},y_{U},z_{L},z_{U}$, necessary for computing bounds for $x$. Indeed, note that such choice is completely driven by the sign of the resulting product. Analogously, the correct interval extrema $y_{L},y_{U},z_{L},z_{U}$ related to the case $\mathop{\mathrm{sgn}}\nolimits(z_{\ell})=\mathop{\mathrm{sgn}}\nolimits(z_{u})$ can be determined by applying function $\sigma$ of Figure 7, but swapping the role of $y$ and $z$. Hence, if the sign of $y$ or of $z$ is constant (see the second part of Algorithm 9) function $\sigma$ of Figure 7 finds the appropriate extrema for $y$ and $z$ to compute the bound for $x$.

Concerning the cases $\mathop{\mathrm{sgn}}\nolimits(y_{\ell})=\mathop{\mathrm{sgn}}\nolimits(z_{\ell})=-1$ and $\mathop{\mathrm{sgn}}\nolimits(y_{u})=\mathop{\mathrm{sgn}}\nolimits(z_{u})=1$ (first part of Algorithm 9), note that we have only two possibilities for the interval extrema $y_{L}$ and $z_{L}$, that are $y_{\ell}$ and $z_{u}$ or $y_{u}$ and $z_{\ell}$. Since the product of $y_{L}$ and $z_{L}$ will have a negative sign in both cases, the right extrema for determining the lower bound $x^{\prime}_{\ell}$ have to be chosen by selecting the smallest product of $y_{L}$ and $z_{L}$. Analogously, for $y_{U}$ and $z_{U}$ there are two possibilities: $y_{\ell}$ and $z_{\ell}$ or $y_{u}$ and $z_{u}$. Since the product of $y_{U}$ and $z_{U}$ will have a positive sign in both cases, the appropriate extrema for determining the upper bound $x^{\prime}_{u}$ have to be chosen as the biggest product of $y_{U}$ and $z_{U}$.

Remember that by Proposition 2, following the same reasoning as in the previous proofs, it suffices to find a lower bound for $y_{L}\boxdot_{r_{\ell}}z_{L}$ and an upper bound for $y_{U}\boxdot_{r_{u}}z_{U}$.

We now comment on some critical case analyses of function $\mathrm{dm}_{\ell}$ of Figure 17. Consider, for example, when $y_{L}=\pm\infty$ and $z_{L}=\pm 0$. In particular, we analyze the case in which $y_{L}=-\infty$ and $z_{L}=\pm 0$. Note that $y_{L}=-\infty$ implies $y_{\ell}=-\infty$. Assume, first, that $z_{L}=+0$. Recall that by the IEEE 754 Standard IEEE-754-2008 $\pm\infty\boxdot\pm 0$ is an invalid operation. However, since $y_{\ell}=-\infty$, we have two cases:

$y_{u}\geq-{f_{\mathrm{max}}}:$

note that, in this case, $-{f_{\mathrm{max}}}\boxdot+0=-0$;

$y_{u}=-\infty:$

in this case, $z_{L}$ must correspond to $z_{u}$ (see the last three cases of function $\sigma$). Since $-\infty\boxdot z$ for $z<0$ results in $+\infty$, we can conclude that $-0$ is a correct lower bound for $x$.

A similar reasoning applies for the cases $y_{L}=+\infty$, $z_{L}=\pm 0$. Dually, the only critical entries of function $\mathrm{dm}_{u}$ of Figure 17 are those in which $y_{U}=\pm\infty$ and $z_{U}=\pm 0$. In these cases we can reason in a similar way, too.

We are left to prove that $\forall X^{\prime\prime}\subset X\mathrel{:}\exists r\in S,y\in Y,z\in Z\mathrel{.}y\boxdot_{r}z\not\in X^{\prime\prime}$. Let us focus on the lower bound $x^{\prime}_{\ell}$, proving that there exist values $r\in S,y\in Y,z\in Z$ such that $y\boxdot_{r}z=x^{\prime}_{\ell}$. Consider the particular values of $y_{L}$, $z_{L}$ and $r_{\ell}$ that correspond to the value of $x^{\prime}_{\ell}$ chosen by Algorithm 9, that is $y_{L}$, $z_{L}$ and $r_{\ell}$ are such that $\mathrm{dm}_{\ell}(y_{L},z_{L},r_{\ell})=x^{\prime}_{\ell}$. By Algorithm 9, such values of $y_{L}$ and $z_{L}$ must exist. First, consider the cases in which $y_{L}\not\in(\mathbb{R}_{-}\cup\mathbb{R}_{+})$ or $z_{L}\not\in(\mathbb{R}_{-}\cup\mathbb{R}_{+})$. In these cases, a brute-force verification was successfully conducted to verify that $y\boxdot_{r_{\ell}}z=x^{\prime}_{\ell}$. For the cases in which $y_{L}\in(\mathbb{R}_{-}\cup\mathbb{R}_{+})$ and $z_{L}\in(\mathbb{R}_{-}\cup\mathbb{R}_{+})$ we have, by definition of $\mathrm{dm}_{\ell}$ of Figure 17, that $x^{\prime}_{\ell}=y_{L}\boxdot_{r_{\ell}}z_{L}$. Remember that, by Proposition 2, there exist $r^{\prime}\in S$ such that $y_{L}\boxdot_{\mathrm{r_{\ell}}}z_{L}=y_{L}\boxdot_{\mathrm{r^{\prime}}}z_{L}$. Since $y_{L}\in Y$ and $z_{L}\in Z$, we can conclude that $x^{\prime}_{\ell}\not\in X^{\prime\prime}$ implies that $y^{\prime}_{L}\boxdot_{\mathrm{r^{\prime}}}z_{L}\not\in X^{\prime\prime}$. An analogous reasoning allows us to conclude that $\exists r\in S$ for which the following holds: $x^{\prime}_{u}\not\in X^{\prime\prime}$ implies $y_{U}\boxdot_{\mathrm{r}}z_{U}\not\in X^{\prime\prime}$.

Proof (of Theorem A.5)

Given the constraint $x=y\boxdot_{S}z$ with $x\in X=[x_{\ell},x_{u}]$, $y\in Y=[y_{\ell},y_{u}]$ and $z\in Z=[z_{\ell},z_{u}]$, Algorithm 10 computes $Y^{\prime}$, a new and refined interval for variable $y$.

First, note that either $Y^{\prime}\mathrel{\mathord{:}\mathord{=}}(Y\cap[y^{-}_{\ell},y^{-}_{u}])\biguplus(Y\cap[y^{+}_{\ell},y^{+}_{u}])$ or $Y^{\prime}=\varnothing$, hence, in both cases, we are sure that $Y^{\prime}\subseteq Y$ holds.

By Proposition 2, we can focus on finding a lower bound for $y\in Y$ by exploiting the constraint $y\boxdot_{\bar{r}_{\ell}}z=x$ and an upper bound for $y\in Y$ by exploiting the constraint $y\boxdot_{\bar{r}_{u}}z=x$.

Now, in order to compute correct bounds for $y$, we first need to split the interval of $z$ into the sign-homogeneous intervals $Z_{-}$ and $Z_{+}$, because it is crucial to be sure of the sign of $z$. As a consequence, for $W=Y_{-}$ (and, analogously, for $W=Y_{+}$) function $\tau$ of Figure 5 picks the appropriate interval extrema of $X$ and $W$ to be used to compute the new lower and upper bounds for $y$. It is easy to verify that the values of $x_{L}$ and $w_{L}$ (resp., $x_{U}$ and $w_{U}$) computed using function $\tau$ of Figure 5 are the boundaries of $X$ and $W$ upon which $y$ touches its minimum (resp., maximum). Functions $\mathrm{im}_{\ell}$ of Figure 18 and $\mathrm{im}_{u}$ of Figure 19 are then employed to find the new bounds for $y$. The so obtained intervals for $y$ are then joined using convex union between intervals, in order to obtain the refining interval for $y$.

Observe that functions $\mathrm{im}_{\ell}$ of Fig 18 and $\mathrm{im}_{u}$ of Fig 19 are dual to each other: every row/column of one table can be found in the other table reversed and changed of sign. This is due to the fact that, for each $r\in R$ and each $D\subseteq\mathbb{F}\times\mathbb{F}$, we have

[TABLE]

Concerning the case analysis of $\mathrm{im}_{\ell}$ marked as $a_{4}$ of Fig 18, we must consider the following cases:

$\bar{r}_{\ell}={\mathord{\uparrow}}:$

we clearly must have $y=+\infty$ in this case;

$\bar{r}_{\ell}={\mathord{\downarrow}}:$

inequality $y\cdot w_{L}<-{f_{\mathrm{max}}}$ must hold and thus, since $w_{L}$ is negative, $y>-{f_{\mathrm{max}}}/w_{L}$ and, by (26) of Proposition 5, $y\succcurlyeq\mathop{\mathrm{succ}}\nolimits(-{f_{\mathrm{max}}}\boxslash_{\mathord{\downarrow}}w_{L})$.

$\bar{r}_{\ell}={\mathrm{n}}:$

since $\mathop{\mathrm{odd}}\nolimits({f_{\mathrm{max}}})$, for $x_{L}=-\infty$ we need $y$ to be greater or equal than \bigl{(}-{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(-{f_{\mathrm{max}}})/2\bigr{)}/w_{L}. If \llbracket\bigl{(}-{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(-{f_{\mathrm{max}}})/2\bigr{)}/w_{L}\rrbracket_{{\mathord{\uparrow}}}=[-{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(-{f_{\mathrm{max}}})/2\bigr{)}/w_{L}]_{{\mathord{\uparrow}}}, by (29) of Proposition 5, we can conclude y\succcurlyeq\llbracket\bigl{(}-{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(-{f_{\mathrm{max}}})/2\bigr{)}/w_{L}\rrbracket_{{\mathord{\uparrow}}}. On the other hand, if \llbracket\bigl{(}-{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(-{f_{\mathrm{max}}})/2\bigr{)}/w_{L}\rrbracket_{{\mathord{\uparrow}}}\neq[-{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(-{f_{\mathrm{max}}})/2\bigr{)}/w_{L}]_{{\mathord{\uparrow}}}, then we can only apply (25) of Proposition 5, obtaining y\succcurlyeq\llbracket\bigl{(}-{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}-}}\nolimits(-{f_{\mathrm{max}}})/2\bigr{)}/w_{L}\rrbracket_{{\mathord{\downarrow}}}.

Regarding the case analysis of $\mathrm{im}_{\ell}$ marked a $a_{5}$ of Fig 18, we have the following cases:

$\bar{r}_{\ell}={\mathord{\downarrow}}:$

in this case, we must have $y=-0$;

$\bar{r}_{\ell}={\mathord{\uparrow}}:$

inequality $y\cdot w_{L}>-{f_{\mathrm{min}}}$ must hold and thus, since $w_{L}$ is positive, $y>-{f_{\mathrm{min}}}/w_{L}$ and, by (28) of Proposition 5, $y\succcurlyeq\mathop{\mathrm{succ}}\nolimits(-{f_{\mathrm{min}}}\boxslash_{\mathord{\downarrow}}w_{L})$.

$\bar{r}_{\ell}={\mathrm{n}}:$

since $\mathop{\mathrm{odd}}\nolimits({f_{\mathrm{min}}})$, for $x_{L}=-0$ we need $y$ to be greater or equal than $-{f_{\mathrm{min}}}/(2\cdot w_{L})$. Since in this case $\llbracket-{f_{\mathrm{min}}}/(2\cdot w_{L})\rrbracket_{{\mathord{\uparrow}}}=[-{f_{\mathrm{min}}}/(2\cdot w_{L})]_{{\mathord{\uparrow}}}=(-{f_{\mathrm{min}}})\boxslash_{\mathord{\uparrow}}(2\cdot w_{L})$, by (29) of Proposition 5, we can conclude $y\succcurlyeq-{f_{\mathrm{min}}}\boxslash_{\mathord{\uparrow}}(2\cdot w_{L})$.

As for the case analysis of $\mathrm{im}_{\ell}$ marked as $a_{6}$ of Figure 18, the following cases must be studied:

$\bar{r}_{\ell}={\mathord{\uparrow}}:$

we must have $y=+0$ in this case;

$\bar{r}_{\ell}={\mathord{\downarrow}}:$

it should be $y\cdot w_{L}<{f_{\mathrm{min}}}$ and thus, since $w_{L}$ is negative, $y>{f_{\mathrm{min}}}/w_{L}$ and, by (28) of Proposition 5, $y\succcurlyeq\mathop{\mathrm{succ}}\nolimits(-{f_{\mathrm{min}}}\boxslash_{\mathord{\downarrow}}w_{L})$.

$\bar{r}_{\ell}={\mathrm{n}}:$

since $\mathop{\mathrm{odd}}\nolimits({f_{\mathrm{min}}})$, for $x_{L}=-0$ we need $y$ to be greater than or equal to \bigl{(}{f_{\mathrm{min}}}/(2\cdot w_{L})\bigr{)}. Since in this case $\llbracket{f_{\mathrm{min}}}/(2\cdot w_{L})\rrbracket_{{\mathord{\uparrow}}}=[{f_{\mathrm{min}}}/(2\cdot w_{L})]_{{\mathord{\uparrow}}}={f_{\mathrm{min}}}\boxslash_{\mathord{\uparrow}}(2\cdot w_{L})$, by (29) of Proposition 5 we can conclude $y\succcurlyeq{f_{\mathrm{min}}}\boxslash_{\mathord{\uparrow}}(2\cdot w_{L})$.

Finally, for the case analysis of $\mathrm{im}_{\ell}$ marked as $a_{7}$ of Fig 18, the following cases must be considered:

$\bar{r}_{\ell}={\mathord{\downarrow}}:$

in this case we must have $y=+\infty$;

$\bar{r}_{\ell}={\mathord{\uparrow}}:$

it should be $y\cdot w_{L}>-{f_{\mathrm{max}}}$ and thus, since $w_{L}$ is positive, $y>{f_{\mathrm{max}}}/w_{L}$ and, by (26) of Proposition 5, $y\succcurlyeq\mathop{\mathrm{succ}}\nolimits({f_{\mathrm{max}}}\boxslash_{\mathord{\downarrow}}w_{L})$.

$\bar{r}_{\ell}={\mathrm{n}}:$

since $\mathop{\mathrm{odd}}\nolimits({f_{\mathrm{max}}})$, for $x_{L}=+\infty$ we need $y$ to be greater than or equal to \bigl{(}{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2\bigr{)}/w_{L}. If \llbracket\bigl{(}{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2\bigr{)}/w_{L}\rrbracket_{{\mathord{\uparrow}}}=[{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2\bigr{)}/w_{L}]_{{\mathord{\uparrow}}}, by (29) of Proposition 5, we can conclude y\succcurlyeq\llbracket\bigl{(}{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2\bigr{)}/w_{L}\rrbracket_{{\mathord{\uparrow}}}. On the other hand, if \llbracket\bigl{(}{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2\bigr{)}/w_{L}\rrbracket_{{\mathord{\uparrow}}}\neq[{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2\bigr{)}/w_{L}]_{{\mathord{\uparrow}}}, then we can only apply (25) of Proposition 5, obtaining y\succcurlyeq\llbracket\bigl{(}{f_{\mathrm{max}}}+\mathop{\nabla_{2}^{\mathrm{n}+}}\nolimits({f_{\mathrm{max}}})/2\bigr{)}/w_{L}\rrbracket_{{\mathord{\downarrow}}}.

Similar arguments can be used to prove the case analyses of $\mathrm{im}_{u}$ of Figure 19 marked as $a_{9}$, $a_{10}$, $a_{11}$ and $a_{12}$.

We now analyze the case analyses of $\mathrm{im}_{\ell}$ of Fig 18 marked as $a^{-}_{3}$ and $a^{+}_{3}$ and the ones of $\mathrm{im}_{u}$ of Fig 19 marked as $a^{-}_{8}$ and $a^{+}_{8}$, for which we can assume $x_{L},w_{L}\in\mathbb{F}\cap\mathbb{R}$ and $x_{U},w_{U}\in\mathbb{F}\cap\mathbb{R}$, and $\mathop{\mathrm{sgn}}\nolimits(w_{\ell})=\mathop{\mathrm{sgn}}\nolimits(w_{u})$. Exploiting $x\preccurlyeq y\boxdot z$ and $x\succcurlyeq y\boxdot z$, by Proposition 4 we have

[TABLE]

Since the case $z=0$ is handled separately by $\mathrm{im}_{\ell}$ of Figure 18 and by $\mathrm{im}_{u}$ of Figure 19, we can assume $z\neq 0$. Thanks to the splitting of $Z$ into a positive and a negative part, the sign of $z$ is determined. In the following, we will prove the case analyses marked as $a_{3}^{+}$ and $a_{8}^{+}$. Hence, assuming $z>0$, the previous case analysis gives us

[TABLE]

Note that the numerator and the denominator of the previous fractions are independent. Therefore, we can find the minimum of the fractions by minimizing the numerator and maximizing the denominator. Since we are analyzing the case in which $W=Z_{+}$, let $(x_{L},w_{L},x_{U},w_{U})$ as the result of function $\tau$ of Figure 5. Hence, by Proposition 3 and the monotonicity of ‘$\mathop{\mathrm{pred}}\nolimits$’ and ‘$\mathop{\mathrm{succ}}\nolimits$ we obtain

[TABLE]

We can now exploit Proposition 5 and obtain:

[TABLE]

Indeed, if $x_{L}\neq 0$, then Proposition 5 applies and we have $y\succcurlyeq x_{L}\boxslash_{\mathord{\uparrow}}w_{L}$. On the other hand, if $x_{L}=0$, since by hypothesis $z>0$ implies $w_{L}>0$, according to IEEE 754 (IEEE-754-2008, , Section 6.3), we have $(x_{L}\boxslash_{\mathord{\uparrow}}w_{L})=\mathop{\mathrm{sgn}}\nolimits(x_{L})\cdot 0$ and, indeed, for each non-NaN, nonzero and finite $w\in\mathbb{F}\cap[+0,+\infty]$, $\mathop{\mathrm{sgn}}\nolimits(x_{L})\cdot 0$ is the least value for $y$ that satisfies $\mathop{\mathrm{sgn}}\nolimits(x_{L})\cdot 0=y\boxdot_{\mathord{\downarrow}}w$.

Analogously, if $x_{L}\neq{f_{\mathrm{min}}}$, then Proposition 5 applies and we have \mathop{\mathrm{succ}}\nolimits\bigl{(}\mathop{\mathrm{pred}}\nolimits(x_{L})\boxslash_{\mathord{\downarrow}}w_{L}\bigr{)}. On the other hand, if $x_{L}={f_{\mathrm{min}}}$, \mathop{\mathrm{succ}}\nolimits\bigl{(}\mathop{\mathrm{pred}}\nolimits(x_{L})\boxslash_{\mathord{\downarrow}}w_{L}\bigr{)}={f_{\mathrm{min}}}, which is consistent with the fact that, for each non-NaN, nonzero and finite $w\in\mathbb{F}\cap[+0,+\infty]$, ${f_{\mathrm{min}}}$ is the lowest value of $y$ that satisfies ${f_{\mathrm{min}}}=y\boxdot_{\mathord{\uparrow}}w$.

A symmetric argument justifies (71).

As before, we will consider both the cases $\llbracket e^{+}_{\ell}\rrbracket_{{\mathord{\uparrow}}}=[e^{+}_{\ell}]_{{\mathord{\uparrow}}}$ and $\llbracket e^{+}_{u}\rrbracket_{{\mathord{\downarrow}}}=[e^{+}_{u}]_{{\mathord{\downarrow}}}$ as well as $\llbracket e^{+}_{\ell}\rrbracket_{{\mathord{\uparrow}}}>[e^{+}_{\ell}]_{{\mathord{\uparrow}}}$ and $\llbracket e^{+}_{u}\rrbracket_{{\mathord{\downarrow}}}<[e^{+}_{u}]_{{\mathord{\downarrow}}}$. Thus, when $\llbracket e^{+}_{u}\rrbracket_{{\mathord{\downarrow}}}<[e^{+}_{u}]_{{\mathord{\downarrow}}}$ by (69) and (27) of Proposition 5 we obtain $y\preccurlyeq\llbracket e^{+}_{u}\rrbracket_{{\mathord{\uparrow}}}$. Instead, when $\llbracket e^{+}_{\ell}\rrbracket_{{\mathord{\downarrow}}}>[e^{+}_{\ell}]_{{\mathord{\downarrow}}}$, by (69) and (25) of Proposition 5 we obtain $y\succcurlyeq\llbracket e^{+}_{\ell}\rrbracket_{{\mathord{\downarrow}}}$. In conclusion, for the case in which $\bar{r}_{\ell}={\mathrm{n}}$, since $e_{u}\neq 0$ and $e_{\ell}\neq 0$, by Proposition 5, we have

[TABLE]

An analogous reasoning with $z<0$ allows us to obtain the case analyses marked as $a_{3}^{-}$ and $a_{8}^{-}$.