¹¹institutetext: University of Perpignan, LAMPS laboratory, 52 Av. P. Alduy, Perpignan, France ²²institutetext: Numalis, Cap Omega, Rond-point Benjamin Franklin, Montpellier, France ²²email: {assale.adje,matthieu.martel,dorra.ben-khalifa}@univ-perp.fr

Fast and Efficient Bit-Level Precision Tuning

Assalé Adjé 11 Dorra Ben Khalifa 11 Matthieu Martel 1122

Abstract

In this article, we introduce a new technique for precision tuning. This problem consists of finding the least data types for numerical values such that the result of the computation satisfies some accuracy requirement. State of the art techniques for precision tuning use a try and fail approach. They change the data types of some variables of the program and evaluate the accuracy of the result. Depending on what is obtained, they change more or less data types and repeat the process. Our technique is radically different. Based on semantic equations, we generate an Integer Linear Problem (ILP) from the program source code. Basically, this is done by reasoning on the most significant bit and the number of significant bits of the values which are integer quantities. The integer solution to this problem, computed in polynomial time by a (real) linear programming solver, gives the optimal data types at the bit level. A finer set of semantic equations is also proposed which does not reduce directly to an ILP problem. So we use policy iteration to find the solution. Both techniques have been implemented and we show that our results encompass the results of state of the art tools.

Keywords:

Static analysis, computer arithmetic, integer linear problems, numerical accuracy, policy iteration.

1 Introduction

Let us consider a program $P$ computing some numerical result $R$ , typically but not necessarily in the IEEE754 floating-point arithmetic [1]. Precision tuning then consists of finding the smallest data types for all the variables and expressions of $P$ such that the result $R$ has some desired accuracy. These last years, much attention has been paid to this problem [6, 8, 10, 11, 12, 18]. Indeed, precision tuning makes it possible to save memory and, by way of consequence, it has a positive impact on the footprint of programs concerning energy consumption, bandwidth usage, computation time, etc.

A common point to all the techniques cited previously is that they follow a try and fail strategy. Roughly speaking, one chooses a subset $S$ of the variables of $P$ , assigns to them smaller data types (e.g. binary32 instead of binary64 [1]) and evaluates the accuracy of the tuned program $P^{\prime}$ . If the accuracy of the result returned by $P^{\prime}$ is satisfying then new variables are included in $S$ or even smaller data types are assigned to certain variables already in $S$ (e.g. binary16). Otherwise, if the accuracy of the result of $P^{\prime}$ is not satisfying, then some variables are removed from $S$ . This process is applied repeatedly, until a stable state is found. Existing techniques differ in their way to evaluate the accuracy of programs, done by dynamic analysis [10, 11, 12, 18] or by static analysis [6, 8] of $P$ and $P^{\prime}$ . They may also differ in the algorithm used to define $S$ , delta debugging being the most widespread method [18].

Anyway all these techniques suffer from the same combinatorial limitation: If $P$ has $n$ variables and if the method tries $k$ different data types then the search space contains $k^{n}$ configurations. They scale neither in the number $n$ of variables (even if heuristics such as delta debugging [18] or branch and bound [6] reduce the search space at the price of optimality) or in the number $k$ of data types which can be tried. In particular, bit level precision tuning, which consists of finding the minimal number of bits needed for each variable to reach the desired accuracy, independently of a limited number $k$ of data types, is not an option.

So the method introduced in this article for precision tuning of programs is radically different. Here, no try and fail method is employed. Instead, the accuracy of the arithmetic expressions assigned to variables is determined by semantic equations, in function of the accuracy of the operands. By reasoning on the number of significant bits of the variables of $P$ and knowing the weight of their most significant bit thanks to a range analysis performed before the tuning phase, we are able to reduce the problem to an Integer Linear Problem (ILP) which can be optimally solved in one shot by a classical linear programming solver (no iteration). Concerning the number $n$ of variables, the method scales up to the solver limitations and the solutions are naturally found at the bit level, making the parameter $k$ irrelevant. An important point is that the optimal solution to the continuous linear programming relaxation of our ILP is a vector of integers, as demonstrated in Section 4.2. By consequence, we may use a linear solver among real numbers whose complexity is polynomial [19] (contrarily to the linear solvers among integers whose complexity is NP-complete [16]). This makes our precision tuning method solvable in polynomial-time, contrarily to the existing exponential methods.

Next, we go one step further by introducing a second set of semantic equations. These new equations make it possible to tune even more the precision by being less pessimistic on the propagation of carries in arithmetic operations. However the problem do not reduce any longer to an ILP problem ( $\min$ and $\max$ operators are needed). Then we use policy iteration (PI) [7] to find efficiently the solution.

Both methods have been implemented inside a tool for precision tuning named XXX¹¹1In this article, XXX hides the actual name of our tool and missing references refer to our previous work for anonymity. [5]. Formerly, XXX was expressing the precision tuning problem as a set of first order logical propositions among relations between linear integer expressions. An SMT solver (Z3 in practice [15]) was used repeatedly to find the existence of a solution with a certain weight expressing the number of significant bits ( $\mathsf{nsb}$ ) of the variables [5]. In the present article, we compare experimentally our new methods to the SMT based method previously used by XXX [5] and to the Precimonious tool [10, 18]. These experiments on programs coming from mathematical libraries or other applicative domains such as IoT [3] show that the technique introduced in this article for precision tuning clearly encompasses the state of the art techniques.

The rest of this article is organized as follows. In the next section, we provide a motivating example. We then present in Section 3 some essential background on the functions needed for the constraint generation and also we detail the set of constraints for both ILP and PI methods. Section 4 presents the proofs of correctness. We end up in Section 5 by showing that our new techniques exhibits very good results in practice before concluding in Section 6.

2 Motivating Example

A motivating example to better explain our method is given by the code snippet of Figure 1. In this example, we aim at modeling the movement of a simple pendulum without damping. Let $l=0.5\ m$ be the length of this pendulum, $m=1\ kg$ its mass and $g=9.81\ m\cdot s^{-2}$ Newton’s gravitational constant. We denote by $\theta$ the tilt angle in radians as shown in Figure 1 (initially $\theta=\frac{\pi}{4}$ ). The Equation describing the movement of the pendulum is given in Equation (1).

Refer to caption — Figure 1: Top left: source program. Top right: pendulum movement for $\theta=\frac{\pi}{4}$ . Bottom left: program annotated with labels. Bottom right: program with inferred accuracies.

\small m\cdot l\cdot\frac{d^{2}\theta}{dt^{2}}=-m\cdot g\cdot\sin\theta

(1)

Equation (1) being a second order differential equation, we need to transform it into a system of two first order differential equations for resolution so we obtain $y_{1}=\theta$ and $y_{2}=\frac{d\theta}{dt}$ . By applying Euler’s method to these last equations, we obtain Equation (2) implemented in Figure 1.

\small\frac{dy_{1}}{dt}=y_{2}\quad\text{and}\quad\frac{dy_{2}}{dt}=-\frac{g}{l}\cdot\sin y_{1}

(2)

The key point of our technique is to generate a set of constraints for each statement of our imperative language introduced further in Section 3. For our example, we suppose that all variables, before XXX analysis, are in double precision (source program in the top left corner of Figure 1) and that a range determination is performed by dynamic analysis on the program variables (we plan to use a static analyzer in the future). XXX assigns to each node of the program’s syntactic tree a unique control point in order to determine easily the number of significant bits of the result as mentioned in the bottom corner of Figure 1. Some notations can be stressed about the structure of XXX source code. For instance, the annotation g ${}^{\ell_{1}}=9.81^{\ell_{0}}$ denotes that g has a unique control point $\ell_{1}$ . As well, we have the statement require_nsb(y2,20) which informs the tool that the user wants to get on variable y2 only 20 bits (we consider that a result has $n$ significants if the relative error between the exact and approximated results is less than $2^{-n}$ ). Finally, the minimal precision needed for the inputs and intermediary results satisfying the user assertion is observed on the bottom right corner of Figure 1. In this code, if we consider for instance Line $6$ , then y1new|20| means that the variable needs $20$ significant bits at this point. Similarly, y1 and y2 need $21$ bits each and the addition requires $20$ bits.

The originality of our method is that we reduce the precision tuning problem to an ILP. For example, taking again Line $6$ of the pendulum code, we generate six constraints as shown in Equation 3 (this is detailed further in Section 3).

\scriptsize C_{1}=\left\{\begin{array}[]{l}\mathsf{nsb}(\ell_{17})\geq\mathsf{nsb}(\ell_{23})+(-1)+\xi(\ell_{23})(\ell_{17},\ell_{22})-(-1),\\ \mathsf{nsb}(\ell_{22}\geq\mathsf{nsb}(\ell_{23}+0+\xi(\ell_{23})(\ell_{17},\ell_{22})-(1),\\ \mathsf{nsb}(\ell_{19})\geq\mathsf{nsb}(\ell_{22})+\xi(\ell_{22})(\ell_{19},\ell_{21})-1,\\ \mathsf{nsb}(\ell_{21})\geq\mathsf{nsb}(\ell_{22})+\xi(\ell_{22})(\ell_{19},\ell_{21})-1,\\ \mathsf{nsb}(\ell_{23})\geq\mathsf{nsb}(\ell_{24}),\ \xi(\ell_{23})(\ell_{17},\ell_{22})\geq 1,\ \ \xi(\ell_{22})(\ell_{19},\ell_{21})\geq 1\\ \end{array}\right\}

(3)

The first two constraints are for the addition. Basically, $\mathsf{nsb}(\ell_{23})$ stands for number of significant bits as described in Section 3.1. It represents the difference between the unit in the first place ( $\mathsf{ufp}$ , see Section 3.1) of the result of the sum and the $\mathsf{ufp}$ of its error denoted $\mathsf{ufp_{e}}(\ell_{23})$ . We have $\mathsf{ufp}(\ell_{17})=\mathsf{ufp}(\ell_{17})-\mathsf{nsb}(\ell_{17})$ . As mentioned previously, the $\mathsf{ufp}$ are computed by a prior range analysis. Then, at constraint generation time, they are constants. For our example, $\mathsf{ufp}(\ell_{17})=-1$ . This quantity occurs in the first constraints. The next two constraints are for the multiplication. The fourth constraint $\mathsf{nsb}(\ell_{23})\geq\mathsf{nsb}(\ell_{24})$ is for the assignment and the last two constraints are for the constant functions $\xi(\ell_{23})(\ell_{17},\ell_{22})$ and $\xi(\ell_{22})(\ell_{19},\ell_{21})$ respectively for the addition and multiplication (see next paragraphs and Section 3 for more details). Note that XXX generates such constraints for all the statements of the program.

For a user requirement of 20 bits on the variable y2 as shown in the original program of the top right corner of Figure 1 (all variables are in double precision initially), XXX succeeds in tuning the majority of variables of the pendulum program into simple precision with a total number of bits at bit level equivalent to 274 (originally the program used 689 bits). The new mixed precision formats obtained for line $6$ are y1new|20| = y1|21| +|20| y2|22| $\times$ |22| h|22|.

Let us now focus on the term $\xi(\ell_{23})(\ell_{17},\ell_{22})$ (for the addition). In our ILP, we always assume that $\xi$ is a constant function equal to $1$ . This corresponds to the carry bit which can be propagated up to the $\mathsf{ufp}$ and $\mathsf{ufp_{e}}$ and increments them by $1$ which is correct but pessimistic. In large codes, this function becomes very costly if we perform several computations at a time and therefore the errors would be considerable, especially that in many cases adding this carry bit is useless because the operands and their errors do not overlap. It is then crucial to use the most precise function $\xi$ . Unfortunately, when we model this optimization the problem is no more linear ( $\min$ and $\max$ operators arise) and we have to use the PI technique [7] to solve it (see Section 3.3).

In this case, by analyzing Line $6$ of our program, we have to add the following new constraints (along with the former ones) as mentioned in Equation (4). In fact, policy iteration makes it possible to break the $\min$ in the $\xi(\ell_{23})(\ell_{17},\ell_{22})$ function by choosing the $\max$ between $\mathsf{ufp}(\ell_{22})-\mathsf{ufp}(\ell_{17})-\mathsf{nsb}(\ell_{17})-\mathsf{nsb}(\ell_{22})-\mathsf{nsb_{e}}(\ell_{17})$ and $0$ , the $\max$ between $\mathsf{ufp}(\ell_{17})-\mathsf{ufp}(\ell_{22})+\mathsf{nsb}(\ell_{22})-\mathsf{nsb}(\ell_{17})-\mathsf{nsb_{e}}(\ell_{22})$ and $0$ and the constant $1$ . Next, it becomes possible to solve the corresponding ILP. If no fixed point is reached, XXX iterates until a solution is found as shown in Section 3.3. By applying this optimization, the new formats of line $6$ are given as y1new|20| = y1|21| +|20| y2|21| $\times$ |22| h|21|. By comparing with the formats obtained above, a gain of precision of $1$ bit is observed on variables y2 and h (total of 272 bits at bit level for the optimized program). The program on the bottom right corner of Figure 1 illustrates the new optimized formats obtained by using the policy iteration technique.

\scriptsize C_{2}=\left\{\begin{array}[]{l}\mathsf{nsb_{e}}(\ell_{23})\geq\mathsf{nsb_{e}}(\ell_{17}),\\ \mathsf{nsb_{e}}(\ell_{23})\geq\mathsf{nsb_{e}}(\ell_{22}),\\ \mathsf{nsb}(\ell_{23})\geq-1-0+\mathsf{nsb}(\ell_{22})-\mathsf{nsb}(\ell_{17})+\mathsf{nsb_{e}}(\ell_{22})+\xi(\ell_{23},\ell_{17},\ell_{22}),\\ \mathsf{nsb_{e}}(\ell_{23})\geq 0-(-1)+\mathsf{nsb}(\ell_{17})-\mathsf{nsb}(\ell_{22})+\mathsf{nsb_{e}}(\ell_{17})+\xi(\ell_{23},\ell_{17},\ell_{22}),\\ \mathsf{nsb_{e}}(\ell_{23})\geq\mathsf{nsb_{e}}(\ell_{24}),\\ \mathsf{nsb_{e}}(\ell_{22})\geq\mathsf{nsb}(\ell_{19})+\mathsf{nsb_{e}}(\ell_{19})+\mathsf{nsb_{e}}(\ell_{21})-2,\\ \mathsf{nsb_{e}}(\ell_{22})\geq\mathsf{nsb}(\ell_{21})+\mathsf{nsb_{e}}(\ell_{21})+\mathsf{nsb_{e}}(\ell_{19})-2,\\ \xi(\ell_{23})(\ell_{17},\ell_{22})=\min\left(\begin{array}[]{l}\max\big{(}0-6+\mathsf{nsb}(\ell_{17})-\mathsf{nsb}(\ell_{22})-\mathsf{nsb_{e}}(\ell_{17}),0\big{)},\\ \max\big{(}6-0+\mathsf{nsb}(\ell_{22})-\mathsf{nsb}(\ell_{17})-\mathsf{nsb_{e}}(\ell_{22}),0\big{)},1\end{array}\right)\\ \end{array}\right\}

(4)

\varepsilon(c^{\ell}\texttt{\#p})\leq 2^{\mathsf{ufp}(c)-\min\big{(}\texttt{p},\mathsf{prec}(\ell)\big{)}}

(5)

\varepsilon(c_{1}^{\ell_{1}}\texttt{\#$p_{1}$}+^{\ell}c_{2}^{\ell_{2}}\texttt{\#$p_{2}$})\leq\varepsilon({c_{1}^{\ell_{1}}\texttt{\#$p_{1}$}})+\varepsilon({c_{2}^{\ell_{2}}\texttt{\#$p_{2}$}})+{2^{\mathsf{ufp}(c_{1}+c_{2})-\mathsf{prec}(\ell)}}

(6)

\varepsilon(c_{1}^{\ell_{1}}\texttt{\#$p_{1}$}-^{\ell}c_{2}^{\ell_{2}}\texttt{\#$p_{2}$})\leq\varepsilon({c_{1}^{\ell_{1}}\texttt{\#$p_{1}$}})-\varepsilon({c_{2}^{\ell_{2}}\texttt{\#$p_{2}$}})+{2^{\mathsf{ufp}(c_{1}-c_{2})-\mathsf{prec}(\ell)}}

(7)

\begin{array}[]{c}\varepsilon(c_{1}^{\ell_{1}}\texttt{\#$p_{1}$}\times^{\ell}c_{2}^{\ell_{2}}\texttt{\#$p_{2}$})\leq\\ c_{1}\cdot\varepsilon({c_{2}^{\ell_{2}}\texttt{\#$p_{2}$}})+c_{2}\cdot\varepsilon({c_{1}^{\ell_{1}}\texttt{\#$p_{1}$}})+\varepsilon({c_{1}^{\ell_{1}}\texttt{\#$p_{1}$}})\cdot\varepsilon({c_{2}^{\ell_{2}}\texttt{\#$p_{2}$}})+{2^{\mathsf{ufp}(c_{1}\times c_{2})-\mathsf{prec}(\ell)}}\end{array}

(8)

\varepsilon(c_{1}^{\ell_{1}}\texttt{\#$p_{1}$}\div^{\ell}c_{2}^{\ell_{2}}\texttt{\#$p_{2}$})\leq\varepsilon(c_{1}^{\ell_{1}}\texttt{\#$p_{1}$}\times^{\ell}c_{2}^{\prime\ell_{2}}\texttt{\#$p_{2}$})\quad\text{with}\quad c^{\prime}_{2}=\frac{1}{c_{2}}

(9)

\varepsilon\left(\phi(c^{\ell_{1}}\texttt{\#p})^{\ell}\right)\leq 2^{\mathsf{ufp}(\phi(c))-\texttt{p}+\varphi}+2^{\mathsf{ufp}(\phi(c))-\mathsf{prec}(\ell)}\ \text{with}\quad\phi\in\{\sin,\cos,\tan,\log,\ldots\}

(10)

\varepsilon\left(\sqrt{(c^{\ell_{1}}\texttt{\#p})}^{\ell}\right)\leq 2^{\mathsf{ufp}(\sqrt{c})-\texttt{p}}+2^{\mathsf{ufp}(\sqrt{c})-\mathsf{prec}(\ell)}

(11)

Figure 2: Numerical error on arithmetic expressions.

3 Constraints Generation for Bit-Level Precision Tuning

In this section, we start by providing essential definitions for understanding the rest of the article. Also, we define a simple imperative language (see Figure 3) from which we generate semantic equations in order to determine the least precision needed for the program numerical values. Then, we will focus on the difference between the two sets of constraints obtained when using the simple ILP and the more complex PI formulations which optimizes the carry bit that can propagate throughout computations. The operational semantics of the language as well as the theorem proving that the solution to the system of constraints gives the desired $\mathsf{nsb}$ when running programs are detailed further in Section 4.

3.1 Elements of Computer Arithmetic

Our technique is independent of a particular computer arithmetic. In fact, we manipulate numbers for which we know their unit in the first place ( $\mathsf{ufp}$ ) and the number of significant digits ( $\mathsf{nsb}$ ) defined as follows.

Unit in the First Place

The unit in the first place of a real number $x$ (possibly encoded up to some rounding mode by a floating-point or a fixpoint number) is given in Equation (12). This function, which is independent of the representation of $x$ , will be used further in this section to describe the error propagation across the computations.

\small\mathsf{ufp}(x)=\min\{i\in\mathbb{Z}:2^{i+1}>x\}=\lfloor\log_{2}(x)\rfloor\enspace.

(12)

Number of Significant Bits

Intuitively, $\mathsf{nsb}(x)$ is the number of significant bits of $x$ . Formally, following Parker [17], if $\mathsf{nsb}(x)=k$ , for $x\not=0$ then the error $\varepsilon(x)$ on $x$ is less than $2^{\mathsf{ufp}(x)-k}$ . If $x=0$ then $\mathsf{nsb}(x)=0$ . For example, if the exact binary value $1.0101$ is approximated by either $x=1.010$ or $x=1.011$ then $\mathsf{nsb}(x)=3$ .

In the following, we also use $\mathsf{ufp_{e}}(x)$ and $\mathsf{nsb_{e}}(x)$ to denote the $\mathsf{ufp}$ and $\mathsf{nsb}$ of the error on $x$ , i.e. $\mathsf{ufp_{e}}(x)=\mathsf{ufp}(\varepsilon(x))$ and $\mathsf{nsb_{e}}(x)=\mathsf{nsb}(\varepsilon(x))$ .

In this article, we consider a finite precision arithmetic, independently of any particular representation (IEEE754 [1], POSIT [7], $\ldots$ ). Nevertheless, the representation being finite, roundoff errors may arise when representing values or performing elementary operation. These errors are defined in Figure 2. First of all, by definition, using the function $\mathsf{ufp}$ of Equation (12), for a number $x$ with p number of significant bits, the roundoff error $\epsilon(x)$ is bounded as shown in Equation (13).

\small\varepsilon(x)\leq 2^{\mathsf{ufp}(x)-\texttt{p}+1}

(13)

Let $\mathsf{prec}(\ell)$ be the precision of the operation at control point $\ell$ . For example, the precision is $53$ bits for the IEEE754 binary64 format. In fact, $\mathsf{prec}(\ell)$ is used to compute the truncation error of the operations. Equations (5) to (11) of Figure 2 define the numerical errors of the arithmetic expressions of our language (presented in Section 3). For constants occurring in the code, the initial precision must be given by the user and we write $c^{\ell}\texttt{\#p}$ a constant $c$ with p significant bits at control point $\ell$ . Then, following Equation (13), $\varepsilon(c^{\ell}\texttt{\#p})$ is defined in Equation (5) of Figure 2: the $\mathsf{nsb}$ of the constant is $\min(\texttt{p},\mathsf{prec}(\ell))$ and consequently the error is bounded by $2^{\mathsf{ufp}(c)-\min\big{(}\texttt{{p}},\mathsf{prec}(\ell)\big{)}}$ . In equations (6) to (9), we propagate the errors on the operands and we add the roundoff error due to the elementary operation itself. For an elementary function $\phi$ , we assume that $\varphi$ bits are lost as shown in Equation (10) (more details are given in Section 3). The last equation is for the square root function. This function being computable exactly, no more error than for the elementary operation is introduced.

3.2 Integer Linear Problem Formulation

First, we define in Figure 3 the simple imperative language in which our input programs are written.

$\ell\in Lab$ $x\in Id$ $\odot$ $\in$ {+, -, $\times$ , $\div$ } $math$ $\in$ { $\sin$ , $\cos$ , $\tan$ , $\arcsin$ , $\log$ , $\ldots$ }

Expr $\ni$ e : e ::= c# $p^{\ell}$ $|$ $x^{\ell}$ $|$ $e_{1}^{\ell_{1}}\odot^{\ell}e_{2}^{\ell_{2}}$ $|$ $math(e^{\ell_{1}})^{\ell}$ $|$ $sqrt(e^{\ell_{1}})^{\ell}$

Cmd $\ni$ c : $c::=c_{1}^{\ell_{1}};c_{2}^{\ell_{2}}$ $|$ $x=^{\ell}e^{\ell_{1}}$ $|$ $\textbf{{while}}^{\ell}\>b^{\ell_{0}}\>\textbf{{do}}\>\>c_{1}^{\ell_{1}}$ $|$ $\textbf{{if}}^{\ell}\>b^{\ell_{0}}\>\textbf{{then}}\>\>c_{1}^{\ell_{1}}\>\textit{{else}}\>\>c$ $|$ $\textbf{{require\_nsb}}(x,n)^{\ell}$

Figure 3: Simple imperative language of constraints.

\mathcal{E}[c\texttt{\#p}^{\ell}]\varrho=\emptyset\quad(\textsc{Const})\hskip 28.45274pt\mathcal{E}[x^{\ell}]\varrho=\big{\{}\mathsf{nsb}(\varrho(x))\geq\mathsf{nsb}(\ell)\big{\}}\quad(\textsc{Id})

\begin{array}[]{c}\mathcal{E}[e_{1}^{\ell_{1}}+^{\ell}e_{2}^{\ell_{2}}]\varrho=\mathcal{E}[e_{1}^{\ell_{1}}]\varrho\ \cup\ \mathcal{E}[e_{2}^{\ell_{2}}]\varrho\\ \cup\\ \left\{\mathsf{nsb}(\ell_{1})\geq\mathsf{nsb}(\ell)+\mathsf{ufp}(\ell_{1})-\mathsf{ufp}(\ell)+\xi(\ell)(\ell_{1},\ell_{2}),\right.\\ \left.\hskip 2.84544pt\mathsf{nsb}(\ell_{2})\geq\mathsf{nsb}(\ell)+\mathsf{ufp}(\ell_{2})-\mathsf{ufp}(\ell)+\xi(\ell)(\ell_{1},\ell_{2})\right\}\end{array}\quad(\textsc{Add})

\begin{array}[]{c}\mathcal{E}[e_{1}^{\ell_{1}}-^{\ell}e_{2}^{\ell_{2}}]\varrho=\mathcal{E}[e_{1}^{\ell_{1}}]\varrho\ \cup\ \mathcal{E}[e_{2}^{\ell_{2}}]\varrho\\ \cup\\ \left\{\mathsf{nsb}(\ell_{1})\geq\mathsf{nsb}(\ell)+\mathsf{ufp}(\ell_{1})-\mathsf{ufp}(\ell)+\xi(\ell)(\ell_{1},\ell_{2}),\right.\\ \left.\hskip 2.84544pt\mathsf{nsb}(\ell_{2})\geq\mathsf{nsb}(\ell)+\mathsf{ufp}(\ell_{2})-\mathsf{ufp}(\ell)+\xi(\ell)(\ell_{1},\ell_{2})\right\}\end{array}\quad(\textsc{Sub})

\begin{array}[]{c}\mathcal{E}[e_{1}^{\ell_{1}}\times^{\ell}e_{2}^{\ell_{2}}]\varrho=\mathcal{E}[e_{1}^{\ell_{1}}]\varrho\ \cup\ \mathcal{E}[e_{2}^{\ell_{2}}]\varrho\\ \cup\\ \left\{\mathsf{nsb}(\ell_{1})\geq\mathsf{nsb}(\ell)+\xi(\ell)(\ell_{1},\ell_{2})-1,\right.\left.\hskip 2.84544pt\mathsf{nsb}(\ell_{2})\geq\mathsf{nsb}(\ell)+\xi(\ell)(\ell_{1},\ell_{2})-1\right\}\end{array}\quad(\textsc{Mult})

\begin{array}[]{c}\mathcal{E}[e_{1}^{\ell_{1}}\div^{\ell}e_{2}^{\ell_{2}}]\varrho=\mathcal{E}[e_{1}^{\ell_{1}}]\varrho\ \cup\ \mathcal{E}[e_{2}^{\ell_{2}}]\varrho\\ \cup\\ \left\{\mathsf{nsb}(\ell_{1})\geq\mathsf{nsb}(\ell)+\xi(\ell)(\ell_{1},\ell_{2})-1,\right.\left.\hskip 2.84544pt\mathsf{nsb}(\ell_{2})\geq\mathsf{nsb}(\ell)+\xi(\ell)(\ell_{1},\ell_{2})-1\right\}\end{array}\quad(\textsc{Div})

\mathcal{E}\left[\sqrt{e^{\ell_{1}}}^{\ell}\right]\varrho=\mathcal{E}[e_{1}^{\ell_{1}}]\varrho\ \cup\ \big{\{}\mathsf{nsb}(\ell_{1})\geq\mathsf{nsb}(\ell)\big{\}}\quad(\textsc{Sqrt})

\mathcal{E}\left[\phi\big{(}{e^{\ell_{1}}}\big{)}^{\ell}\right]\varrho=\mathcal{E}[e_{1}^{\ell_{1}}]\varrho\ \cup\ \big{\{}\mathsf{nsb}(\ell_{1})\geq\mathsf{nsb}(\ell)+\varphi\big{\}}\ \text{with}\ \phi\in\{\sin,\cos,\tan,\log,\ldots\}\quad(\textsc{Math})

\begin{array}[]{c}\mathcal{C}\left[x\texttt{:=}^{\ell}e^{\ell_{1}}\right]\varrho=\big{(}C,\varrho\left[x\mapsto\ell\right]\big{)}\ \text{where}\ C=\mathcal{E}[e_{1}^{\ell_{1}}]\varrho\cup\left\{\mathsf{nsb}(\ell_{1})\geq\mathsf{nsb}(\ell)\right\}\end{array}\quad(\textsc{Assign})

\begin{array}[]{c}\mathcal{C}\left[c_{1}^{\ell_{1}}\texttt{;}c_{2}^{\ell_{2}}\right]\varrho=\big{(}C_{1}\cup C_{2},\varrho_{2}\big{)}\par\\ \text{where}\ \big{(}C_{1},\varrho_{1}\big{)}=\mathcal{C}\left[c_{1}^{\ell_{1}}\right]\varrho\ \text{and}\ \big{(}C_{2},\varrho_{2}\big{)}=\mathcal{C}\left[c_{2}^{\ell_{2}}\right]\varrho_{1}\end{array}\quad(\textsc{Seq})

\begin{array}[]{c}\mathcal{C}[\texttt{if}^{\ell}\ e^{\ell_{0}}\ \texttt{then}\ c^{\ell_{1}}\ \texttt{else}\ c^{\ell_{2}}]\ \varrho=(C_{1}\cup C_{2}\cup C_{3},\varrho^{\prime})\\ \text{where}\ \left|\begin{array}[]{l}\forall x\in\text{Id},\ \varrho^{\prime}(x)=\ell,\ (C_{1},\varrho_{1})=\mathcal{C}[c_{1}^{\ell_{1}}]\ \varrho,\ (C_{2},\varrho_{2})=\mathcal{C}[c_{2}^{\ell_{2}}]\ \varrho,\\ C_{3}=\underset{x\in\text{Id}}{\bigcup}\left\{\mathsf{nsb}(\varrho_{1}(x))\geq\mathsf{nsb}(\ell),\ \mathsf{nsb}(\varrho_{2}(x))\geq\mathsf{nsb}(\ell)\right\}\end{array}\right.\end{array}\quad(\textsc{Cond})

\begin{array}[]{c}\mathcal{C}[\texttt{while}^{\ell}\ e^{\ell_{0}}\ \texttt{do}\ c^{\ell_{1}}]\ \varrho=(C_{1}\cup C_{2},\varrho^{\prime})\\ \text{where}\ \left|\begin{array}[]{l}\forall x\in\text{Id},\ \varrho^{\prime}(x)=\ell,\ (C_{1},\varrho_{1})=\mathcal{C}[c_{1}^{\ell_{1}}]\ \varrho^{\prime}\\ C_{2}=\underset{x\in\text{Id}}{\bigcup}\left\{\mathsf{nsb}(\varrho(x))\geq\mathsf{nsb}(\ell),\ \mathsf{nsb}(\varrho_{1}(x))\geq\mathsf{nsb}(\ell)\right\}\end{array}\right.\end{array}\quad(\textsc{While})

\mathcal{C}[\texttt{require\_nsb}(x,\texttt{p})^{\ell}]\varrho=\big{\{}\mathsf{nsb}(\varrho(x))\geq\texttt{p}\big{\}}\quad(\textsc{Req})

$\xi(\ell)(\ell_{1},\ell_{2})=1$

Figure 4: ILP constraints with pessimistic carry bit propagation

\xi=1

We denote by $Id$ the set of identifiers and by $Lab$ the set of control points of the program as a means to assign to each element $e\in Expr$ and $c\in Cmd$ of our language a unique control point $\ell$ $\in$ $Lab$ . First, in $c\texttt{\#p}$ , p indicates the number of significant bits of the constant $c$ in the source code. The parameter p is computed by our tool XXX, when solving the constraints. Next, the statement require_nsb(x,n)^ℓ indicates the number of significant bits $n$ that a variable $x$ must have at a control point $\ell$ . The rest of the grammar is standard.

As we have mentioned, we are able to reduce the problem of determining the lowest precision on variables and intermediary values in programs to an Integer Linear Problem (ILP) by reasoning on their unit in the first place ( $\mathsf{ufp}$ ) and the number of significant bits. In addition, we assign to each control point $\ell$ an integer variable $\mathsf{nsb}(\ell)$ corresponding to the $\mathsf{nsb}$ of the arithmetic expressions. $\mathsf{nsb}(\ell)$ is determined by solving the ILP generated by the rules of Figure 4.

Let us now focus on the rules of Figure 4 where $\varrho\ :\ \text{Id}\rightarrow\text{Id}\times\text{Lab}$ is an environment which relates each identifier $x$ to its last assignment $x^{\ell}$ : Assuming that $x\ \texttt{:=}^{\ell}e^{\ell_{1}}$ is the last assignment of $x$ , the environment $\varrho$ maps $x$ to $x^{\ell}$ . Then, $\mathcal{E}[e]\ \varrho$ generates the set of constraints for an expression $e\in Expr$ in the environment $\varrho$ . We now formally define these constraints for each element of our language. No constraint is generated for a constant $c\texttt{\#p}$ as mentioned in Rule (Const) of Figure 4. For Rule (Id) of a variable $x^{\ell}$ , we require that the $\mathsf{nsb}$ at control point $\ell$ is less than its $\mathsf{nsb}$ in the last assignment of $x$ given in $\varrho(x)$ . For a binary operator $\odot$ $\in$ {+, -, $\times$ , $\div$ }, we first generate the set of constraints $\mathcal{E}[e_{1}^{\ell_{1}}]\varrho$ and $\mathcal{E}[e_{2}^{\ell_{2}}]\varrho$ for the operands at control points $\ell_{1}$ and $\ell_{2}$ . Considering Rule (ADD), the result of the addition of two numbers is stored in control point $\ell$ . Recall that a range determination is performed before the accuracy analysis, $\mathsf{ufp}(\ell)$ , $\mathsf{ufp}(\ell_{1})$ and $\mathsf{ufp}(\ell_{2})$ are known at constraint generation time.

Now, before going further in the explanation of the constraints generation for binary operations, we introduce the function $\xi$ which computes the carry bit that can occur throughout an addition (similar reasoning will be done for the other elementary operations). In the present ILP of Figure 4, we over-approximate the function $\xi$ by $\xi(\ell)(\ell_{1},\ell_{2})=1$ for all $\ell,\ \ell_{1}$ and $\ell_{2}$ , thus assuming the worst case, i.e. a carry bit is added at each operation. We will optimize $\xi$ in Section 3.3 but the problem will not remain an ILP any longer. To wrap up, for the addition (Rule (Add)), the $\mathsf{nsb}(\ell)$ of the exact result is the number of bits between $\mathsf{ufp}(\ell_{1}+\ell_{2})$ and the $\mathsf{ufp}$ of the error $e$ which is:

\small e=\max\big{(}\mathsf{ufp}(\ell_{1})-\mathsf{nsb}(\ell_{1}),\mathsf{ufp}(\ell_{2})-\mathsf{nsb}(\ell_{2})\big{)}-\xi(\ell)(\ell_{1},\ell_{2})

(14)

Hence, the error on the addition in precision $\mathsf{prec}(\ell)$ is

\small e^{\prime}=\max\big{(}\mathsf{ufp}(\ell_{1})-\mathsf{nsb}(\ell_{1}),\mathsf{ufp}(\ell_{2})-\mathsf{nsb}(\ell_{2}),\mathsf{prec}(\ell)\big{)}-\xi(\ell)(\ell_{1},\ell_{2})

(15)

Equation (14) is obtained from Equation (6):

$\displaystyle\small\varepsilon(c_{1}^{\ell_{1}}\texttt{\#$p_{1}$}+^{\ell}c_{2}^{\ell_{2}}\texttt{\#$p_{2}$})$	$\displaystyle\leq$	$\displaystyle\varepsilon({c_{1}^{\ell_{1}}\texttt{\#$p_{1}$}})+\varepsilon({c_{2}^{\ell_{2}}\texttt{\#$p_{2}$}})+{2^{\mathsf{ufp}(c_{1}+c_{2})-\mathsf{prec}(\ell)}}$
	$\displaystyle=$	$\displaystyle\max\big{(}\varepsilon({c_{1}^{\ell_{1}}\texttt{\#$p_{1}$}}),\varepsilon({c_{2}^{\ell_{2}}\texttt{\#$p_{2}$}})\big{)}+{2^{\mathsf{ufp}(c_{1}+c_{2})-\mathsf{prec}(\ell)}}-\xi(\ell)(\ell_{1},\ell_{2})$
	$\displaystyle=$	$\displaystyle\max\big{(}\mathsf{ufp}(\ell_{1})-\mathsf{nsb}(\ell_{1}),\mathsf{ufp}(\ell_{1})-\mathsf{nsb}(\ell_{1})\big{)}$
		$\displaystyle+{2^{\mathsf{ufp}(c_{1}+c_{2})-\mathsf{prec}(\ell)}}-\xi(\ell)(\ell_{1},\ell_{2})$

Since $\mathsf{nsb}(\ell)\leq\mathsf{prec}(\ell)$ , we may get rid of the last term of $e^{\prime}$ in Equation (15) and the two equations generated for Rule (ADD) are derived from Equation (16).

\small\mathsf{nsb}(\ell)\leq\mathsf{ufp}(\ell)-\max\big{(}\mathsf{ufp}(\ell_{1})-\mathsf{nsb}(\ell_{1}),\mathsf{ufp}(\ell_{2})-\mathsf{nsb}(\ell_{2})\big{)}-\xi(\ell)(\ell_{1},\ell_{2})

(16)

For example, let us consider the piece of code hereafter.

⬇

1 x

{}^{\ell_{2}}

= 5.0

{}^{\ell_{1}}

; y

{}^{\ell_{4}}

= 3.0

{}^{\ell_{3}}

;

2 z

{}^{\ell_{8}}

= x

{}^{\ell_{6}}

{}^{\ell_{5}}

{}^{\ell_{7}}

;

3 require_nsb(z,15)

{}^{\ell_{9}}

;

Wrapping up our constraints, we have $\ 15\leq\mathsf{nsb}(\ell_{9})\leq\mathsf{nsb}(\ell_{8})\leq\mathsf{nsb}(\ell_{5}),$ and $\mathsf{nsb}(\ell_{6})\geq\mathsf{nsb}(\ell_{8})+\mathsf{ufp}(\ell_{6})-\mathsf{ufp}(\ell_{7})+\xi(\ell_{8})(\ell_{6},\ell_{7})$ . Since $\mathsf{ufp}(\ell_{6})=2$ and $\mathsf{ufp}(\ell_{7})=1$ we have $\mathsf{nsb}(\ell_{6})\geq\mathsf{nsb}(\ell_{8})+1+\xi(\ell_{8})(\ell_{6},\ell_{7})$ and, consequently, $\mathsf{nsb}(\ell_{6})\geq 15+1+1=17$ . Rule (Sub) for the subtraction is obtained similarly to the addition case. For Rule (Mult) of multiplication (and in the same manner Rule(Div)), the reasoning mimics the one of the addition. Let $x$ and $y$ be two floating point numbers and $z$ the result of their product, $z=x^{\ell_{1}}\times^{\ell}y^{\ell_{2}}$ . We denote by $\varepsilon(x)$ , $\varepsilon(y)$ and $\varepsilon(z)$ the errors on $x$ , $y$ and $z$ , respectively. The error $\varepsilon(z)$ of this multiplication is $\varepsilon(z)=x\cdot\varepsilon(y)+y\cdot\varepsilon(x)+\varepsilon(x)\cdot\varepsilon(y)$ . These numbers are bounded as shown in Equation (17).

\small\begin{array}[]{rcl}2^{\mathsf{ufp}(x)}\leq x\leq 2^{\mathsf{ufp}(x)+1}\quad\text{and}\quad 2^{\mathsf{ufp}(x)-\mathsf{nsb}(x)}&\leq&\varepsilon(x)\leq 2^{\mathsf{ufp}(x)-\mathsf{nsb}(x)+1}\\ 2^{\mathsf{ufp}(y)}\leq y\leq 2^{\mathsf{ufp}(y)+1}\quad\text{and}\quad 2^{\mathsf{ufp}(y)-\mathsf{nsb}(y)}&\leq&\varepsilon(y)\leq 2^{\mathsf{ufp}(y)-\mathsf{nsb}(y)+1}\\ 2^{\mathsf{ufp}(x)+\mathsf{ufp}(y)-\mathsf{nsb}(y)}+2^{\mathsf{ufp}(y)+\mathsf{ufp}(x)-\mathsf{nsb}(x)}&\leq&\varepsilon(z)\leq 2^{\mathsf{ufp}(z)-\mathsf{nsb}(z)+1}\\ +2^{\mathsf{ufp}(x)+\mathsf{ufp}(y)-\mathsf{nsb}(x)-\mathsf{nsb}(y)}&&\end{array}

(17)

By getting rid of the last term $2^{\mathsf{ufp}(x)+\mathsf{ufp}(y)-\mathsf{nsb}(x)-\mathsf{nsb}(y)}$ of the error $\epsilon(z)$ which is strictly less than the former two ones, assuming that $\mathsf{ufp}(x+y)=\mathsf{ufp}(z)$ and, finally, by reasoning on the exponents, we obtain the equations of Rule (Mult).

\small\mathsf{nsb}(x)\geq\mathsf{nsb}(z)+\xi(\ell)(\ell_{1},\ell_{2})-1\quad\text{and}\quad\mathsf{nsb}(y)\geq\mathsf{nsb}(y)+\xi(\ell)(\ell_{1},\ell_{2})-1

Although the square root (Sqrt) is included, e.g, in the IEEE754 Standard, it is not the case for the other elementary functions such as the natural logarithm, the exponential functions and the hyperbolic and trigonometric functions gathered in Rule (Math). Also, each implementation of these functions has its own $\mathsf{nsb}$ which we have to know to model the propagation of errors in our analyses. To cope with this limitation, we consider that each elementary function introduces a loss of precision of $\varphi$ bits, where $\varphi\in\mathbb{N}$ is a parameter of the analysis and consequently of our tool, XXX.

The rules of commands are rather classical, we use control points to distinguish many assignments of the same variable and also to implement joins in conditions and loops. Given a command $c$ and an environment $\varrho$ , $\mathcal{C}[c]\ \varrho$ returns a pair $(C,\varrho^{\prime})$ made of a set $C$ of constraints and of a new environment $\varrho^{\prime}$ . The function $\mathcal{C}$ is defined by induction on the structure of commands in figures 4 and 5. For conditionals, we generate the constraints for the then and else branches plus additional constraints to join the results of both branches. For loops, we relate the number of significants bits at the end of the body to the $\mathsf{nsb}$ of the same variables and the beginning of the loop.

\mathcal{E}^{\prime}[c\texttt{\#p}^{\ell}]\varrho=\big{\{}\mathsf{nsb_{e}}(\ell)=0\big{\}}\quad(\textsc{Const}^{\prime})\hskip 28.45274pt\mathcal{E}^{\prime}[x^{\ell}]\varrho=\big{\{}\mathsf{nsb_{e}}(\varrho(x))\geq\mathsf{nsb_{e}}(\ell)\big{\}}\quad(\textsc{Id}^{\prime})

\begin{array}[]{c}\mathcal{E}^{\prime}[e_{1}^{\ell_{1}}+^{\ell}e_{2}^{\ell_{2}}]\varrho=\mathcal{E}^{\prime}[e_{1}^{\ell_{1}}]\varrho\ \cup\ \mathcal{E}^{\prime}[e_{2}^{\ell_{2}}]\varrho\quad(\textsc{Add}^{\prime})\\ \cup\\ \left\{\begin{array}[]{c}\mathsf{nsb_{e}}(\ell)\geq\mathsf{nsb_{e}}(\ell_{1}),\ \mathsf{nsb_{e}}(\ell)\geq\mathsf{nsb_{e}}(\ell_{2}),\\ \mathsf{nsb_{e}}(\ell)\geq\mathsf{ufp}(\ell_{1})-\mathsf{ufp}(\ell_{2})+\mathsf{nsb}(\ell_{2})-\mathsf{nsb}(\ell_{1})+\mathsf{nsb_{e}}(\ell_{2})+\xi(\ell)(\ell_{1},\ell_{2}),\\ \mathsf{nsb_{e}}(\ell)\geq\mathsf{ufp}(\ell_{2})-\mathsf{ufp}(\ell_{1})+\mathsf{nsb}(\ell_{1})-\mathsf{nsb}(\ell_{2})+\mathsf{nsb_{e}}(\ell_{1})+\xi(\ell)(\ell_{1},\ell_{2})\end{array}\right\}\end{array}

\begin{array}[]{c}\mathcal{E}^{\prime}[e_{1}^{\ell_{1}}-^{\ell}e_{2}^{\ell_{2}}]\varrho=\mathcal{E}^{\prime}[e_{1}^{\ell_{1}}]\varrho\ \cup\ \mathcal{E}^{\prime}[e_{2}^{\ell_{2}}]\varrho\quad(\textsc{Sub}^{\prime})\\ \cup\\ \left\{\begin{array}[]{c}\mathsf{nsb_{e}}(\ell)\geq\mathsf{nsb_{e}}(\ell_{1}),\ \mathsf{nsb_{e}}(\ell)\geq\mathsf{nsb_{e}}(\ell_{2}),\\ \mathsf{nsb_{e}}(\ell)\geq\mathsf{ufp}(\ell_{1})-\mathsf{ufp}(\ell_{2})+\mathsf{nsb}(\ell_{2})-\mathsf{nsb}(\ell_{1})+\mathsf{nsb_{e}}(\ell_{2})+\xi(\ell)(\ell_{1},\ell_{2}),\\ \mathsf{nsb_{e}}(\ell)\geq\mathsf{ufp}(\ell_{2})-\mathsf{ufp}(\ell_{1})+\mathsf{nsb}(\ell_{1})-\mathsf{nsb}(\ell_{2})+\mathsf{nsb_{e}}(\ell_{1})+\xi(\ell)(\ell_{1},\ell_{2})\end{array}\right\}\end{array}

\begin{array}[]{c}\mathcal{E}^{\prime}[e_{1}^{\ell_{1}}\times^{\ell}e_{2}^{\ell_{2}}]\varrho=\mathcal{E}^{\prime}[e_{1}^{\ell_{1}}]\varrho\ \cup\ \mathcal{E}^{\prime}[e_{2}^{\ell_{2}}]\varrho\quad(\textsc{Mult}^{\prime})\\ \cup\\ \left\{\begin{array}[]{c}\mathsf{nsb_{e}}(\ell)\geq\mathsf{nsb}(\ell_{1})+\mathsf{nsb_{e}}(\ell_{1})+\mathsf{nsb_{e}}(\ell_{2})-2,\ \mathsf{nsb_{e}}(\ell)\geq\mathsf{nsb}(\ell_{2})+\mathsf{nsb_{e}}(\ell_{2})+\mathsf{nsb_{e}}(\ell_{1})-2\end{array}\right\}\end{array}

\begin{array}[]{c}\mathcal{E}^{\prime}[e_{1}^{\ell_{1}}\div^{\ell}e_{2}^{\ell_{2}}]\varrho=\mathcal{E}^{\prime}[e_{1}^{\ell_{1}}]\varrho\ \cup\ \mathcal{E}^{\prime}[e_{2}^{\ell_{2}}]\varrho\quad(\textsc{Div}^{\prime})\\ \cup\\ \left\{\begin{array}[]{c}\mathsf{nsb_{e}}(\ell)\geq\mathsf{nsb}(\ell_{1})+\mathsf{nsb_{e}}(\ell_{1})+\mathsf{nsb_{e}}(\ell_{2})-2,\ \mathsf{nsb_{e}}(\ell)\geq\mathsf{nsb}(\ell_{2})+\mathsf{nsb_{e}}(\ell_{2})+\mathsf{nsb_{e}}(\ell_{1})-2\end{array}\right\}\end{array}

\mathcal{E}^{\prime}\left[\sqrt{e^{\ell_{1}}}^{\ell}\right]\varrho=\mathcal{E}^{\prime}[e_{1}^{\ell_{1}}]\varrho\ \cup\ \big{\{}\mathsf{nsb_{e}}(\ell)\geq\mathsf{nsb_{e}}(\ell_{1})\big{\}}\quad(\textsc{Sqrt}^{\prime})

\mathcal{E}^{\prime}\left[\phi\big{(}{e^{\ell_{1}}}\big{)}^{\ell}\right]\varrho=\mathcal{E}^{\prime}[e_{1}^{\ell_{1}}]\varrho\ \cup\ \big{\{}\mathsf{nsb_{e}}(\ell)\geq+\infty\big{\}}\ \text{with}\ \phi\in\{\sin,\cos,\tan,\log,\ldots\}\quad(\textsc{Math}^{\prime})

\begin{array}[]{c}\mathcal{C}^{\prime}\left[x\texttt{:=}^{\ell}e^{\ell_{1}}\right]\varrho=\big{(}C,\varrho\left[x\mapsto\ell\right]\big{)}\ \text{where}\ C=\mathcal{E}^{\prime}[e_{1}^{\ell_{1}}]\varrho\cup\left\{\mathsf{nsb_{e}}(\ell_{1})\geq\mathsf{nsb_{e}}(\ell)\right\}\end{array}\quad(\textsc{Assign}^{\prime})

\mathcal{C}^{\prime}\left[c_{1}^{\ell_{1}}\texttt{;}c_{2}^{\ell_{2}}\right]\varrho=\big{(}C_{1}\cup C_{2},\varrho_{2}\big{)}\ \text{with}\ \big{(}C_{1},\varrho_{1}\big{)}=\mathcal{C}^{\prime}\left[c_{1}^{\ell_{1}}\right]\varrho\ \text{and}\ \big{(}C_{2},\varrho_{2}\big{)}=\mathcal{C}^{\prime}\left[c_{2}^{\ell_{2}}\right]\varrho_{1}\quad(\textsc{Seq}^{\prime})

\begin{array}[]{c}\mathcal{C}^{\prime}[\texttt{if}^{\ell}\ e^{\ell_{0}}\ \texttt{then}\ c^{\ell_{1}}\ \texttt{else}\ c^{\ell_{2}}]\ \varrho=(C_{1}\cup C_{2}\cup C_{3},\varrho^{\prime})\\ \text{where}\ \left|\begin{array}[]{l}\forall x\in\text{Id},\ \varrho^{\prime}(x)=\ell,\ (C_{1},\varrho_{1})=\mathcal{C}^{\prime}[c_{1}^{\ell_{1}}]\ \varrho,\ (C_{2},\varrho_{2})=\mathcal{C}^{\prime}[c_{2}^{\ell_{2}}]\ \varrho,\\ C_{3}=\underset{x\in\text{Id}}{\bigcup}\left\{\mathsf{nsb_{e}}(\varrho_{1}(x))\geq\mathsf{nsb_{e}}(\ell),\ \mathsf{nsb_{e}}(\varrho_{2}(x))\geq\mathsf{nsb_{e}}(\ell)\right\}\end{array}\right.\end{array}\quad(\textsc{Cond}^{\prime})

\begin{array}[]{c}\mathcal{C}^{\prime}[\texttt{while}^{\ell}\ e^{\ell_{0}}\ \texttt{do}\ c^{\ell_{1}}]\ \varrho=(C_{1}\cup C_{2},\varrho^{\prime})\\ \text{where}\ \left|\begin{array}[]{l}\forall x\in\text{Id},\ \varrho^{\prime}(x)=\ell,\ (C_{1},\varrho_{1})=\mathcal{C}^{\prime}[c_{1}^{\ell_{1}}]\ \varrho^{\prime}\\ C_{2}=\underset{x\in\text{Id}}{\bigcup}\left\{\mathsf{nsb_{e}}(\varrho(x))\geq\mathsf{nsb_{e}}(\ell),\ \mathsf{nsb_{e}}(\varrho_{1}(x))\geq\mathsf{nsb_{e}}(\ell)\right\}\end{array}\right.\end{array}\quad(\textsc{While}^{\prime})

\mathcal{C}^{\prime}[\texttt{require\_nsb}(x,\texttt{p})^{\ell}]\varrho=\emptyset\quad(\textsc{Req}^{\prime})

$\xi(\ell)(\ell_{1},\ell_{2})=\min\left(\begin{array}[]{l}\max\big{(}\mathsf{ufp}(\ell_{2})-\mathsf{ufp}(\ell_{1})+\mathsf{nsb}(\ell_{1})-\mathsf{nsb}(\ell_{2})-\mathsf{nsb_{e}}(\ell_{2}),0\big{)},\\ \max\big{(}\mathsf{ufp}(\ell_{1})-\mathsf{ufp}(\ell_{2})+\mathsf{nsb}(\ell_{2})-\mathsf{nsb}(\ell_{1})-\mathsf{nsb_{e}}(\ell_{1}),0\big{)},1\end{array}\right)$

Figure 5: Constraints solved by PI with

\min

and

\max

carry bit formulation.

3.3 Policy Iteration for Optimized Carry Bit Propagation

The policy iterations algorithm is used to solve nonlinear fixpoint equations when the function is written as the infimum of functions for which a fixpoint can be easily computed. The infimum formulation makes the function not being differentiable in the classical sense. The one proposed in [7] to solve smallest fixpoint equations in static analysis requires the fact that the function is order-preserving to ensure the decrease of the intermediate solutions provided by the algorithm. In this article, because of the nature of the semantics, we propose a policy iterations algorithm for a non order-preserving function.

More precisely, let $F$ be a map from a complete $L$ to itself such that $F=\inf_{\pi\in\Pi}f^{\pi}$ . Classical policy iterations solve $F(\mathbf{x})=\mathbf{x}$ by generating a sequence $(\mathbf{x}^{k})_{k}$ such that $f^{\pi^{k}}(\mathbf{x}^{k})=\mathbf{x}^{k}$ and $\mathbf{x}^{k+1}<\mathbf{x}^{k}$ . The set $\Pi$ is called the set of policies and $f^{\pi}$ a policy map (associated to $\pi$ ). The set of policy maps has to satisfy the selection property meaning that for all $\mathbf{x}\in L$ , there exists $\pi\in\Pi$ such that $F(\mathbf{x})=f^{\pi}(\mathbf{x})$ . This is exactly the same as for each $\mathbf{x}\in L$ , the minimization problem $\operatorname*{Min}_{\pi\in\Pi}f^{\pi}(\mathbf{x})$ has an optimal solution. If $\Pi$ is finite and $F$ is order-preserving, policy iterations converge in finite time to a fixpoint of $F$ . The number of iterations is bounded from above by the number of policies. Indeed, a policy cannot be selected twice in the running of the algorithm. This is implied by the fact that the smallest fixpoint of a policy map is computed. In this article, we adapt policy iterations to the problem of precision tuning. The function $F$ here is constructed from inequalities depicted in Figure 4 and Figure 5. We thus have naturally constraints of the form $F(\mathbf{x})\leq\mathbf{x}$ . We will give details about the construction of $F$ at Proposition 1. Consequently, we are interested in solving:

\small\begin{array}[]{lll}\displaystyle{\operatorname*{Min}_{\mathsf{nsb},\mathsf{nsb_{e}}}}&&\displaystyle{\sum_{\ell}\mathsf{nsb}(\ell)}\\ &\operatorname*{s.t.}&F\left(\begin{array}[]{c}\mathsf{nsb}\\ \mathsf{nsb_{e}}\end{array}\right)\leq\left(\begin{array}[]{c}\mathsf{nsb}\\ \mathsf{nsb_{e}}\end{array}\right)\quad\mathsf{nsb}\in\mathbb{N}^{Lab},\ \mathsf{nsb_{e}}\in\mathbb{N}^{Lab}\end{array}

(18)

Let $\xi:Lab\to\{0,1\}$ . We will write $S_{\xi}^{1}$ the system of inequalities depicted in Figure 4 and $S_{\xi}^{2}$ the system of inequalities presented at Figure 5. Note that the final system of inequalities is $S_{\xi}=S_{\xi}^{1}\cup S_{\xi}^{2}$ meaning that we add new constraints to $S_{\xi}^{1}$ . If the system $S_{\xi}^{1}$ is used alone, $\xi$ is the constant function equal to 1. Otherwise, $\xi$ is defined by the formula at the end of Figure 5.

Proposition 1

The following results hold:

1.

Let $\xi$ the constant function equal to 1. The system $S_{\xi}^{1}$ can be rewritten as $\{\mathsf{nsb}\in\mathbb{N}^{Lab}\mid F(\mathsf{nsb})\leq(\mathsf{nsb})\}$ where $F$ maps $\mathbb{R}^{Lab}\times\mathbb{R}^{Lab}$ to itself, $F(\mathbb{N}^{Lab}\times\mathbb{N}^{Lab})\subseteq(\mathbb{N}^{Lab}\times\mathbb{N}^{Lab})$ and has coordinates which are the maximum of a finite family of affine order-preserving functions.
2.

Let $\xi$ the function such that $\xi(\ell)$ equals the function defined at Figure 5. The system $S_{\xi}$ can be rewritten as $\{(\mathsf{nsb},\mathsf{nsb_{e}})\in\mathbb{N}^{Lab}\times\mathbb{N}^{Lab}\mid F(\mathsf{nsb},\mathsf{nsb_{e}})\leq(\mathsf{nsb},\mathsf{nsb_{e}})\}$ where $F$ maps $\mathbb{R}^{Lab}\times\mathbb{R}^{Lab}$ to itself, $F(\mathbb{N}^{Lab}\times\mathbb{N}^{Lab})\subseteq(\mathbb{N}^{Lab}\times\mathbb{N}^{Lab})$ and all its coordinates are the min-max of a finite family of affine functions.

From Proposition 1, when $S_{\xi}$ is used, we can write $F$ as $F=\min_{\pi\in\Pi}f^{\pi}$ , where $f^{\pi}$ is the maximum of a finite family of affine functions and thus used a modified policy iterations algorithm. The set of policies here is a map $\pi:Lab\mapsto\{0,1\}$ . A choice is thus a vector of 0 or 1. A policy map $f^{\pi}$ is a function $\mathbb{N}^{Lab}$ to itself such that the coordinates are $f_{\ell}^{\pi}(\ell)$ . If the coordinate $f_{\ell}^{\pi}(\ell)$ depends on $\xi$ then $\xi(\ell)=\pi(\ell)$ . Otherwise, the function is the maximum of affine functions and a choice is not required.

Result: An over-approximation of an optimal solution of Equation (18)

1 Let

k:=0

S:=+\infty

;

2 Choose

\pi^{0}\in\Pi

;

3 Select an optimal solution of

(\mathsf{nsb}^{k},\mathsf{nsb_{e}}^{k})

the integer linear program:

\operatorname*{Min}\left\{\sum_{\ell\in Lab}\mathsf{nsb}(\ell)\mid f^{\pi^{k}}(\mathsf{nsb},\mathsf{nsb_{e}})\leq(\mathsf{nsb},\mathsf{nsb_{e}}),\ \mathsf{nsb}\in\mathbb{N}^{Lab},\ \mathsf{nsb_{e}}\in\mathbb{N}^{Lab}\right\}\enspace;

if $\sum_{\ell\in Lab}\mathsf{nsb}^{k}(\ell)<S$ then

S:=\sum_{\ell\in Lab}\mathsf{nsb}^{k}(\ell)

;

5 Choose

\pi^{k+1}\in\Pi

such that

F(\mathsf{nsb}^{k},\mathsf{nsb_{e}}^{k})=f^{\pi^{k+1}}(\mathsf{nsb}^{k},\mathsf{nsb_{e}}^{k})

;

k:=k+1

and go to 3;

8else

9 Return

S

and

\mathsf{nsb}^{k}

10 end if

Algorithm 1 Non-monotone Policy Iterations Algorithm

Proposition 2 (Algorithm correctness)

The sequence $(\sum_{\ell\in Lab}\mathsf{nsb}^{k}(\ell))_{0\leq k\leq K}$ generated by Algorithm 1 is of finite length (i.e. $K\in\mathbb{N}$ ) and satisfies a strict decrease before convergence: $\sum_{\ell\in Lab}\mathsf{nsb}^{k+1}(\ell)<\sum_{\ell\in Lab}\mathsf{nsb}^{k}(\ell)$ if $k<K-1$ and $\sum_{\ell\in Lab}\mathsf{nsb}^{K}(\ell)=\sum_{\ell\in Lab}\mathsf{nsb}^{K-1}(\ell)$ . The number of terms is smaller than the number of policies.

Proof

Let $\sum_{\ell\in Lab}\mathsf{nsb}^{k}(\ell)$ be a term of the sequence and $(\mathsf{nsb}^{k},\mathsf{nsb_{e}}^{k})$ be the optimal solution of $\operatorname*{Min}\{\sum_{\ell\in Lab}\mathsf{nsb}(\ell)\mid f^{\pi^{k}}(\mathsf{nsb},\mathsf{nsb_{e}})\leq(\mathsf{nsb},\mathsf{nsb_{e}}),\ \mathsf{nsb}\in\mathbb{N}^{Lab},\ \mathsf{nsb_{e}}\in\mathbb{N}^{Lab}\}$ . Then $F(\mathsf{nsb}^{k},\mathsf{nsb_{e}}^{k})\leq f^{\pi^{k}}(\mathsf{nsb}^{k},\mathsf{nsb_{e}}^{k})$ by definition of $F$ . Moreover,
$F(\mathsf{nsb}^{k},\mathsf{nsb_{e}}^{k})=f^{\pi^{k+1}}(\mathsf{nsb}^{k},\mathsf{nsb_{e}}^{k})$ and $f^{\pi^{k}}(\mathsf{nsb}^{k},\mathsf{nsb_{e}}^{k})\leq(\mathsf{nsb}^{k},\mathsf{nsb_{e}}^{k})$ . It follows that $f^{\pi^{k+1}}(\mathsf{nsb}^{k},\mathsf{nsb_{e}}^{k})\leq(\mathsf{nsb}^{k},\mathsf{nsb_{e}}^{k})$ and $(\mathsf{nsb}^{k},\mathsf{nsb_{e}}^{k})$ is feasible for the minimisation problem for which $(\mathsf{nsb}^{k+1},\mathsf{nsb_{e}}^{k+1})$ is an optimal solution. We conclude that $\sum_{\ell\in Lab}\mathsf{nsb}^{k+1}(\ell)\leq\sum_{\ell\in Lab}\mathsf{nsb}^{k}(\ell)$ and the Algorithm terminates if the equality holds or continues as the criterion strictly decreases. Finally, from the strict decrease, a policy cannot be selected twice without terminating the algorithm. In conclusion, the number of iterations is smaller than the number of policies. $\Box$

Figure 5 displays the new rules that we add to the global system of constraints in the case where we optimize the carry bit of the elementary operations. Before introducing the optimized function $\xi$ , the definition of the unit in the last place $\mathsf{ulp}$ of a real number $x$ is defined in Equation(19).

\small\mathsf{ulp}(x)=\mathsf{ufp}(x)-\mathsf{nsb}(x)+1\enspace.

(19)

Indeed, during an operation between two numbers, the $\xi$ function is computed as follows: If the $\mathsf{ulp}$ of one of the two operands is greater than the $\mathsf{ufp}$ of the other one (or conversely) then the two numbers are not aligned and no carry bit can be propagated through the operation (otherwise $\xi=1$ ). This idea is presented in Figure 5 in which $\xi$ is formulated by $\min$ and $\max$ operators. Winning one bit may seen a ridiculous optimization at first sight. However, when many operations are done in a program which has to compute with some tens of $\mathsf{nsb}$ , this is far from being negligible. Formally, let $x$ and $y$ be the operands of some operation whose result is $z$ . The errors are $\varepsilon(x)$ , $\varepsilon(y)$ and $\varepsilon(z)$ respectively on $x$ , $y$ and $z$ . Note that the behaviour of $\xi$ is also valid in the case of propagation of the carry bit on the errors. Consequently, using the definitions of $\mathsf{ufp_{e}}$ and $\mathsf{nsb_{e}}$ introduced in Section 3.1, we have $\mathsf{ufp_{e}}(x)=\mathsf{ufp}(x)-\mathsf{nsb}(x)$ and $\mathsf{ulp_{e}}(x)=\mathsf{ufp_{e}}(x)-\mathsf{nsb_{e}}(x)+1$ with $\mathsf{ufp_{e}}$ and $\mathsf{ulp_{e}}$ are the unit in the first and in the last place of the errors, respectively, for $x$ (same reasoning for $y$ and $z$ ). The optimized function $\xi$ of Figure 5 derives from Equation (20).

\small\xi(z,x,y)=\begin{cases}0&\mathsf{ufp_{e}}(x)-\mathsf{nsb_{e}}(x)\geq\mathsf{ufp}(y)-\mathsf{nsb}(y)~{}~{}or~{}~{}\text{conversely},\\ 1&\text{otherwise.}\end{cases}

(20)

As mentioned in Equation (20), to compute the $\mathsf{ulp}$ of the errors on the operands, we need to estimate the number of bits of the error $\mathsf{nsb_{e}}$ for each operand which explains the new rules and constraints given in Figure 5. Hence, these constraints are complementary to the rules of Figure 4, already explained in Section 3.2, in which the only difference is that we activate the new function $\xi$ instead of its over-approximation of Figure 4. Let us concentrate on the rules of Figure 5. The function $\mathcal{E}^{\prime}[e]\ \varrho$ generates the new set of constraints for an expression $e\in Expr$ in the environment $\varrho$ . For Rule $(\textsc{Const}^{\prime})$ , the number of significant bits of the error $\mathsf{nsb_{e}}=0$ whereas we impose that the $\mathsf{nsb_{e}}$ of a variable $x$ at control point $\ell$ is less than the last assignment of $\mathsf{nsb_{e}}$ in $\varrho(x)$ as shown in Rule $(\textsc{Id}^{\prime})$ of Figure 5. Considering Rule $(\textsc{Add}^{\prime})$ , we start by generating the new set of constraints $\mathcal{E}^{\prime}[e_{1}^{\ell_{1}}]\varrho$ and $\mathcal{E}^{\prime}[e_{2}^{\ell_{2}}]\varrho$ on the operands at control points $\ell_{1}$ and $\ell_{2}$ . Then, we require that $\mathsf{nsb_{e}}(\ell)\geq\mathsf{nsb_{e}}(\ell_{1})$ and $\mathsf{nsb_{e}}(\ell)\geq\mathsf{nsb_{e}}(\ell_{2})$ where the result of the addition is stored at control point $\ell$ . Additionally, the number of significant bits of the error on the result $\mathsf{nsb_{e}}(\ell)$ at control point $\ell$ is computed as shown hereafter.

\small\mathsf{nsb_{e}}(\ell)\geq\max\left(\begin{array}[]{l}\mathsf{ufp}(\ell_{1})-\mathsf{nsb}(\ell_{1})\\ \mathsf{ufp}(\ell_{2})-\mathsf{nsb}(\ell_{2})\end{array}\right)-\min\left(\begin{array}[]{l}\mathsf{ufp}(\ell_{1})-\mathsf{nsb}(\ell_{1})-\mathsf{nsb_{e}}(\ell_{1})\\ \mathsf{ufp}(\ell_{2})-\mathsf{nsb}(\ell_{2})-\mathsf{nsb_{e}}(\ell_{2})\end{array}\right)+\xi(\ell)(\ell_{1},\ell_{2})

By breaking the $\min$ and $\max$ operators, we obtain the constraints on $\mathsf{nsb_{e}}(\ell)$ of Rule $(\textsc{Add}^{\prime})$ . For the subtraction, the constraints generated are similar to the addition case. Considering now Rule $(\textsc{Mult}^{\prime})$ , as we have defined in Section 3.2, $\varepsilon(z)=x\cdot\varepsilon(y)+y\cdot\varepsilon(x)+\varepsilon(x)\cdot\varepsilon(y)$ where $z$ is the result of the product of $x$ and $y$ . By reasoning on the $\mathsf{ulp}$ of the error, we bound $\varepsilon(z)$ by

	$\displaystyle\varepsilon(z)=2^{\mathsf{ufp}(x)}\cdot 2^{\mathsf{ufp}(y)-\mathsf{nsb}(y)-\mathsf{nsb_{e}}(y)+1}+2^{\mathsf{ufp}(y)}\cdot 2^{\mathsf{ufp}(x)-\mathsf{nsb}(x)-\mathsf{nsb_{e}}(x)+1}$
	$\displaystyle+2^{\mathsf{ufp}(y)+\mathsf{ufp}(x)-\mathsf{nsb}(x)-\mathsf{nsb}(y)-\mathsf{nsb_{e}}(x)-\mathsf{nsb_{e}}(y)+2}$

By selecting the smallest term $\mathsf{ufp}(y)+\mathsf{ufp}(x)-\mathsf{nsb}(x)-\mathsf{nsb}(y)-\mathsf{nsb_{e}}(x)-\mathsf{nsb_{e}}(y)+2$ , we obtain that

\small\mathsf{nsb_{e}}(\ell)\geq\max\left(\begin{array}[]{l}\mathsf{ufp}(\ell_{1})+\mathsf{ufp}(\ell_{2})-\mathsf{nsb}(\ell_{1})\\ \mathsf{ufp}(\ell_{1})+\mathsf{ufp}(\ell_{2})-\mathsf{nsb}(\ell_{2})\end{array}\right)-\begin{array}[]{l}\mathsf{ufp}(\ell_{1})+\mathsf{ufp}(\ell_{2})-\mathsf{nsb}(\ell_{1})-\\ \mathsf{nsb}(\ell_{2})-\mathsf{nsb_{e}}(\ell_{1})-\mathsf{nsb_{e}}(\ell_{2})\\ +2\end{array}\enspace.

Finally, by simplifying the equation above we found the constraints of Rule $(\textsc{Mult}^{\prime})$ in Figure 5 (same for Rule $(\textsc{Div}^{\prime})$ ). For Rule $(\textsc{Sqrt}^{\prime})$ , we generate the constraints on the expression $\mathcal{E}^{\prime}[e_{1}^{\ell_{1}}]\varrho$ and we require that $\mathsf{nsb_{e}}$ of the result stored at control point $\ell$ is greater than the $\mathsf{nsb_{e}}$ of the expression a control point $\ell_{1}$ . For Rule $(\textsc{Math}^{\prime})$ , we assume that $\mathsf{nsb_{e}}(\ell)$ is unbounded. Concerning the commands, we define the set $\mathcal{C}^{\prime}[c]\ \varrho$ which has the same function as $\mathcal{C}$ defined in Figure 4. The reasoning on the commands also remains similar except that this time we reason on the number of bits of the errors $\mathsf{nsb_{e}}$ . The only difference is in Rule $(\textsc{Req}^{\prime})$ where the set of constraints is empty. Let us recall that the constraints of Figure 5 are added to the former constraints of Figure 4 and are sent to a linear solver (GLPK in practice).

\frac{\varrho(x)=c\texttt{\#p}}{\langle x^{\ell},\varrho\rangle\longrightarrow\langle c^{\ell}\texttt{\#p},\varrho\rangle}

\frac{c=c_{1}\odot c_{2},\ \texttt{p}=\mathsf{ufp}(c)-\mathsf{ufp_{e}}\big{(}c^{\ell}\texttt{\#p}\big{)}}{\langle c_{1}^{\ell_{1}}\texttt{\#$p_{1}$}\odot^{\ell}c_{2}^{\ell_{2}}\texttt{\#$p_{2}$},\varrho\rangle\longrightarrow\langle c\texttt{\#p},\varrho\rangle}\quad\odot\in\{+,-,\times,\div\}

\frac{\langle e_{1}^{\ell_{1}},\varrho\rangle\longrightarrow\langle e_{1}^{\prime\ell_{1}},\varrho\rangle}{\langle e_{1}^{\ell_{1}}\odot^{\ell}e_{2}^{\ell_{2}},\varrho\rangle\longrightarrow\langle e_{1}^{\prime\ell_{1}}\odot^{\ell}e_{2}^{\ell_{2}},\varrho\rangle}\qquad\frac{\langle e_{2}^{\ell_{2}},\varrho\rangle\longrightarrow\langle e_{2}^{\prime\ell_{2}},\varrho\rangle}{\langle c_{1}^{\ell_{1}}\texttt{\#p}\odot^{\ell}e_{2}^{\ell_{2}},\varrho\rangle\longrightarrow\langle c_{1}^{\ell_{1}}\sharp\texttt{p}\odot^{\ell}e_{2}^{\prime\ell_{2}},\varrho\rangle}

\frac{\langle e^{\ell_{1}},\varrho\rangle\longrightarrow\langle e^{\prime\ell_{1}},\varrho\rangle}{\langle\phi(e^{\ell_{1}})^{\ell},\varrho\rangle\longrightarrow\langle\phi(e^{\prime\ell_{1}})^{\ell},\varrho\rangle}\qquad\frac{c=\phi(c_{1})\quad\texttt{q}=\mathsf{ufp}(c)-\mathsf{ufp_{e}}(c^{\ell}\texttt{\#p})}{\langle\phi(c_{1}^{\ell_{1}}\texttt{\#p})^{\ell},\varrho\rangle\longrightarrow\langle c^{\ell}\texttt{\#q},\varrho\rangle}\quad\phi\in\{\sin,\cos,\ldots\}

\frac{\langle e^{\ell_{1}},\varrho\rangle\longrightarrow\langle e^{\prime\ell_{1}},\varrho\rangle}{\langle\sqrt{e^{\ell_{1}}}^{\ell},\varrho\rangle\longrightarrow\langle\sqrt{e^{\prime\ell_{1}}}^{\ell},\varrho\rangle}\qquad\frac{c=\sqrt{c_{1}}\quad\texttt{q}=\mathsf{ufp}(c)-\mathsf{ufp_{e}}(c^{\ell}\texttt{\#p})}{\langle\sqrt{c_{1}^{\ell_{1}}\texttt{\#p}}^{\ell},\varrho\rangle\longrightarrow\langle c^{\ell}\texttt{\#q},\varrho\rangle}

Figure 6: Small Step Operational semantics of arithmetic expressions.

4 Correctness

In this section, we present proofs of correctness concerning the soundness of the analysis (Section 4.1) and the integer nature of the solutions (Section 4.2).

4.1 Soundness of the Constraint System

Let $\equiv$ denote the syntactic equivalence and let $e^{\ell}\in Expr$ be an expression. We write $\mathsf{Const}(e^{\ell})$ the set of constants occurring in the expression $e^{\ell}$ . For example, $\mathsf{Const}(18.0^{\ell_{1}}\times^{\ell_{2}}x^{\ell_{3}}+^{\ell_{4}}12.0^{\ell_{5}}\times^{\ell_{6}}y^{\ell_{7}}+^{\ell_{8}}z^{\ell_{9}})=\{18.0^{\ell_{1}},12.0^{\ell_{5}}\}$ . Also, we denote by $\tau\ :\ \text{Lab}\rightarrow\mathbb{N}$ a function mapping the labels of an expression to a $\mathsf{nsb}$ . The notation $\tau\models\mathcal{E}[e^{\ell}]\varrho$ means that $\tau$ is the minimal solution to the ILP $\mathcal{E}[e^{\ell}]\varrho$ . We write $\varrho_{\perp}$ the empty environment ( $dom(\varrho_{\perp})=\emptyset$ ).

The small step operational semantics of our language is displayed in Figure 6. It is standard, the only originality being to indicate explicitly the $\mathsf{nsb}$ of constants. For the result of an elementary operation, this $\mathsf{nsb}$ is computed in function of the $\mathsf{nsb}$ of the operands. Lemma 1 below asses the soundness of the constraints for one step of the semantics.

Lemma 1

Given an expression $e^{\ell}\in Expr$ , if $e^{\ell}\rightarrow e^{\prime\ell}$ and $\tau\models\mathcal{E}[e^{\ell}]\varrho_{\perp}$ then for all $c^{\ell_{c}}$ #p $\in$ $\mathsf{Const}(e^{\prime\ell})$ we have p = $\tau(\ell_{c})$ .

Proof

By case examination of the rules of Figure 4. Hereafter, we focus on the most interesting case of addition of two constants. Recall that $\mathsf{ufp_{e}}(\ell)=\mathsf{ufp}(\ell)-\mathsf{nsb}(\ell)$ for any control point $\ell$ . Assuming that $e^{\ell}\equiv c_{1}^{\ell_{1}}+^{\ell}c_{2}^{\ell_{2}}$ then by following the reduction rule of Figure 6, we have $e^{\ell}\rightarrow c^{\ell}\texttt{\#p}$ with $\texttt{p}=\mathsf{ufp}(c)-\mathsf{ufp_{e}}\big{(}c\big{)}$ . On the other side, by following the set of constraints of Rule (Add) in Figure 4 we have $\small\mathcal{E}[e^{\ell}]\varrho=\{\mathsf{nsb}(\ell_{1})\geq\mathsf{nsb}(\ell)+\mathsf{ufp}(\ell_{1})-\mathsf{ufp}(\ell)+\xi(\ell)(\ell_{1},\ell_{2}),\hskip 2.84544pt\mathsf{nsb}(\ell_{2})\geq\mathsf{nsb}(\ell)+\mathsf{ufp}(\ell_{2})$ $-\mathsf{ufp}(\ell)+\xi(\ell)(\ell_{1},\ell_{2})\}$ . These constraints can be written as

	$\displaystyle\small\mathsf{nsb}(\ell)$	$\displaystyle\leq$	$\displaystyle\mathsf{ufp}(\ell)-\mathsf{ufp}(\ell_{1})+\mathsf{nsb}(\ell_{1})-\xi(\ell)(\ell_{1},\ell_{2})$
	$\displaystyle\mathsf{nsb}(\ell)$	$\displaystyle\leq$	$\displaystyle\mathsf{ufp}(\ell)-\mathsf{ufp}(\ell_{2})+\mathsf{nsb}(\ell_{2})-\xi(\ell)(\ell_{1},\ell_{2})$

and may themselves be rewritten as Equation (16), i.e.

\small\mathsf{nsb}(\ell)\leq\mathsf{ufp}(\ell)-\max\big{(}\mathsf{ufp}(\ell_{1})-\mathsf{nsb}(\ell_{1}),\mathsf{ufp}(\ell_{2})-\mathsf{nsb}(\ell_{2})\big{)}-\xi(\ell)(\ell_{1},\ell_{2})\enspace.

Since, obviously, $\mathsf{ufp}(c)=\mathsf{ufp}(\ell)$ and since the solver finds the minimal solution to the ILP, it remains to show that

\small\mathsf{ufp_{e}}(\ell)=\max\big{(}\mathsf{ufp}(\ell_{1})-\mathsf{nsb}(\ell_{1}),\mathsf{ufp}(\ell_{2})-\mathsf{nsb}(\ell_{2}),\mathsf{ufp}(\ell)-\mathsf{prec}(\ell)\big{)}-\xi(\ell)(\ell_{1},\ell_{2})

which corresponds to the assertion of Equation(15). Consequently, $\mathsf{nsb}(\ell)=\texttt{p}$ as required, for this case, in Figure 6. $\Box$

Theorem 4.1

Given an expression $e^{\ell}\rightarrow e^{\prime\ell}$ . If $e^{\ell}\rightarrow^{*}e^{\prime\ell}$ and if $\tau\models\mathcal{E}[e^{\ell}]\varrho_{\perp}$ , then $\forall$ $c^{\ell_{c}}\texttt{\#p}$ $\in$ $\mathsf{Const}(e^{\prime\ell})$ we have $\texttt{p}=\tau(\ell_{c})$ .

Proof

By recurrence on the length of the reduction path. $\Box$

4.2 ILP Nature of the Problem

In this section, we give insights about the complexity of the problem. The computation relies on integer linear programming. Integer linear programming is known to belong to the class of NP-Hard problems. A lower bound of the optimal value in a minimization problem can be furnished by the continuous linear programming relaxation. This relaxation is obtained by removing the integrity constraint. Recall that a (classical) linear program can be solved in polynomial-time. Then, we can solve our problem in polynomial-time if we can show that the continuous linear programming relaxation of our ILP has an unique optimal solution with integral coordinates. Proposition 3 presents a situation where a linear program has a unique optimal solution which is a vector of integers.

Proposition 3

Let $G:[0,+\infty)^{d}\mapsto[0,+\infty)^{d}$ be an order-preserving function such that $G(\mathbb{N}^{d})\subseteq\mathbb{N}^{d}$ . Suppose that the set $\{y\in\mathbb{N}^{d}\mid G(y)\leq y\}$ is non-empty. Let $\varphi:\mathbb{R}^{d}\mapsto\mathbb{R}$ a strictly monotone function such that $\varphi(\mathbb{N}^{d})\subseteq\mathbb{N}$ . Then, the minimization problem:

\displaystyle{\operatorname*{Min}_{y\in[0,+\infty)^{d}}}\ \varphi(y)\ \operatorname*{s.t.}G(y)\leq y

has an unique optimal solution which is integral.

Proof

Let $L:=\{x\in[0,+\infty)^{d}\mid G(x)\leq x\}$ and $u=\inf L$ . It suffices to prove that $u\in\mathbb{N}^{d}$ . Indeed, as $\varphi$ is stricly monotone then $\varphi(u)<\varphi(x)$ for all $x\in[0,+\infty)^{d}$ s.t. $G(x)\leq x$ and $x\neq u$ . The optimal solution is thus $u$ . If $u=0$ , the result holds. Now suppose that $0<u$ , then $0\leq G(0)$ . Let $M:=\{y\in\mathbb{N}^{d}\mid y\leq G(y),y\leq u\}$ . Then $0\in M$ and we write $v:=\sup M$ . As $M$ is a complete lattice s.t. $G(M)\subseteq M$ , from Tarski’s theorem, $v$ satisfies $G(v)=v$ and $v\leq u$ . Moreover, $v\in\mathbb{N}^{d}$ and $v\leq u$ . Again, from Tarski’s theorem, $u$ is the smallest fixpoint of $G$ then it coincides with $v$ . We conclude that $u\in\mathbb{N}^{d}$ . $\Box$

Theorem 4.2

Assume that the system of inequalities depicted in Figure 4 has a solution. The smallest amount of memory $\sum_{\ell\in Lab}\mathsf{nsb}(\ell)$ for the system of inequalities depicted in Figure 4 can be computed in polynomial-time by linear programming.

Proof

The function $\sum_{\ell\in Lab}\mathsf{nsb}(\ell)$ is strictly monotone and stable on integers. From the first statement of Proposition 1, the system of constraints is of the form $F(\mathsf{nsb})\leq\mathsf{nsb}$ where $F$ is order-preserving and stable on integers. By assumption, there exists an vector of integers $\mathsf{nsb}$ s.t. $F(\mathsf{nsb})\leq\mathsf{nsb}$ . We conclude from Proposition 3. $\Box$

For the second system, in practice, we get integral solutions to the continuous linear programming relaxation of our ILP of Equation (18). However, because of the lack of monotonicity of the functions for the rules (ADD) and (SUB), we cannot exploit Proposition 3 to prove the polynomial-time solvability.

5 Experimental Results

Program TH BL IEEE ILP-time BL IEEE PI-time H S D LD $10^{-4}$ 61% 43% 0.9s 62% 45% 1.5s 8 88 25 0 $10^{-6}$ 50% 21% 0.9s 51% 21% 1.4s 2 45 74 0 arclength $10^{-8}$ 37% 3% 0.8s 38% 4% 1.6s 2 6 113 0 $10^{-10}$ 24% -1% 1.0s 25% -1% 1.7s 2 0 116 3 $10^{-12}$ 12% -17% 0.3s 14% -8% 1.5s 2 0 109 10 $10^{-4}$ 64% 45% 0.1s 67% 56% 0.5s 6 42 1 0 $10^{-6}$ 53% 30% 0.2s 56% 31% 0.5s 1 27 21 0 simpson $10^{-8}$ 40% 4% 0.1s 43% 7% 0.3s 1 5 43 0 $10^{-10}$ 27% 1% 0.1s 28% 1% 0.4s 1 0 48 0 $10^{-12}$ 16% 1% 0.1s 16% 1% 0.3s 0 1 48 0 $10^{-4}$ 73% 61% 0.2s 76% 62% 1.0s 53 69 0 0 $10^{-6}$ 62% 55% 0.2s 65% 55% 1.0s 2 102 0 0 accelerometer $10^{-8}$ 49% 15% 0.2s 52% 18% 1.0s 2 33 69 0 $10^{-10}$ 36% 1% 0.2s 39% 1% 1.0s 2 0 102 0 $10^{-12}$ 25% 1% 0.2s 28% 1% 1.0s 2 0 102 0 $10^{-4}$ 78% 66% 0.08s 79% 68% 1.3s 46 38 0 0 $10^{-6}$ 67% 53% 0.08s 68% 56% 0.5s 12 70 2 0 rotation $10^{-8}$ 53% 29% 0.07s 54% 29% 0.4s 0 46 38 0 $10^{-10}$ 40% 0% 0.1s 41% 0% 0.5s 0 0 84 0 $10^{-12}$ 29% 0% 0.09s 30% 0% 0.5s 0 0 48 0 $10^{-4}$ 68% 46% 1.8s 69% 46% 10.7s 260 581 0 0 $10^{-6}$ 57% 38% 1.8s 58% 45% 11.0s 258 580 3 0 lowPassFilter $10^{-8}$ 44% -7% 2.0s 45% -7% 11.4s 258 2 581 0 $10^{-10}$ 31% -7% 1.7s 32% -7% 10.9s 258 0 583 0 $10^{-12}$ 20% -7% 1.8s 21% -7% 11.3s 258 0 583 0 $10^{-4}$ 71% 54% 0.15s 71% 54% 0.4s 0 13 0 0 $10^{-6}$ 60% 50% 0.2s 60% 50% 0.5s 0 12 1 0 Pendulum $10^{-8}$ 47% 0% 012s 47% 0% 0.4s 0 0 13 0 $10^{-10}$ 33% 0% 0.16s 34% 0% 0.5 0 0 13 0 $10^{-12}$ 22% 0% 0.11s 22% 0% 0.4s 0 0 13 0

Table 1: Precision tuning results for XXX for the ILP and PI methods.

In this section, we aim at evaluating the performance of our tool XXX implementing the techniques of Section 3. Recall that XXX reduces the problem of precision tuning to an ILP by generating the constraints defined in Section 3.2 which can be solved by a linear solver. We use GLPK²²2https://www.gnu.org/software/glpk/ in practice with continuous variables since we know from Section 3.3 that it is possible. Indeed, GLPK’s constant col-kind is set to cv and not to iv as for integer problems. Alternatively, a second optimized set of constraints may be generated by applying the PI technique introduced in Section 3.3.

We have evaluated XXX on several numerical programs. Two of them were used as a benchmark for precision tuning in prior work [18] and are coming from the GNU scientific library (GSL): arclength which is a program first introduced in [2] and the simpson program which corresponds to an implementation of the widely used Simpson’s rule for integration in numerical analysis [13]. In addition, we have experimented the techniques introduced in this article on three programs used as benchmarks for XXX in its former version [3, 4, 5]. The program rotation performs a matrix-vector product to rotate a vector around the $z$ axis by an angle of $\theta$ [5]. The second program accelerometer is an application to measure an inclination angle with an accelerometer [3] and, finally, lowPassFilter program [4] which is taken from a pedometer application [14] that counts the number of footsteps of the user. The last two programs come from the IoT field. Also, we take the pendulum program already introduced in Section 2.

The experiments shown in Table 1 present the results of precision tuning returned by XXX for each error threshold $10^{-4}$ , $10^{-6}$ , $10^{-8}$ and $10^{-10}$ . The error threshold represents the $\mathsf{nsb}$ . For instance, the result is required to be correct up to 4 digits for an error threshold of $10^{-4}$ . Obviously, XXX counting the number of significant bits, these threshold are converted in binary. In Table 1, we represent by ”TH” the error threshold given by the user. ”BL” is the percentage of optimization at bit level. ”IEEE” denotes the percentage of optimized variables in IEEE754 formats (binary16, binary32, $\ldots$ The $\mathsf{nsb}$ obtained at bit level is approximated by the immediately the upper number of bits corresponding to a IEEE754 format). ”ILP-time” is the total time of XXX analysis in the case of ILP formulation. We have also ”PI-time” to represent the time passed by XXX to find the right policy and to resolve the precision tuning problem. ”H”, ”S”, ”D” and ”LD” denote respectively the number of variables obtained in, half, float, double and long-double precision when using the PI formulation that clearly displays better results.

Program	Tool	#Bits saved - Time in seconds
		Threshold $10^{-4}$	Threshold $10^{-6}$	Threshold $10^{-8}$	Threshold $10^{-10}$
arclength	XXX ILP (28)	2464b. - 1.8s.	2144b. - 1.5s.	1792b. - 1.7s.	1728b. - 1.8s.
	XXX SMT (22)	1488b. - 4.7s.	1472b. - 3.04s.	864b. - 3.09s.	384b. - 2.9s.
	Precimonious (9)	576b. - 146.4s.	576b. - 156.0s.	576b. - 145.8s.	576b. - 215.0s.
simpson	XXX ILP (14)	1344b. - 0.4s.	1152b. - 0.5s.	896b. - 0.4s.	896b. - 0.4s.
	XXX SMT (11)	896b. - 2.9s.	896b. - 1.9s.	704b. - 1.7s.	704b. - 1.8s.
	Precimonious (10)	704b. - 208.1s.	704b. - 213.7s.	704b. - 207.5s.	704b. - 200.3s.
rotation	XXX ILP (25)	2624b. - 0.47s.	2464b. - 0.47s.	2048b. - 0.54s.	1600b. - 0.48s.
	XXX SMT (22)	1584b. - 1.85s.	2208b. - 1.7s.	1776b. - 1.6s.	1600b. - 1.7s.
	Precimonious (27)	2400b. - 9.53s.	2592b. - 12.2s.	2464b. - 10.7s.	2464b. - 7.4s.
accel.	XXX ILP (18)	1776b. - 1.05s.	1728b. - 1.05s.	1248b. - 1.04s.	1152b. - 1.03s.
	XXX SMT (15)	1488b. - 2.6s.	1440b. - 2.6s.	1056 - 2.4s.	960b. - 2.4s.
	Precimonious (0)	-	-	-	-

Table 2: Comparison between XXX ILP, XXX SMT and Precimonious: number of bits saved by the tool and time in seconds for analyzing the programs.

Let us focus on the first ”TH”, ”BL”, ”IEEE” and ”ILP-time” columns of Table 1. These columns have been measured for the case of ILP by being pessimistic on the propagation of carry bit as described in Section 3.2. We recall that all variables of our programs are in double precision before analysis. For the arclength program, the final percentage of bits after optimization reaches $61\%$ at bit level while it achieves $43\%$ in IEEE formats ( $100\%$ is the percentage of all variables initially in double precision, $121$ variables for the original arclength program that used $7744$ bits). This is obtained in only $0.9$ second by applying the ILP formulation. As another option, when we refine the solution by applying the policy iteration method (from the sixth column), XXX displays better results attaining $62\%$ at bit level and $43\%$ for the IEEE formats. Although XXX needs more time of analysis to find and iterate between policies, the time of analysis remain negligible, not exceeding $1.5$ seconds. For a total of $121$ variables for the arclength original program, XXX succeeds in tuning $8$ variables to half precision (H), the majority with $88$ variables passes to simple precision (S) whereas $25$ variables remain in double precision (D) for an error threshold of $10^{-4}$ . We remark that our second method displays better results also for the other user error thresholds. For the simpson, accelerometer, rotation and lowPassFilter, the improvement is also more important when using the PI technique than when using the ILP formulation. For instance, for an error threshold of $10^{-6}$ for the simpson program, only one variable passes to half precision, $27$ variables turns to simple precision while $21$ remains in double precision with $56\%$ of percentage of total number of bits at bit level using the policy iteration method. Concerning the pendulum code, the two techniques return the same percentage at bit level and IEEE754 format for the majority of error thresholds (except $10^{-10}$ where XXX reaches $34\%$ at bit level when using the PI method).

Now, we stress on the negative percentage that we obtain in Table 1, especially for the arclength program with $10^{-10}$ and $10^{-12}$ for the columns IEEE and for the lowPassFilter program for errors of $10^{-8}$ , $10^{-10}$ and $10^{-12}$ . In fact, XXX being able to return new formats for any precision required by the user without additional cost nor by increasing the complexity even if it fails to have a significant improvement on the program output. To be specific, taking again the arclength program, for an error of $10^{-12}$ , XXX fulfills this requirement by informing the user that this precision is achievable only if 10 variables passes to the long double precision (LD) which is more than the original program whose variables are all in double precision. By doing so, the percentage of IEEE formats for both ILP and PI formulations reaches $-17\%$ and $-8\%$ , respectively. Same reasoning is adopted for the lowPassFilter which spends more time, nearly $12$ seconds, with the policy iteration technique to find the optimized formats (total of 841 variables). Note that in these cases, other tools like Precimonious [18] fail to propose formats.

Table 2 shows a comparison between the new version of XXX combining both ILP and PI formulations, the former version of XXX that uses the Z3 SMT solver coupled to binary search to find optimal solution [5] and the prior state-of the-art Precimonious [18]. The results of the mixed precision tuning are shown for the arclength, simpson, rotation and accelerometer programs. Since XXX and Precimonious implement two different techniques, we have adjusted the criteria of comparison in several points. First, we mention that XXX optimizes much more variables than Precimonious. While it disadvantages XXX, we only consider in the experiments of Table 2 the variables optimized by Precimonious to estimate the quality of the optimization. Second, let us note that the error thresholds are expressed in base $2$ in XXX and in base $10$ in Precimonious. For the relevance of comparisons, all the thresholds are expressed in base $10$ in tables 1 and 2. In practice, XXX will use the base $2$ threshold immediately lower than the required base $10$ threshold.

In Table 2, we indicate in bold the tool that exhibits better results for each error threshold and each program. Starting with the arclength program, XXX ILP displays better results than the other tools by optimizing $28$ variables. For an error threshold of $10^{-4}$ , $2464$ bits are saved by XXX ILP in $1.8$ seconds while XXX Z3 saved only 1488 bits in more time ( $11$ seconds). Precimonious were the slowest tool on this example, more than $2$ minutes with $576$ bits for only $9$ variables optimized. For the simpson program, XXX ILP do also better than both other tools. However, for the rotation program, XXX ILP saves more bits than the other tools only for an error of $10^{-4}$ while Precimonious do well for this program for the rest of error thresholds. Finally, Precimonious fails to tune the accelerometer program (0 variables) at the time that XXX ILP do faster (only 1 second) to save much more bits than XXX SMT for the four given error thresholds.

6 Conclusion and Perspectives

In this article, we have introduced a new technique for precision tuning, totally different from the existing ones. Instead of changing more or less randomly the data types of the numerical variables and running the programs to see what happens, we propose a semantical modelling of the propagation of the numerical errors throughout the code. This yields a system of constraints whose minimal solution gives the best tuning of the program, furthermore, in polynomial time. Two variants of this system are proposed. The first one corresponds to a pure ILP. The second one, which optimizes the propagation of carries in the elementary operations can be solved using policy iterations [7]. Proofs of correctness concerning the soundness of the analysis and the integer nature of the solutions have been presented in Section 4 and experimental results showing the efficiency of our method have been introduced in Section 5.

Compared to other approaches, the strength of our method is to find directly the minimal number of bits needed at each control point to get a certain accuracy on the results. Consequently, it is not dependant of a certain number of data types (e.g. the IEEE754 formats) and its complexity does not increase as the number of data types increases. The information provided may also be used to generate computations in the fixpoint arithmetic with an accuracy guaranty on the results. Concerning scalability, we generate a linear number of constraints and variables in the size of the analyzed program. The only limitation is the size of the problem accepted by the solver. Note that the number of variables could be reduced by assigning the same precision to a whole piece of code (for example an arithmetic expression, a line of code, a function, etc.) Code synthesis for the fixpoint arithmetic and assigning the same precision to pieces of code are perspectives we aim at explore at short term.

At longer term, other developments of the present work are planned. First we wish to adapt the techniques developed in this article to the special case of Deep Neural Networks for which it is important to save memory usage and computational resources. Second, we aim at using our precision tuning method to guide lossy compression techniques for floating-point datasets [9]. In this case, the bit-level accuracy inferred by our method would determine the compression rate of the lossy technique.

References

[1] ANSI/IEEE: IEEE Standard for Binary Floating-point Arithmetic, std 754-2008 edn. (2008)
[2] Bailey, D.: Resolving numerical anomalies in scientific computation (03 2008)
[3] Ben Khalifa, D., Martel, M.: Precision tuning and internet of things. In: International Conference on Internet of Things, Embedded Systems and Communications, IINTEC 2019. pp. 80–85. IEEE (2019)
[4] Ben Khalifa, D., Martel, M.: Precision tuning of an accelerometer-based pedometer algorithm for iot devices. In: International Conference on Internet of Things and Intelligence System, IOTAIS 2020. pp. 113–119. IEEE (2020)
[5] Ben Khalifa, D., Martel, M., Adjé, A.: POP: A tuning assistant for mixed-precision floating-point computations. In: Formal Techniques for Safety-Critical Systems - 7th International Workshop, FTSCS 2019. Communications in Computer and Information Science, vol. 1165, pp. 77–94. Springer (2019)
[6] Chiang, W., Baranowski, M., Briggs, I., Solovyev, A., Gopalakrishnan, G., Rakamaric, Z.: Rigorous floating-point mixed-precision tuning. In: Proceedings of the 44th ACM SIGPLAN Symposium on Principles of Programming Languages, POPL. pp. 300–315. ACM (2017)
[7] Costan, A., Gaubert, S., Goubault, E., Martel, M., Putot, S.: A policy iteration algorithm for computing fixed points in static analysis of programs. In: Computer Aided Verification, 17th International Conference, CAV 2005, Edinburgh, Scotland, UK, July 6-10, 2005, Proceedings. Lecture Notes in Computer Science, vol. 3576, pp. 462–475. Springer (2005)
[8] Darulova, E., Horn, E., Sharma, S.: Sound mixed-precision optimization with rewriting. In: Proceedings of the 9th ACM/IEEE International Conference on Cyber-Physical Systems, ICCPS. pp. 208–219. IEEE Computer Society / ACM (2018)
[9] Diffenderfer, J., Fox, A., Hittinger, J.A.F., Sanders, G., Lindstrom, P.G.: Error analysis of ZFP compression for floating-point data. SIAM J. Sci. Comput. 41(3), A1867–A1898 (2019)
[10] Guo, H., Rubio-González, C.: Exploiting community structure for floating-point precision tuning. In: Proceedings of the 27th ACM SIGSOFT International Symposium on Software Testing and Analysis, ISSTA 2018. pp. 333–343. ACM (2018)
[11] Kotipalli, P.V., Singh, R., Wood, P., Laguna, I., Bagchi, S.: AMPT-GA: automatic mixed precision floating point tuning for GPU applications. In: Proceedings of the ACM International Conference on Supercomputing, ICS. pp. 160–170. ACM (2019)
[12] Lam, M.O., Hollingsworth, J.K., de Supinski, B.R., LeGendre, M.P.: Automatically adapting programs for mixed-precision floating-point computation. In: International Conference on Supercomputing, ICS’13. pp. 369–378. ACM (2013)
[13] McKeeman, W.M.: Algorithm 145: Adaptive numerical integration by simpson’s rule. Commun. ACM 5(12), 604 (1962)
[14] Morris, D., Saponas, T., Guillory, A., Kelner, I.: Recofit: Using a wearable sensor to find, recognize, and count repetitive exercises. Conference on Human Factors in Computing Systems - Proceedings (2014)
[15] de Moura, L.M., Bjørner, N.: Z3: an efficient SMT solver. In: Tools and Algorithms for the Construction and Analysis of Systems. LNCS, vol. 4963, pp. 337–340. Springer (2008)
[16] Papadimitriou, C.H.: On the complexity of integer programming. Journal of the ACM (JACM) 28(4), 765–768 (1981)
[17] Parker, D.S.: Monte carlo arithmetic: exploiting randomness in floating-point arithmetic. Tech. Rep. CSD-970002, University of California (Los Angeles) (1997)
[18] Rubio-González, C., Nguyen, C., Nguyen, H.D., Demmel, J., Kahan, W., Sen, K., Bailey, D.H., Iancu, C., Hough, D.: Precimonious: tuning assistant for floating-point precision. In: International Conference for High Performance Computing, Networking, Storage and Analysis, SC’13. pp. 27:1–27:12. ACM (2013)
[19] Schrijver, A.: Theory of linear and integer programming. John Wiley & Sons (1998)

7 Appendix

We need a lemma on some algebraic operations stable on the set of functions written as the min-max of a finite family of affine functions. The functions are defined on $\mathbb{R}^{d}$ .

Lemma 1

The following statements hold:

•

The sum of two min-max of a finite family of affine functions is a min-max of a finite family of affine functions.
•

The maximum of two min-max of a finite family of affine functions is a min-max of a finite family of affine functions.

Proof

Let $g$ and $h$ be two min-max of a finite family of affine functions and $f=g+h$ . We have $g=\min_{i}\max_{j}g^{ij}$ and $h=\min_{k}\max_{l}h^{kl}$ . Let $x\in\mathbb{R}^{d}$ . There exist $i,k$ such that $f(x)\geq\max_{j}g^{ij}(x)+\max_{l}h^{kl}(x)=\max_{j,l}g^{ij}(x)+h^{kl}(x)$ . We have also, for all $i,k$ , $f(x)\leq\max_{j}g^{ij}(x)+\max_{l}h^{kl}(x)=\max_{j,l}g^{ij}(x)+h^{kl}(x)$ . We conclude that $f(x)=\min_{i,k}\max_{j,l}g^{ij}(x)+h^{kl}(x)$ for all $x$ . We use the same argument for the max. $\Box$

Proposition 4

The following results hold:

1.

Let $\xi$ the constant function equal to 1. The system $S_{\xi}^{1}$ can be rewritten as $\{\mathsf{nsb}\in\mathbb{N}^{Lab}\mid F(\mathsf{nsb})\leq(\mathsf{nsb})\}$ where $F$ maps $\mathbb{R}^{Lab}\times\mathbb{R}^{Lab}$ to itself, $F(\mathbb{N}^{Lab}\times\mathbb{N}^{Lab})\subseteq(\mathbb{N}^{Lab}\times\mathbb{N}^{Lab})$ and has coordinates which are the maximum of a finite family of affine order-preserving functions.
2.

Let $\xi$ the function such that $\xi(\ell)$ equals the function defined at Fig. 5. The system $S_{\xi}$ can be rewritten as $\{(\mathsf{nsb},\mathsf{nsb_{e}})\in\mathbb{N}^{Lab}\times\mathbb{N}^{Lab}\mid F(\mathsf{nsb},\mathsf{nsb_{e}})\leq(\mathsf{nsb},\mathsf{nsb_{e}})\}$ where $F$ maps $\mathbb{R}^{Lab}\times\mathbb{R}^{Lab}$ to itself, $F(\mathbb{N}^{Lab}\times\mathbb{N}^{Lab})\subseteq(\mathbb{N}^{Lab}\times\mathbb{N}^{Lab})$ and all its coordinates are the min-max of a finite family of affine functions.

Proof

We only give details about the system $S_{\xi}^{1}$ (Figure 4). By induction on the rules. We write $L=\{\ell\in Lab\mid F_{\ell}\text{ is constructed }\}$ . This set is used in the proof to construct $F$ inductively.

For the rule (CONST), there is nothing to do. For the rule (ID), if the label $\ell^{\prime}=\rho(x)\in L$ then we define $F_{\ell^{\prime}}(\mathsf{nsb})=\max(F_{\ell^{\prime}}(\mathsf{nsb}),\mathsf{nsb}(\ell))$ . Otherwise $F_{\ell^{\prime}}(\mathsf{nsb})=\mathsf{nsb}(\ell)$ . As $\mathsf{nsb}\mapsto\mathsf{nsb}(\ell)$ is order-preseving and the maximum of one affine function, $F_{\ell^{\prime}}$ is the maximum of a finite family of order-preserving affine functions since $\max$ preserves order-preservation.

For the rules (ADD), (SUB), (MULT), (DIV), (MATH) and (ASSIGN), by induction, it suffices to focus on the new set of inequalities. If $\ell_{1}\in L$ , we define $F_{\ell_{1}}$ as the max with old definition and $RHS(\mathsf{nsb})$ i.e. $F_{\ell_{1}}(\mathsf{nsb})=\max(RHS(\mathsf{nsb}),F_{\ell_{1}}(\mathsf{nsb}))$ where $RHS(\mathsf{nsb})$ is the right-hand side part of the new inequality. If $\ell_{1}\notin L$ , we define $F_{\ell_{1}}(\mathsf{nsb})=RHS(\mathsf{nsb})$ . In the latter rules, $RHS(\mathsf{nsb})$ are order-preserving affine functions. It follows that $F_{\ell_{1}}$ is the maximum of a finite family of order-preserving affine functions.

The result follows by induction for the rule (SEQ).

The rules (COND) and (WHILE) are treated as the rules (ADD), (SUB), (MULT), (DIV), (MATH) and (ASSIGN), by induction and the consideration of the new set of inequalities.

The last rule (REQ) constructs $F_{\rho(x)}$ either as the constant function equal to $p$ at label $\rho(x)$ or the maximum of the old definition of $F_{\rho(x)}$ and $p$ if $\rho(x)\in L$ .

The proof for the system $S_{\xi}$ uses the same arguments and Lemma 1. $\Box$