Automatic Modular Abstractions for Template Numerical Constraints

We propose a method for automatically generating abstract transformers for static analysis by abstract interpretation. The method focuses on linear constraints on programs operating on rational, real or floating-point variables and containing linear assignments and tests. Given the specification of an abstract domain, and a program block, our method automatically outputs an implementation of the corresponding abstract transformer. It is thus a form of program transformation. In addition to loop-free code, the same method also applies for obtaining least fixed points as functions of the precondition, which permits the analysis of loops and recursive functions. The motivation of our work is data-flow synchronous programming languages, used for building control-command embedded systems, but it also applies to imperative and functional programming. Our algorithms are based on quantifier elimination and symbolic manipulation techniques over linear arithmetic formulas. We also give less general results for nonlinear constraints and nonlinear program constructs.


Introduction
Program analysis consists in deriving properties of the possible executions of a program from an algorithmic processing of its source or object code.Example of interesting properties include: "the program always terminates"; "the program never executes a division by zero"; "the program always outputs a well-formed XML document"; "variable x always lies between 1 and 3".There has been a considerable amount of work done since the late 1970s on sound methods of program analysis -that is, methods that produce results that are guaranteed to hold for all program executions, as opposed to bug finding methods such as program testing, which cannot provide such guarantees in general.
Static analysis by abstract interpretation is one of the various approaches to sound program analysis.Grossly speaking, abstract interpretation casts the problem of obtaining supersets of the set of reachable states of programs into a problem of finding fixed points of certain monotone operators on certain ordered sets, known as abstract domains.When dealing with programs operating on arithmetic values (integer or real numbers, or, more realistically, bounded integers and floating-point values), these sets are often defined by numerical constraints, and ordered by inclusion.One may, for instance, attempt to compute, for each program point and each variable, an interval that is guaranteed to contain all possible values of that variable at that point.The problem is of course how to compute these fixed points.Obviously, the smaller the intervals, the better, so we would like to compute them as small as possible.Ideally, we would like to compute the least fixed point, that is, the least inductive invariant that can be expressed using intervals.
The purpose of this article is to expose how to compute such least fixed points exactly, at least for certain classes of programs and certain abstract domains.Specifically, we consider programs operating over real variables using only linear comparisons (e.g.x + 2y ≤ 3 but not x 2 ≤ y), and abstract domains defined using a finite number of linear constraints i a i v i ≤ C, where the a i are fixed coefficients, C is a parameter (whose computation is the goal of the analysis) and the v i are the program variables.Such domains evidently include the intervals, where the constraints are of the form v i ≤ C and −v i ≤ C.
Not only can we compute such least fixed points exactly if all parameters are known, but we can also deal with the case where some of the parameters are unknowns, in which case we obtain the parameters of the least fixed point as explicit, algorithmic, functions of the unknowns.We can thus generate, once and for all, the abstract transformers for blocks of code: that is, those that map the parameters in the precondition of the block of code to a suitable postcondition.For instance, in the case of interval analysis, we derive explicit functions mapping the bounds on variables at the entrance of a loop-free block to the tightest bounds at the outcome of the block, or to the least inductive interval invariant of a loop.This allows modular analysis: it is possible to analyse functions or other kind of program modules (such as nodes in synchronous programming) in isolation of the rest of the code.
A crucial difference with other methods, even on loop-free code, is that we derive the optimal -that is, the most precise -abstract transformer for the whole sequence.In contrast, most static analyses only have optimal transformers for individual instructions; they build transformers for whole sequences of code by composition of the transformers for the individual instructions.Even on very simple examples, the optimal transformer for a sequence is not the composition of the optimal transformers of the individual instructions.Furthermore, most other methods are forced to merge the information from different execution traces at the end of "if-then-else" and other control flow constructs in a way that loses precision.In contrast, our method distinguishes different paths in the control flow graph as long as they do not go through loop iteration.
Our approach is based on quantifier elimination, a technique operating on logical formulas that transforms a formula with quantifiers into an equivalent quantifier-free formula: for instance, ∀x (x ≥ y ⇒ x ≥ 1) is turned into the equivalent y ≥ 1. Essentially, we define the result that we would like to obtain using a quantified logical formula, and, using quantifier elimination and further formula manipulations, we obtain an algorithm that computes this result.Thus, one can see our method as automatically transforming a non executable specification of the result of the analysis into an algorithm computing that result.
In §2, we shall provide background about abstract interpretation and quantifier elimination.In §3, we shall provide the main results of the article, that is, how to derive optimal abstract transformers for template linear constraint domains.In §4, we shall investigate various extensions to that framework, still based on quantifier elimination in linear real arithmetic; e.g.how to deal with floating-point values.In §5, we shall explain how to move from the single loop case to more complex control flow.In §6, we shall describe our implementation of the main algorithm and some limited extensions.In §7, we shall investigate extensions of the same framework using quantifier elimination on other theories, namely Presburger arithmetic and the theory of real closed fields (nonlinear real arithmetic).In §8, we shall list some related work and compare our method to other relevant approaches.

Background
In this section, we shall recall a few definitions, notations and results on program analysis by abstract interpretation, then quantifier elimination.
2.1.Abstract interpretation.It is well-known that, in the general case, fully automatic program analysis is impossible for any nontrivial property. 1Thus, all analysis methods must have at least one of the following characteristics: • They may bound the memory size of the studied program, which then becomes a finite automaton, on which most properties are decidable.Explicit-state modelchecking works by enumerating all reachable states, which is tractable only while implicit state model-checking represents sets of states using data structures such as binary decision diagrams [27].• They may restrict the programming language used, making it not Turing-complete, so that properties become decidable.For instance, reachability in pushdown automata is decidable even though their memory size is unbounded [17].• They may restrict the class of properties expressed to properties of bounded executions; e.g., "within the first 10000 steps of execution, there is no division by zero", as in bounded model checking [13].• They may be unsound as proof methods: they may fail to detect that the desired property is violated.Typically, bug-finding and testing programs are in that category, because they may fail to detect a bug.Some such analysis techniques are not based on program semantics, but rather on finding patterns in the program syntax [38].• They may be incomplete as proof methods: they may fail to prove that a certain property holds, and report spurious violations.Methods based on abstraction fall in that category.Abstraction by over-approximation consists in replacing the original problem, undecidable or very difficult to decide, by a simpler "abstract" problem whose behaviours are guaranteed to include all behaviours of the original problem.An example of an abstraction is to erase from the program all constructs dealing with numerical and pointer types (or replacing them with nondeterministic choices, if their value is used within a test), keeping only Boolean types (Fig. 1).Obviously, the behaviours of the resulting program encompass all the behaviours of the original program, plus maybe some extra ones.
Further abstraction can be applied to this Boolean program: for instance, the "3-value logic" abstraction [91] which maps any input or output variable to an abstract parameter taking its value in a 3-element set: "is 0", "is 1", "can be 0 or 1"; 2 for practical purposes it may be easier to encode these values using a couple of Booleans, respectively meaning "can be 0" and "can be 1", thus the abstract values (1, 0), (0, 1) and (1,1).The abstract value (0, 0) obtained for any variable at a program point means that this program point is unreachable.Given a vector of input abstract parameters, one for each input variable of the program, the forward abstract transfer function gives a correct vector of output abstract parameters, one for each output variable of the program.Quite obviously, in the absence of loops, it is possible to obtain a suitable forward abstract transfer function by applying simple logical rules to the Boolean function defined by the Boolean program.An effective implementation of the forward abstract transfer function can thus be obtained by a transformation of the source program.For programs operating over numerical quantities, a common abstraction is intervals.[33,34] To each input x, one associates an interval [x min , x max ], to each output x ′ an interval [x ′ min , x ′ max ].How can one compute the (x ′ min , x ′ max ) x∈V bounds from the (x min , x max ) x∈V ?The most common method is interval arithmetic: to each elementary arithmetic operation, one attaches an abstract counterpart that gives bounds on the output of the operation given bounds on the inputs.For instance, if one knows [a min , a max ] and [b min , b max ], and c = a + b, then one computes [c min , c max ] as follows: c min = a min + b min and c max = a max + b max .If a program point can be reached by several paths (e.g. at the end of an if-then-else construct), then a suitable interval [x min , x max ] can be obtained by a join of all the intervals for x at the end of these paths: The abstract transfer function defined by interval arithmetic is always correct, but is not necessarily the most precise.For instance, on example in Fig. 2, the best abstract transfer function maps any input range to z min = z max = 0, since the output z is always zero, 2 For brevity, we identify "false" with 0 and "true" with 1.For each statement (on the left) we use a corresponding optimal transformer (on the right), but the composition of these transformers is not optimal.For the sake of simplicity, all variables are considered to be real numbers.
while the one obtained by applying interval arithmetic to all program statements yields, in general, larger intervals.The weakness of the interval domain on this example is evidently due to the fact that it does not keep track of relationships between variables (here, that x = y).Relational abstract domains such as the octagons [74,76,77] or convex polyhedra [4,36,58] address this issue.Yet, neither octagons nor polyhedra provide analyses that are guaranteed to give optimal results.Consider the following program: Listing 5: Undistinguished paths i n t x ; i f ( x > 0 ) x= 1 ; e l s e x= −1; i f ( x == 0 ) x= 2 ; Obviously, the second test can never be taken, since x can only be ±1; however an interval, octagon or polyhedral analysis will conclude, after joining the informations for both branches of the first test, that x lies in [−1, 1] and will not be able to exclude the case x = 0.The final invariant will therefore be In contrast, our method, applied to this example, will correctly conclude that the optimal output invariant for x is [−1, 1].In fact, our method yields the same result as considering the (potentially exponentially large) set of paths between the beginning and the end of the program, and for each path, computing the least output interval, then computing the join of all these intervals.
We have so far left out programs containing loops; when programs contain loops or recursive functions, a central problem of program analysis is to find inductive invariants.In the case of Boolean programs, given constraints on the input parameters, the set of reachable states can be computed exactly by model-checking algorithms; yet, these algorithms do not give a closed-form representation of the abstract transfer function mapping input parameters to output parameters for the 3-value abstraction.In the case of numerical abstractions such as the intervals, octagons or polyhedra, the most common way to find invariants is through the use of a widening operator [34,35].
Intuitively, widening operators observe the sets of reachable states after N and N + 1 loop iterations and extrapolate them to a "candidate invariant".For instance, the widening The standard widening on convex polyhedra [36,58], here demonstrated on polyhedra in dimension 2 (polygons).The widening operator observes the sets of reachable states P 0 and P 1 at two consecutive iterations, and keeps only the constraints (polyhedral faces, here polygon edges) that are stable across iterations.The resulting P ∞ polyhedron is then proposed as an invariant.
operator, observing a sequence of intervals [0, 1], [0, 2], [0, 3] may wish to try [0, +∞).See Fig. 3 for an example with the standard widening operator on convex polyhedra.Let u 0 be the set of initial states of a loop, and let → τ be transition relation for this loop (σ → τ σ ′ means that σ ′ is reachable in one loop step from σ).The set of reachable states at the head of the loop is the least fixed point of f : , which is obtained as the limit of the ascending sequence defined by u n+1 = f (u n ).By abstract interpretation, we replace this sequence by an abstract sequence u ♯ n defined by N is an abstraction of the least fixed point of f and thus of the least invariant of the transition relation τ containing u 0 .
When one uses a widening operator ▽, the function p is an abstraction of f .The design of the widening operator ensures convergence of u ♯ n in finite time.The exceptions to the use of widening operators are static analysis domains that satisfy the ascending condition, such as the domain of linear equality constraints [65] and that of linear congruence constraints [55]: with f ♯ = f ♯ p the sequence u ♯ n always becomes stable within finite time.Listing 6: Circular buffer i = 0 ; while ( t r u e ) { i f ( nondet ( ) Again, widening operators provide correct results, but these results can be grossly over-approximated.Much of the literature on applied analysis by abstract interpretation discusses workarounds that give precise results in certain cases: narrowing iterations [33,34], widening "up to" [57, §3.2], "delayed" or with "thresholds" [15], etc.
A simple example of a program with one single variable where narrowing iterations fail to improve precision is Listing 6.This program is a much simplified version of a piece of code maintaining a circular buffer inside a large reactive control loop, where some piece of data is inserted only at certain loop iterations.We only kept the instructions relevant to the array index i and abstracted away the choice of the loop iterations where data is inserted as nondeterministic nondet().
If we analyse this loop using intervals with the standard widening with a widening point at the head of the loop, we obtain the sequence [0, 0], [0, 1], [0, 2] and then widening to [0, +∞).Narrowing iterations then fail to increase the precision.The reason is that the transition relation for this loop includes the identity function (when nondet() is false), thus the concrete function whose least fixed point defines the set of reachable states [35] satisfies X ⊆ φ(X) for all X (in other word, φ is expansive).Thus, once the widening operator overshoots the least fixed point, it can never recover it.
A similar problem is posed by: The usual widening iterations overshoot to [0, +∞), and narrowing does not recover from there.Furthermore, this example illustrates how widening makes analysis non monotonic: contrary to one could expect, having extra precision on the precondition of a program can result in worse precision for the inferred invariants.For instance, consider the above problem and replace the first line by assume(i>=0 && i<=9).Clearly, the resulting program is an abstraction of the above example, since it has strictly more behaviours (we allow 1, . . ., 9 as initial values for i).Yet, the analysis of the loop will yield a more precise behaviour: the interval [0, 9] is stable and the analysis terminates immediately.
Both these very simple examples can be precisely analysed using widening up to[57, §3.2], also known as widening with thresholds [15, Sec.7.1.2]:a preliminary phase collects all constants to which i is compared, and instead of widening to +∞, one widens to the next larger such constant if it exists -in this case, since x < 10 stands for x ≤ 9, widening with threshold would widen to 9, which is the correct value.This approach is not generalit fails if instead of the constant 10 we have some computed value.Of course, improvements are possible: for instance, one could analyse all the program up to this loop in order to prove that certain variables are constant, then use this information for setting thresholds for further loops.Yet, again, this is not a general approach.This is the second problem that this article addresses: how to obtain, in general, optimal invariants for certain classes of programs and numerical constraints.Furthermore, our methods provide these invariants as functions of the parameters of the precondition of the loop; this is one difference with our proposal, which computes the best, thus, again, they provide effective, optimal abstract transfer functions for loop constructs.
2.2.Quantifier elimination.Consider a set A of atomic formulas.The set U (A) of quantifier-free formulas is the set of formulas constructed from A using operators ∧, ∨ and ¬; the set Q(A) of quantified formulas is the set of formulas constructed from A using the above operators and the ∃ and ∀ quantifiers.Such formulas are thus trees whose leaves are the atomic formulas.A literal is an atomic formula or the negation thereof.The set of free variables F V (F ) of a formula F is defined as usual.A quantifier-free formula without variables is said to be ground.A formula without free variables is said to be closed; the existential closure of a formula F is F with existential quantifiers for all free variables prepended; the universal closure is the same with universal quantifiers.A quantifier-free formula is said to be in disjunctive normal form (DNF) if it is a disjunction of conjunctions, that is, is of the form (l ), and is said to be in conjunctive normal form (CNF) if it is a conjunction of disjunctions.Any quantifierfree formula can be converted into CNF or DNF by application of the distributivity laws though better algorithms exist, such as ALL-SAT modulo theory [81].
2.2.1.Linear real inequalities.Let A be the set of linear inequalities with integer or rational coefficients over a set of variables V .By elementary calculus, such inequalities can be equivalently written in the following forms: l(v 1 , v 2 , . . . ) ≥ C or l(v 1 , v 2 , . . . ) > C, with l a linear expression with integer coefficients over V and C a constant.Let us first consider the theory of linear real arithmetic (LRA): models of a formula F are mappings from F to the real field R, and notions of equivalence and satisfiability follow.Note that satisfiability and equivalence are not affected by taking models to be mappings from F to the rational field Q. Deciding whether a LRA formula is satisfiable is, again, a NP-complete problem known as satisfiability modulo theory (SMT) of real linear arithmetic.Again, practical implementations, known as SMT-solvers , are capable of dealing with rather large formulas; examples include Yices [45] and Z3 [41]. 3  The theory of linear real arithmetic admits quantifier elimination.For instance, the quantified formula ∀x (x ≥ y =⇒ x ≥ 3) is equivalent to the quantifier-free formula y ≥ 3. Quantified linear real arithmetic formulas are thus decidable: by quantifier elimination, one can convert the existential closure of the formula to an equivalent ground formula, the truth of which is trivially decidable by evaluation.The decision problem for quantified formulas over rational linear inequalities requires at least exponential time, thus quantifier elimination is at least exponential.[19, §7.4].[48,Th. 3] Weispfenning [107] discusses complexity issues in more detail.
Again, most quantifier elimination algorithms proceed by induction over the structure of the formula, and thus begin by eliminating the innermost quantifiers, progressively replacing branches of the formula containing quantifiers by quantifier-free equivalent branches.By application of the equivalence ∀x F ≡ ¬∃x ¬F , one can reduce the problem to eliminating existential quantifiers only.Consider now the problem of eliminating the existential quantifiers from ∃ x 1 . . .x n F where F is quantifier-free.We can first convert into DNF: where the C i are conjunctions, then to the equivalent formula (∃ x 1 . . .
We thus have reduced the quantifier elimination problem for general formula to the problem of quantifier elimination for conjunctions of linear inequalities.Remark that, geometrically, this quantifier elimination amounts to computing the projection of a convex polyhedron along the dimensions associated with the variables x 1 , . . ., x n , with the original polyhedron and its projection being defined by systems of linear inequalities.The Fourier-Motzkin elimination procedure [66, §5.4] converts ∃x 1 . . .x n C into an equivalent conjunction.This is what we refer to the "conversion to DNF followed by projection approach".This approach is good for quickly proving that the theory admits quantifier elimination, but it is very inefficient.We shall now see better methods.
Ferrante and Rackoff [47] proposed a substitution method [19, §7.3.1][87,§4.2][107]: a formula of the form ∃x F where F is quantifier-free is replaced by an equivalent disjunction where the e i are affine linear expressions built from the free variables of ∃x F .Note the similarity to the naive elimination procedure we described for Boolean variables: even though the existential quantifier ranges over an infinite set of values, it is in fact only necessary to test the formula F at a finite number of points (see Fig. 4).The drawback of this algorithm is that n is proportional to the square of the number of occurrences of x in the formula; thus, the size of the formula can be cubed for each quantifier eliminated.Loos and Weispfenning [69] proposed a virtual substitution algorithm4 [87, §4.4] that works along the same general ideas but for which n is proportional to the number of occurrences of x in the formula.Our benchmarks show that Loos and Weispfenning's algorithm is generally much more efficient than Ferrante and Rackoff's, despite the latter method being better known.[80,81] The main drawback of substitution algorithms is that the size of the formulas generally grows very fast as successive variables are eliminated.There are few opportunities for simplifications, save for replacing inequalities equivalent to false (e.g.0 < 0) by false and similarly for true, then applying trivial rewrites (e.g.false ∨ x x, false ∧ x false).Our experience is that these algorithms tend to terminate with out-of-memory [80,81].Scholl et al. [100] have proposed using and-inverter graphs (AIGs) and SAT modulo theory (SMT) solving in order to simplify formulas and keeping their size manageable.
Another class of algorithms improve on the conversion to DNF then projection approach, by combining both phases: we proposed an eager algorithm based on this idea [81], x y Figure 4: The gray zone S is the set of (x, y) solutions of formula F , whose atoms are the linear inequalities corresponding to the lines ∆ drawn with dashes.For fixed y = y 0 , the set of x such that F (x, y) is true is made of intervals whose ends lie within the set I of intersections of the y = y 0 line with the lines in ∆, drawn with a small circle.y = y 0 therefore has an intersection with S if and only if a point in I, or an interval with both ends in I ∪ {−∞, +∞}, lies within S. This condition can be tested using x → ±∞ and all midpoints to intervals with both ends in I, as per Ferrante and Rackoff [47], or, in addition to I ∪ {−∞}, for any element of I a point infinitesimally close to the right of it, as per Loos and Weispfenning [69].Both methods exploit the fact that the coordinates of all points from I (intersection of y = y 0 and a line from ∆) can be expressed as affine linear functions of y 0 .then a lazy version, which computes parts of formulas as needed [80].Instead of syntactic conversion to DNF, we use a SMT solver to point to successive elements of the DNF, and instead of using Fourier-Motzkin elimination, which tends to needlessly blow up the size of the formulas, we use libraries for computations over convex polyhedra, which can compute a minimal constraint representation of the projection of a polyhedron given by constraints (that is, inequalities) -which is the geometrical counterpart of performing elimination of a block of existentially quantified variables.Experiments have shown that such methods are generally competitive with substitution approaches, with different classes of benchmarks showing a preference for one of the two families of methods [80]; for the kinds of problems we consider in this article, it seems that the SMT+projection methods are more efficient [81].
The practical complexity of Presburger arithmetic is high.In particular, the formulas produced tend to be very complex, even when there exists a considerably simpler and "understandable" equivalent formula, as seen with experiments in §7.1.
Cooper and Pugh's procedures are very geometrical in nature.Integers, however, can also be seen as words over the {0, 1} alphabet, and sets of integers can thus be recognised by finite automata [88, §8].Addition is encoded as a 3-track automaton recognising that the number on the third track truly is the sum of the numbers on the first two tracks; equivalently, this encodes subtraction.Existential quantifier elimination just removes some of the tracks, making transitions depending on bits read on that track nondeterministic.Negation is complementation (which can be costly, thus explaining the high cost of quantifier alternation).Multiplication by powers of two is also easily encoded, and multiplications by arbitrary constants can be encoded by a combination of additions and multiplications by powers of two. 5 The same idea can be extended to real numbers written as their binary expansion, using automata on infinite words.
This leads to an interesting arithmetic theory, with two sorts of variables: reals (or rationals) and integers.This could be used to model computer programs, with integers for integer variables and reals for floating-point variables (if necessary by using the semantic transformations described in § 4.5).Boigelot et al. [16] described a restricted class of ωautomata sufficient for quantifier elimination.Becker et al. [11] implemented the Lira tool based on such ideas.Unfortunately, this approach suffers from two major drawbacks: the practical performances are very bad for purely real problems [81], and it is impossible to recover an arithmetical formula from almost all these automata.We therefore did not pursue this direction.

Nonlinear real arithmetic.
What happens if we do not limit ourselves to linear arithmetic, but also allow polynomials?Over the integers, the resulting theory is known as Peano arithmetic.It is well known that there can exist no decision procedure for quantified Peano arithmetic formulas. 6Since a quantifier elimination algorithm would turn a quantified formula without free variables into an equivalent ground formula, and ground formulas are trivially decidable, it follows that there can exist no quantifier elimination algorithm for this theory.
The situation is wholly different over the real numbers.The satisfiability or equivalence of polynomial formulas does not change whether the models are taken over the real numbers, the real algebraic numbers, or, for the matter, any real closed field; this theory is thus known as the theory of real closed fields, or elementary algebra.Tarski [105,106] and Seidenberg [101] showed that this theory admits quantifier elimination, but their algorithms had impractically high complexity.Collins [28] introduced a better algorithm based on cylindrical algebraic decomposition.For instance, quantifier elimination on ∃x ax 2 +bx+c = 0 by cylindrical algebraic decomposition yields Note the cylindrical decomposition: first, there is a case disjunction according to the values of c, then, for each disjunct for c, a case disjunction for the value of b; more generally, cylindrical algebraic decomposition builds a tree of case disjunctions over a sequence of variables v 1 , v 2 , . . ., where the guard expressions defining the cases for v i can only refer to v 1 , . . ., v i .This decomposition only depends on the polynomials inside the formula and not on its Boolean structure, and computing it may be very costly even if the final result is simple.This is the intuition why despite various improvements [10,24] the practical complexity of quantifier elimination algorithms for the theory of real closed fields remain high.The theoretical space complexity is doubly exponential [21,40].This is why our results in §7.2 are of a theoretical rather than practical interest.Minimal extensions to this formula language may lead to undecidability.This is for instance the case when one adds trigonometric functions: it is possible to define π as the least positive zero of the sine, then define the set of integers as the numbers k such that sin(kπ) = 0, and thus one can encode Peano arithmetic formulas into that language [2].Also, naive restrictions of the language do not lead to lower complexity.For instance, limiting the degree of the polynomials to two does not make the problem simpler, since formulas with polynomials of arbitrary degrees can be encoded as formulas with polynomials of degree at most two, simply by introducing new variables standing for subterms of the original polynomials.For instance, ∃x ax 3 +bx 2 +cx+d = 0 can be encoded, using Horner's form for the polynomial, as ∃x∃y∃z z = ax + b ∧ y = zx + c ∧ yx + d = 0. Certain stronger restrictions may however work; for instance, if the variables to be eliminated occur only linearly, then one can adapt the substitution methods described in §2.2.1.

Optimal Abstraction over Template Linear Constraint Domains
When applying abstract interpretation over domains of linear constraints, such as octagons [74,76,77], one generally applies a widening operator, which may lead to imprecisions.In some cases, acceleration techniques leading to precise results can be applied [52,53]: instead of attempting to extrapolate a sequence of iterates to its limit, as does widening, the exact limit is computed.In this section, we describe a class of constraint domains and programs for which abstract transfer functions of loop-free codes and of loops can be exactly computed; thus the optimality.Furthermore, the analysis outputs these functions in closed form (as explicit expressions combining linear expressions and functional if-then-else constructs), so the result of the analysis of a program fragment can be stored away for future use; thus the modularity.Our algorithms are based on quantifier elimination over the theory of real linear arithmetic ( §2.2).

3.1.
Template Linear Constraint Domains.Let F be a formula over linear inequalities.We call F a domain definition formula if the free variables of F split into n parameters p 1 , . . ., p n and m state variables s 1 , . . ., s m .We note γ F : As an example, the interval abstract domain for 3 program variables In this section, we focus on the case where linear inequalities whose left-hand side is fixed and the righthand sides are parameters.Such conjunctions define the class of template linear constraint domains [99].Particular examples of abstract domains in this class are: • the intervals (for any variable s, consider the linear forms s and −s); because of the inconvenience of talking intervals of the form [−a, b], we shall often taking them of the form [x min , x max ], with the optimal value for x min being a greatest lower bound instead of a least upper bound; • the difference bound matrices (for any variables s 1 and s 2 , consider the linear form s 1 − s 2 ); • the octagon abstract domain (for any variables s 1 and s 2 , distinct or not, consider the linear forms ±s 1 ± s 2 ) [74] • the octahedra (for any tuple of variables s 1 , . . ., s n , consider the linear forms ±s 1 • • •± s n ).[25] Remark that γ F is in general not injective, and thus one should distinguish the syntax of the values of the abstract domain (the vector of parameters p) and their semantics γ F ( p).As an example, if one takes F to be s and (1, 1, 3) define the same set for state variables s 1 and s 2 .If u ≤ v coordinate-wise, then γ F ( u) ⊆ γ F ( v), but the converse is not true due to the non-uniqueness of the syntactic form.
Take any nonempty set of states W ⊆ Q m .Take for all i = 1, . . ., m: We have therefore defined a map α F : P(R m ) → {⊥} ∪ (R ∪ {+∞}) n , and (α F , γ F ) form a Galois connection: α F maps any set to its best upper-approximation. 7 The fixed points of α F • γ F are the normal forms; the normal form of x ♯ is the minimal abstract element that stands for the same 7 See e.g.[35] for more information on Galois connection and their use in static analysis.Not all abstract interpretation techniques can be expressed within Galois connections.Indeed, there are abstract domains where there is not necessarily a best abstraction of a set of concrete states, e.g. the domain of convex polyhedra, which has no best abstraction for a disc, for instance.In this article, all abstractions, except the non-convex ones of § 4.3, are through Galois connections.
concrete set as x ♯ .8For instance, 3.2.Optimal Abstract Transformers for Program Semantics.We shall consider the input-output relationships of programs with rational or real variables.We first narrow the problem to programs without loops and consider programs built from linear arithmetic assignments, linear tests, and sequential composition.Noting a, b, . . . the values of program variables a, b . . .at the beginning of execution and a ′ , b ′ , . . . the output values, the semantics of a program P is defined as a formula P such that (a, b, . . ., a ′ , b ′ , . . . ) |= P if and only if the memory state (a ′ , b ′ , . . . ) can be reached at the end of an execution starting in memory state (a, b, . . .): where K is a real (in practice, rational) constant and L is a linear form, and b, c, d . . .are all the variables except a; where f 1 is P 1 F where a ′ has been replaced by a ′′ , b ′ by b ′′ etc., f 2 is P 2 F where a has been replaced by a ′′ , b by b ′′ etc.In addition to linear inequalities and conjunctions, such formulas contain disjunctions (due to tests and multiple branches) and existential quantifiers (due to sequential composition).
Note that so far, we have represented the concrete denotational semantics exactly.This representation of the transition relation using existentially quantified formulas is evidently as expressive as a representation by a disjunction of convex polyhedra (the latter can be obtained from the former by quantifier elimination and conversion to disjunctive normal form), but is more compact in general.
Consider now a domain definition formula . . ) ≤ p n on the program inputs, with parameters p and free variables s, and another on the program outputs, with parameters p ′ and free variables s ′ .Sound forward program analysis consists in deriving a safe post-condition from a precondition: starting from any state verifying the precondition, one should end up in the post-condition.Using our notations, the soundness condition is written The free variables of this relation are p and p ′ : the formula links the value of the parameters of the input constraints to admissible values of the parameters for the output constraints.
Note that this soundness condition can be written as a universally quantified formula, with no quantifier alternation.Alternatively, it can be written as a conjunction of correctness conditions for each output constraint parameter: Let us take a simple example: if P is the program instruction . Remark that this soundness condition is equivalent to a formula without quantifiers p ′ 1 ≥ p 1 + p 2 , which may be obtained through quantifier elimination.Remark also that while any value for p ′ 1 fulfilling this condition is sound (for instance, p ′ 1 = 1000 for p 1 = p 2 = 1), only one value is optimal (p ′ 1 = 2 for p 1 = p 2 = 1).An optimal value for the output parameter p ′ i is defined by . Again, quantifier elimination can be applied; on our simple example, it yields p ′ 1 = p 1 + p 2 .If there are n input constraint parameters p 1 , . . ., p n , then the optimal value for each output constraint parameter p ′ i is defined by a formula O ′ i with n+1 free variables p 1 , . . ., p n , p ′ i : Lemma 1.The formula O ′ i defined at Eq. 3.3, using the correctness subformula C ′ i from Eq. 3.2, defines p ′ i as the least possible value for the parameter of the constraint L i after executing the transition p from a state verifying constraints F .

Proof. O ′
i explicitly defines the least possible value of C ′ i .C ′ i explicitly defines all the acceptable values for parameter p ′ i in the postcondition constraint.This formula defines a partial function from Q n to Q, in the mathematical sense: for each choice of p 1 , . . ., p n , there exist at most a single p ′ i .The values of p 1 , . . ., p n for which there exists a corresponding p ′ i make up the domain of validity of the abstract transfer function.Indeed, this function is in general not defined everywhere; consider for instance the program: i f ( x >= 1 0 ) { y = n o n d e t e r m i n i s t i c c h o i c e i n a l l r e a l s ; } e l s e { y = 0 ; } and the function is defined only for p 1 < 10, with constant value 0.
At this point, we have a characterisation of the optimal abstract transformer corresponding to a program fragment P and the input and output domain definition formulas as n formulas (where n is the number of output parameters) O ′ i each defining a function (in the mathematical sense) mapping the input parameters p to the output parameter p ′ i .Another example: the absolute value function y := |x|, again with the interval abstract domain.The semantics of the operation is (x ≥ 0 ∧ y = x) ∨ (x < 0 ∧ y = −x); the precondition is x ∈ [x min , x max ] and the post-condition is y ∈ [y min , y max ].Acceptable values for (y min , y max ) are characterised by formula The optimal value for y max is defined by Quantifier elimination over this last formula gives as characterisation for the least, optimal, value for y max : We shall see in the next sub-section that such a formula can be automatically compiled into code (Listing 7).

3.3.
Generation of the Implementation of the Abstract Domain.Consider formula 3.5, defining an abstract transfer function.On this disjunctive normal form we see that the function we have defined is piecewise linear : several regions of the range of the input parameters are distinguished (here, x min + x max < 0 and x min + x max ≥ 0), and on each of these regions, the output parameter y max is a linear function of the input parameters.Given a disjunct (such as y max = −x min ∧ x min + x max < 0), the domain of validity of the disjunct can be obtained by existential quantifier elimination over the result variable (here ∃y max (y max = −x min ∧ x min + x max < 0)).The union of the domains of validity of the disjuncts is the domain of validity of the full formula.The domains of validity of distinct disjuncts can overlap, but in this case, since O ′ i defines a partial function in the mathematical sense, that is, a relation R(x, y) such that for any x there is at most one y such that R(x, y), the functions defined by such disjuncts coincide on their overlapping domains of validity.
This suggests a first algorithm for conversion to an executable form: For each disjunct C j , obtain the validity domain V j as a conjunction of linear inequalities; this corresponds to a projection of the polyhedron defined by C j onto variables p 1 , . . ., p n , parallel to p ′ i .Then, solve C j for p ′ i (obtain p ′ i as a linear function v i of the p 1 , . . ., p n ).
(3) Output the result as a cascade of if-then-else and assignments.Consider now the example at the end of §3.2 and especially Formula 3.5, defining y max as a function of x min and x max : ( Its validity domain is ∃y max C 1 , that is, x min + x max ≥ 0. Furthermore, on this validity domain, we can solve for y max as a function of x min and x max , and we obtain y max = x max .We therefore can print out the following test: . It can similarly be turned into another test, which we put in the "else" branch of the preceding one.We thus obtain: Observe that in the above program, the second test is redundant.If we had been more clever, we would have realized that in the "else" branch of the first test, x min + x max < 0 always holds, and thus it is useless to test this condition.Furthermore, in an if-then-else cascade obtained with the above method, the same condition may have to be tested several times.We shall now propose an algorithm that constructs an if-then-else tree such that no useless tests are generated. Algorithm 1 ToITEtree(F, p ′ , T ): turn a formula defining p ′ as a function of the other free variables of F into a tree of if-then-else constructs, assuming that T holds.
The idea of the algorithm is as follows: • Each path in the if-then-else tree corresponds to a conjunction C of conditions (if one goes through the "if" branch of if (a) and the "else" branch of if (b), then the path corresponds to a ∧ ¬b).• The formula O ′ i is simplified relatively to C, a process that prunes out conditions that are always or never satisfied when C holds.
• If the path is deep enough, then the simplified formula becomes a conjunction.This conjunction, as a relation, is a partial function from the p 1 , . . ., p n to the p ′ i we wish to compute.Again, we solve this conjunction to obtain p ′ i as an explicit function of the p 1 , . . ., p n .For instance, in the above example, we obtain y max as a function of x min and x max .We say that two formulas A and B are equivalent, denoted by A ≡ B, if they have the same models, and we say that they are equivalent modulo a third formula T , denoted by A ≡ T B, if A ∧ T ≡ B ∧ T .Intuitively, this means that we do not care what happens where F is false, thus the terminology "don't care set" sometimes found for ¬F .
QElimDNFModulo( v, F, T ) is a function that, given a possibly empty vector of variables v, a formula F and a formula T , outputs a quantifier-free formula F ′ in disjunctive normal form such that F ′ ≡ T ∃ v F and no useless predicates are used.In [81], we have presented such a function as a variant of quantifier elimination.
We need some more auxiliary functions.Predicates(F ) returns the set of atomic predicates of F .Solve(C, p ′ ) solves a conjunction of inequalities C for variable p ′ .If C does not contain redundant inequalities, then it is sufficient to look for inequalities of the form p ′ ≥ L or p ′ ≤ L and output p ′ = L. 9 Finally Choose(P ) chooses any predicate in the set of predicates P ; one good heuristic seems to be to choose the most frequent atomic predicate where p ′ i does not occur.Let us take, as a simple example, Formula 3.5.We wish to obtain y max as a function of x min and x max , so in the algorithm ToITEtree we set p ′ △ = y max .C 1 is the first disjunct We project C 1 and C 2 parallel to y max , obtaining respectively P 1 = (x min + x max ≥ 0) and P 2 = (x min + x max < 0).We choose K to be the predicate x min + x max ≥ 0 (in this case, the choice does not matter, since P 1 and P 2 are the negation of each other).
• The first recursive call to ToITEtree is made with QElimDNFModulo(y max , F, T ) will then simply output the formula "true".It then suffices to solve for y max in y max = x max .This yields the formula for computing the correct value of y max in the cases where x min + x max ≥ 0.
• The second recursive call is made with T △ = (x min + x max < 0. The result is y max = −x min , the formula for computing the correct value of y max in the cases where x min + x max < 0. These two results are then reassembled into a single if-then-else statement, yielding the program at the end of §3.2.
The algorithm terminates because paths of depth d in the tree of recursive calls correspond to truth assignments to d atomic predicates among those found in the domains of validity of the elements of the disjunctive normal form of F .Since there is only a finite number of such predicates, d cannot exceed that number.A single predicate cannot be assigned truth values twice along the same path because the simplification process in QElimDNFModulo erases this predicate from the formula.This algorithm seems somewhat unnecessarily complex.It is possible that techniques based upon AIGs, performing simplification with respect to "don't care sets" [100], could also be used.
3.4.Least Inductive Invariants.We have so far considered programs without loops.The analysis of program loops, as well as proofs of correctness of programs with loops using Floyd-Hoare semantics, is based upon the notion of inductive invariants.A set of states I is deemed to be an inductive invariant for a loop if it contains the initial state and it is stable by the loop iteration -in other words, if it is impossible to move from a state inside I to a state outside I by one iteration of the loop.The intersection of all inductive invariants is also an inductive invariant -in fact, it is the least inductive invariant.
Any property true over the least inductive invariant of a loop is true throughout the execution of that loop; for this reason, some authors call such properties invariants (without the "inductive" qualifier), while some other authors call invariants what we call inductive invariants.
It would be interesting to be able to compute the least invariant (inductive or noninductive) within our chosen abstract domain; in other words, compute the least element in our abstract domain that contains the least inductive invariant of the loop or program.Unfortunately, this is in general impossible; indeed, doing so would entail solving the halting problem.Just take any program P , create a fresh variable, and consider the program that initialises x to 0, runs P , and then sets x to 1. Clearly, the least invariant interval for this program and variable x is [0, 0] if P does not terminate, and [0, 1] if it does terminate.
We thus settle for a simpler problem: find the least inductive invariant within our abstract domain, that is, the least element in our abstract domain that is an inductive invariant.Note that even on very simple examples, this can be vastly different from computing the least invariant within the abstract domain.Take for instance the simple program of Fig. 5, which has a couple of real variables (x, y) and rotates them by 45 • at each iteration.The (x, y) couple always stays within the square [−1, 1] × [−1, 1], so this square is an invariant within the interval domain.Yet this square is not an inductive invariant, for it is not stable if rotated by 45 • ; in fact, the only inductive invariant within the interval domain is (−∞, +∞) × (−∞, +∞), which is rather uninteresting!Note that on the same figure, the octagon abstract domain would succeed (and produce a regular octagon centered on (0, 0)).  ) can be reached from the state (s 1 , . . ., s m ) in exactly one iteration of the loop.A set W ⊆ Q m is said to be an inductive invariant for the head of the loop if ∀ s ∈ W, ∀ s ′ ( s, s ′ ) |= G =⇒ s ′ ∈ W .We seek inductive invariants of the shape defined by F ′ , thus solutions for p ′ of the stability condition: (3.6)Not only do we want an inductive invariant, but we also want the initial states of the program to be included in it.The condition then becomes This is an invariant condition for the head A of a loop written as:   , the least invariant is an 8point star, and its best approximation inside the abstract domain of intervals is the square [−1, 1] 2 .However, this square is not an inductive invariant: no rectangle (product of intervals) is stable under the iterations, thus there is no abstract inductive invariant within the interval domain.Using the domain of convex polyhedra, one would obtain a regular octagon.
Alternatively, one can consider an invariant condition for location B. The condition then becomes This alternate condition is very similar to the previous one (a universally quantified formula with no alternation, since the ∃ is negated).For the sake of simplicity, we shall only discuss the treatment of formula H; formula H alt can be treated in the same way.This formula links the values of the input constraint parameters p 1 , . . ., p n to acceptable values of the invariant constraint parameters p ′ 1 , . . ., p ′ n ′ .In the same way that our soundness or correctness condition on abstract transformers allowed any sound post-condition, whether optimal or not, this formula allows any inductive invariant of the required shape as long as it contains the precondition, not just the least one.
The intersection of sets defined by p ′ 1 and p ′ 2 is defined by min( p ′ 1 , p ′ 2 ).More generally, the intersection of a family of sets, unbounded yet closed under intersection, defined by p ′ ∈ Z is defined by min{p ′ | p ′ ∈ Z}.We take for Z the set of acceptable parameters p ′ such that p ′ defines an inductive invariant and ∀ s, F =⇒ F ′ ; that is, we consider only inductive invariants that contain the set I = { s | F } of precondition states.We deduce that p ′ i is uniquely defined by: p ′ i = min(∃p ′ 1 , . . ., p ′ i−1 , p ′ i+1 , . . ., p ′ n ′ H) which can be rewritten as The free variables of this formula are p 1 , . . ., p n , p ′ i .This formula defines a function (in the mathematical sense) defining p ′ i from p 1 , . . ., p n .As before, this function can be compiled to an executable version using cascades or trees of tests.

Lemma 2. The formula O ′
i defined at Eq. 3.9 defines the least value of p ′ i for an inductive invariant of the shape defined by F for the transition relation defined by G.
Proof.Similarly as for lemma 1, the formula H defined at Eq. 3.7 defines a set Y of admissible p ′ 1 , . . ., p ′ n ′ such that Thus the overall operation of the analysis method: We start from quantified formulas O ′ i defining the least inductive invariant of the loop in the abstract domain (Lemma 2).We eliminate quantifiers from these formulas; since this does not change their models, the resulting formulas without quantifiers also define the least inductive invariant in the abstract domain.
• Either the problem had no precondition parameters p 1 , . . ., p n , and thus each formula O ′ i has only one variable p ′ i .It consists of linear (in)equalities, and has a single model for p ′ i , which is straightforward to extract.Collect these values, one obtains the values defining the least invariant (p ′ 1 , . . ., p ′ n ′ ) in the abstract domain.
• Either the problem has precondition parameters and one employs one of the methods in §3.3 to obtain equivalent executable code.The overall correctness of the method is quite tautological.We start from a nonexecutable specification of the least inductive invariant in the abstract domain in the form of quantified formulas.We eliminate the quantifiers from these formulas, then process them into equivalent executable code.At all steps, we have preserved logical equivalence with the original definition.In short, we have synthetized the implementation of the best transformer from its specification.

Simple Loop Example.
To show how the method operates in practice, let us consider first a very simple example (nondet() is a nondeterministic choice): Let us abstract i at the head of the loop using an interval [i min , i max ].For simplicity, we consider the case where the loop is at least entered once, and thus i = 0 belongs to the invariant.For better precision, we model each comparison x = y over the integers as x >= y + 1 ∨ x <= y − 1; similar transformations apply for other operators.The formula expressing that such an interval is an inductive invariant is: Quantifier elimination produces: The formulas defining optimal i min and i max are: We note that this invariant is only valid for n > 0, which is unsurprising given that we specifically looked for invariants containing the precondition i = 0.The output abstract transfer function is therefore: The case disjunction n < 2 looks unnecessary, but is a side effect of the use of rational numbers to model a problem over the integers.The resulting abstract transfer function is optimal, but on such a simple case, one could have obtained the same using polyhedra [36] or octagons [74].
We have already noted ( §2.1) that even with we replace n by the constant 10, the classical widening/narrowing approach will fail to identify the least invariant of this loop, and that extra techniques such have widening with thresholds have to be used.This block has three input streams e1, e2, and e3, and one output stream s1.In intuitive terms, s1 is the same as e1 but where the maximal slope of the signal between two successive clock ticks is bounded by e2, thus the name rate limiter.At some points in time ; } e l s e { s1max = e3max ; } (modelled by nondeterministic choice), the value of the signal is reset to that of the third input e3.
We are interested in the input-output behaviour of that block: obtain bounds on the output s1 of the system as functions of bounds on the inputs (e1, e2, e3).One difficulty is that the s1 output is memorised, so as to be used as an input to the next computation step.The semantics of such a block is therefore expressed as a fixed point.
We wish to know the least inductive invariant of the form s 1min ≤ s 1 ≤ s 1max under the assumption that e 1min ≤ e 1max ∧ e 2min ≤ e 2max ∧ e 3min ≤ e 3max .The stability condition yields, after quantifier elimination and projection on s 1max of the condition s 1max ≥ e 1max ∧ s 1max ≥ e 3max .Minimisation then yields an expression that can be compiled to an if-thenelse tree (Listing 9).This result, automatically obtained, coincides with the intuition that a rate limiter (at least, one implemented with exact arithmetic) should not change the range of the signal that it processes.
This program fragment has a rather more complex behaviour if all variables and operations are IEEE-754 floating-point, since rounding errors introduce slight differences of regimes between ranges of inputs.Rounding errors in the program to be analysed introduce difficulties for analyses using widenings, since invariant candidates are likely to be "almost stable", but not truly stable, because of these errors.Again, there exist workarounds so that widening-based approaches can still operate [15,Sec. 7.1.4].We shall see in §4.5 how to correctly abstract floating-point behaviour within our framework; unfortunately, the formulas produced tend to be very large due to many case disjunctions.The implementation of the abstract transformer produced for the above rate limiter in floating-point does not fit within one page of article, this is why we omitted it.

Extensions of the framework using real linear arithmetic
We shall describe here a few extensions to the class of programs and domains that we can handle, all of which are still based on quantifier elimination over real linear arithmetic.(In §7, we shall investigate extensions using other arithmetic theories.)4.1.Emptiness.We have so far supposed that the statement where the inductive invariant is computed is reachable, and thus that there exists some nonempty set of initial states that constrain the inductive invariant from below.More generally, and especially for the constructs described in §4.4 and §5.1, it may be necessary to encode the bottom element ⊥ of the abstract domain, which represents the empty set of states.This can be done using one Boolean variable per element that might be empty: instead of template parameters (p 1 , . . ., p n ), we will have (b, p 1 , . . ., p n ), with the semantics that γ(false, p 1 , . . ., p n ) = ∅ and γ(true, p 1 , . . ., p n ) = { s | ∀i L i ( s) ≤ p i }.
Sankaranarayanan et al. [99] use p i = −∞ for this, but in our framework, infinities themselves have to be encoded using Booleans, as we'll see in the next subsection.Furthermore, if we have an abstract element in normal form with a constraint L i (v 1 , . . . ) ≤ −∞, it means that all p j are −∞ and thus it is sufficient to have a single Boolean variable for all of them.We would like to obtain as a result of the analysis that x lies in [0, ∞).Yet, what will happen is that the formula describing the couples (x min , x max ) defining inductive invariants x min ≤ x ≤ x max will be simplified to false.
More annoyingly, with the following program: x = y ; while ( t r u e ) { x = x + 1 ; } we would at least hope to infer that x min = y min .Yet, if we look for an invariant of the form x min ≤ x ≤ x max , there is no solution for any value of (y min , y max )!In §7.1 we shall see an example where it is actually interesting to infer the domain of values of the precondition where the least inductive invariant interval is finite, and that this domain can simply be obtained by existential quantification on the parameters of the inductive invariant followed by elimination.But in general, this is not what we want; instead, we would prefer to allow infinite values in the p and p ′ .
This can be easily achieved by a minor alteration to our definitions.Each parameter p i is replaced by two parameters p b i and p ∞ i .p ∞ i is constrained to be in {0, 1} (if the quantifier elimination procedure in use allows Boolean variables, then p ∞ i can be taken as a Boolean variable); p ∞ i = 0 means that p i is finite and equal to After this rewriting, all formulas are formulas of the theory of linear inequalities without infinities and are amenable to the appropriate algorithms.
Unfortunately, the added combinatorial complexity induced by these Boolean variables tends to lead to intolerable computation times in the quantifier elimination procedures.Further work is needed, probably in the direction of better quantifier elimination procedures for combinations of Boolean and real quantified variables.Alternatively, one can envision directly including reasoning about infinities inside these procedures, though this is of course a delicate matter because of the possibility of generation of indeterminate forms ∞ − ∞ if formulas are handled without special care.

4.3.
Non-Convex Domains.Section 3.1 constrains formulas to be conjunctions of inequalities of the form L i ≤ p i .What happens if we consider formulas that may contain disjunctions?
The template linear constraint domains of section 3.1 have a very important property: they are closed under (infinite) intersection; that is, if we have a family p ∈ W , then there exists p 0 such that p∈W γ F ( p) = γ F ( p 0 ) (besides, p 0 = inf{ p | p ∈ W }).This is what enables us to request the least element that contains the exact post-condition, or the least inductive invariant in the domain: we take the intersection of all acceptable elements.
Yet, if we allow non-convex domains, there does not necessarily exist a least element γ F ( p) such that S ⊆ γ F ( p).Consider for instance S = {0, 1, 2} and F representing unions of two intervals (( We consider formulas F built out of linear inequalities L i (s 1 , . . ., s n ) ≤ p i as atoms, conjunctions, and disjunctions.By induction on the structure of F , we can show that γ F : (R ∪ {−∞}) n → P(R n ) is inf-continuous; that is, for any descending chain ( p i ) i∈I such that lim i p i = p ∞ , then γ F ( p i ) is decreasing and i∈I γ F ( p i ) = γ F ( p ∞ ).The property is trivial for atomic formulas, and is conserved by greatest lower bounds (∧) as well as binary least upper bounds (∨).
Let us consider a set S ⊆ P(R n ), stable under arbitrary intersection (or at least, greatest lower bounds of descending chains).S can be for instance the set of invariants of a relation, or the set of over-approximations of a set W . γ −1 F (S) is the set of suitable domain parameters; for instance, it is the set of parameters representing inductive invariants of the shape specified by F , or the set of representable over-approximations of W . γ −1 F (S) is stable under greatest lower bounds of descending chains: take a descending chain ( p i ) i∈I , For instance, consider p = (a, b, c, d), so as not to obtain the same solutions twice (by exchange of (a, b) and (c, d)) or solutions with empty intervals.By quantifier elimination over G, we obtain (a = 0 4.4.Domain Partitioning.Non-convex domains, in general, are not stable under intersections and thus "best abstraction" problems admit multiple solutions as minimal elements of the set of correct abstractions.As explained in the preceding subsection, is not very satisfactory nor efficient for analysis.There are, however, non-convex abstract domains that are stable under intersection and thus admit least elements as well as the template linear constraint domains of Sec.3.1: those defined by partitioning of the state space. If, for instance, we know that a program or system behaves differently according to the sign of x, then we can decide in advance to partition the state space into x ≥ 0 and x < 0. On each element of the partition, we can have a separate template (example Fig. 6).We can accommodate this within our framework.
More formally, consider pairwise disjoint subsets (C i ) i∈I of the state space Q m , and abstract domains stable under intersection (S i ) i∈I , S i ⊆ P(C i ).Elements of the partitioned abstract domain are unions i∈I s i where s i ∈ S i .If ( i s i,j ]) j∈J is a family of elements of the domain, then j∈J i∈I s i,j = i∈I j∈J s i,j ; that is, intersections are taken separately in each C i . in other words, the parameters of the templates on each element of the partition can be dealt with independently of each other.
Note the difference with the general disjunctive domains of §4.3: in the general disjunctive domains, there is no partition fixed a priori, this is why we may have several incomparable minimal elements [0, 0] ∪ [1, 2] and [0, 1] ∪ [2,2] in the domain of disjunctions of two intervals representing the same set {0, 1, 2}.For a fixed partition C i and corresponding domains S i , for any set X, for every i, there is a least element representing X ∩ C i in domain S i .This motivates the following construct: Take a family (F i [ p]) i∈I of formulas defining template linear constraint domains (conjunctions of linear inequalities L i (s 1 , . . ., s n ) ≤ p i ) and a family (C i ) i∈I of formulas such that for all i and i ′ , All the techniques of §3.1 then apply.
For instance, by choosing , we can obtain Figure 6.The above constructions are equivalent to assigning a separate control point to each element in the partition, with guards leading to these points according to the C i , and then performing as described in §5.1.

4.5.
Floating-Point Computations.Real-life programs do not operate on real numbers; they operate on fixed-point or floating-point numbers.Floating point operations have few of the good algebraic properties of real operations; yet, they constitute approximations of these real operations, and the rounding error introduced can be bounded.
In IEEE floating-point [61], each atomic operation (noting ⊕, ⊖, ⊗, ⊘, √ f for operations so as to distinguish them from the operations +, −, ×, /, √ over the reals) is mathematically defined as the image of the exact operation over the reals by a rounding function. 10This rounding function, depending on user choice, maps each real x to the nearest floating-point value r n (x) (round to nearest mode, with some resolution mechanism for non representable values exactly in the middle of two floating-point values), r −∞ (x) the greatest floating-point value less or equal to x (round toward −∞), r +∞ (x) the least floating-point value greater or equal to x (round toward +∞), r 0 (x) the floating-point value of the same sign as x but whose magnitude is the greatest floating-point value less or equal to |x| (round toward 0).If x is too large to be representable, r(x) = ±∞ depending on the size of x 10 We leave aside the peculiarities of some implementations, such as those of most C compilers over the 32-bit Intel platform where there are "extended precisions" types used for some temporary variables and expressions can undergo double rounding [84].
The semantics of the rounding operation cannot be exactly represented inside the theory of linear inequalities. 11As a consequence, we are forced to use an axiomatic overapproximation of that semantics: a formula linking a real number x to its rounded version r(x).
Miné [75] uses an inequality |r(x) − x| ≤ ε rel • |x| + ε abs , where ε rel is a relative error and ε abs is an absolute error, leaving aside the problem of overflows.The relative error is due to rounding at the last binary digit of the significand, while the absolute error is due to the fact that the range of exponents is finite and thus that there exists a least positive floating-point number and some nonzero values get rounded to zero instead of incurring a relative error (or get rounded to a denormal, see below).
Because our language for axioms is richer than the interval linear forms used by Miné, we can express more precise properties of floating-point rounding.We recall briefly the characteristics of IEEE-754 floating-point numbers.Nonzero floating point numbers are represented as follows: x = ±2 e m where 1 ≤ m < 2 is the mantissa or significand, which has a fixed number p of bits, and e (E min ≤ e ≤ E max is the exponent).The difference introduced by changing the last binary digit of the mantissa is ±s.ε last where ε last = 2 −(p−1) : the unit in the last place or ulp.Such a decomposition is unique for a given number if we impose that the leftmost digit of the mantissa is 1 -this is called a normalised representation.Except in the case of numbers of very small magnitude, IEEE-754 always works with normalised representations.There exists a least positive normalised number m normal and a least positive denormalized number m denormal , and the denormals are the multiples of m denormal less than m normal .All representable numbers are multiples of m denormal .
We shall now attempt to define an imprecise axiomatisation of the relationship between the result of a floating-point operation and the result of the corresponding ideal operation over real numbers.For this, we shall distinguish operations plus and minus on the one hand, multiplication and division on the other hand, since the former have properties of which we can take advantage.
Let us consider addition or subtraction x = ±a ± b. Suppose that 0 ≤ x ≤ m normal .a and b are multiples of m denormal and thus a − b is exactly represented as a denormalized number; therefore r(x) = x.If x > m normal , then |r(x) − x| ≤ ε rel .x.In other words, if the result of an addition or subtraction is denormal, then it is exact. 12The cases for x ≤ 0 are symmetrical.
We therefore obtain the following axiomatisation of the rounding of a positive real number x, result of a floating-point addition or subtraction, to a floating-point value r, using round-to-nearest: 11 To be pedantic, since IEEE floating-point formats are of a finite size, the rounding operation could be exactly represented by enumeration of all possible cases; this would anyway be impossible in practice due to the enormous size of such an enumeration.
12 William Kahan has written extensively about the advantages of the existence and proper handling of denormals, otherwise known as gradual underflow, as opposed to flushing to zero all numbers too small to be approximated by normal numbers, a process known as flush to zero.For instance, with gradual underflow, a ⊖ b = 0 is equivalent to a = b, quite a desirable property.See e.g.[63, p. 6].Unfortunately, certain architectures do not implement gradual underflow, for the sake of efficient; then one has to use the Round and not the Round ± predicate.
We then use this predicate to construct an axiomatisation for rounding of numbers of any sign: 2) Equivalently, this last formula can be simplified into the equivalent: Consider now multiplication or division x = a ⊗ b or x = a ⊘ b.Here, we cannot assume that if the result is denormal, then it is exact.A correct axiomatization of rounding for positive numbers is: where ε abs = m denormal /2.Now for rounding for any sign: To each floating-point expression e, we associate a "rounded-off" variable r e , the value of which we constrain using Round ± (r e , e) or Round(r e , e).For instance, a expression e = a⊕ b is replaced by a variable r e , and the constraint Round ± (r e , a + b) is added to the semantics.In the case of a compound expression e = ab + c, we introduce e 1 = ab, and we obtain Round ± (r e , r e 1 + c) ∧ Round(r e 1 , ab).If we know that the compiler uses a fused multiplyadd operator, we can use Round(r e , ab + c) instead.
The drawbacks of such axiomatization of floating-point operations are that they introduce case disjunctions into formulas, leading to extra work for quantifier elimination procedures, and also that they make very unnatural coefficients appear -that is, rational numbers with large numerator and denominator, such as 1 + ε rel or 1 − ε rel , or very large integers such as 1/m normal .To go back to our very simple rate limiter running example ( §3.4.3), analysis times are multiplied by 12 if one stops assuming that floating-point variables behave like reals, and instead uses some axiomatization [81, p. 9].Furthermore, the formula produced is very large and hardly readable, doing many case disjunctions.Perhaps it is possible to simplify such large formulas by allowing some limited overapproximation -for instance, we can replace the function mapping p to p ′ as defined by (0 ≤ p ≤ 1 ≤ p ′ = p)∨(p > 1 ≤ p ′ = (1+ǫ)p) by the function defined by 0 ≤ p∧p ′ = (1+ǫ)p, since the latter formula always gives a solution for p ′ greater or equal to that of the former.Even better, such simplifications could be performed during the elimination procedure.Further investigations are needed in that respect.4.6.Integers.We have mentioned in §2.2.2 that Presburger arithmetic admits quantifier elimination.We therefore could apply quantifier elimination in Presburger arithmetic, similarly as we do with respect to real linear arithmetic.Unfortunately, as we shall see in §7.1, such an approach suffers from explosion in the size of the formulas and the cost of the algorithms.
Instead, we used a relaxation approach: all integers are treated as reals; strict inequalities a < b where both sides are integers are recoded a ≤ b − 1.For instance, if the program contains a if-then-else over x ≤ y, then x ≤ y is a precondition for the "then" branch, and x ≥ y + 1 is a precondition of the "else" branch.Note that this means that traces of execution such that y < x < y + 1 are considered to fail.
In some cases, such as the McCarthy 91 function example from §5.2, it is necessary to constraint the reasoning procedures so that they consider that the negation of a ≤ b is a ≥ b + 1.We hope that improvements of quantifier elimination algorithms will be able to allow a more elegant approach.
Another issue is that in many programming languages, integers are bounded and that arithmetic operations are actually performed modulo 2 n (with n typically 8, 16, 32 or 64).The problem then lies within an enormous, but finite, state space.Clever techniques for reasoning about bit-vector arithmetic are being investigated by the SMT-solving community.Again, we hope that future work will provide good quantifier elimination techniques for this arithmetic, or combinations thereof with the linear theory of reals.4.7.Nonlinear constructs.In §7.2 we explain how to fully deal with nonlinear program constructs and templates, at the expense of very high computational complexity.
A more practical approach is to linearise the program expressions [78][76, ch. 6].If we encounter an assignment z := x * y and we know, perhaps through rough interval analysis, an interval y ∈ [y min , y max ], we can use the following abstraction of the semantics of this assignment: If we know intervals for both x and y, then we can apply Eq. 4.6 to both (x, [y min , y max ]) and (y, [x min , x max ]), and take the conjunction of the resulting formulas.

Complex control flow
We have so far assumed no procedure call, and at most one single loop.We shall see here how to deal with arbitrary control flow graphs and call graph structures.5.1.Arbitrary control graph and loop nests.In Sec.3.4, we have explained how to abstract a single fixed point.The method can be applied to multiple nested fixed points by replacing the inner fixed point by its abstraction.For instance, assume the rate limiter of Sec.3.4.3 is placed inside a larger loop.One may replace it by its abstraction: / * s o m e t h i n g * / } } / * and s i m i l a r f o r s1min * / Alternatively, we can extend our framework to an arbitrary control flow graph with nested loops, the semantics of which is expressed as a single fixed point.We may use the same method as proposed by Gulwani et al. [56, §2] and other authors.First, a cut set of program locations is identified; any cycle in the control flow graph must go through at least one program point in the cut set.In widening-based fixed point approximations, one classically applies widening at each point in the cut set.A simple method for choosing a cut set is to include all targets of back edges in a depth-first traversal of the control-flow graph, starting from the start node; in the case of structured program, this amounts to choosing the head node of each loop.This is not necessarily the best choice with respect to precision, though [56, §2.3]; Bourdoncle [18,Sec. 3.6] discusses methods for choosing such as cut-set.
To each point in the cut set we associate an element in the abstract domain, parameterised by a number of variables.The values of these variables for all points in the cut-set define an invariant candidate.Since paths between elements of the cut sets cannot contain a cycle, their denotational semantics can be expressed simply by an existentially quantified formula.Possible paths between each source and destination elements in the cut-set defined a stability condition (Formula 3.6).The conjunction of all these stability conditions defines acceptable inductive invariants.As above, the least inductive invariant is obtained by writing a minimisation formula (Sec.3.4).
Let us consider the loop nest in Listing 10.We choose program points A and B as cutset.At program point A, we look for an invariant of the form I A (i, j) △ = i min,A ≤ i ≤ i max,A , and at program point B, for an invariant of the form I B (i, j) The (somewhat edited for brevity) stability formula is written: where If[b, e 1 , e 2 ] = (b ∧ e 1 ) ∨ (¬b ∧ e 2 ).
Replacing I A and I B into this formula, then applying quantifier elimination, we obtain a formula defining all acceptable tuples (i min,A , i max,A , i min,B , i max,B , j min , j max , δ min , δ max ).Optimal values are then obtained by further quantifier elimination: i min,A = i min,B = j min = 0, i max,A = i max,B = 19, j max = 20, δ min = 1, δ max = 19.
The same example can be solved by replacing 20 by another variable n as in Sec.3.4.2.

5.2.
Procedures and Recursive Procedures.We have so far considered abstractions of program blocks with respect to sets of program states.A program block is considered as a transformer from a state of input program states to the corresponding set of output program states.The analysis outputs a sound and optimal (in a certain way) abstract transformer, mapping an abstract set of input states to an abstract set of output states.Assuming there are no recursive procedures, procedure calls can be easily dealt with.We can simply inline the procedure at the point of call, as done in e.g.Astrée [14,15,37].Because inlining the concrete procedure may lead to code blowup, we may also inline its abstraction, considered as a nondeterministic program.Consider a complex procedure P with input variable x and output variable x.We abstract the procedure automatically with respect to the interval domain for the postcondition (z min ≤ z ≤ z max ); suppose we obtain z max := 1000; z min := x then we can replace the function call by z <= 1000 && z >= x.This is a form of modular interprocedural analysis: considering the call graph, we can abstract the leaf procedures, then those calling the leaf procedures and so on.We can also do likewise for nested loops: abstract the innermost loop, then the next innermost one, etc.This method is however insufficient for dealing with recursive procedures.
In order to analyse recursive procedures, we need to abstract not sets of states, but sets of pairs of states, expressing the input-output relationships of procedures.In the case of recursive procedures, these relationships are the least solution of a system of equations.
To take a concrete example, let us consider McCarthy's famous "91 function" [70,71], which, non-obviously, returns 91 for all inputs less than 101: The concrete semantics of that function is a relationship R between its input n and its output r.It is the least solution of We look for a inductive invariant of the form 4.3).By replacing R by I into inclusion 5.2, and by universal quantification over n, r, n 2 , we obtain the set of admissible parameters for invariants of this shape.By quantifier elimination, we obtain (C = 91)∧(δ = ∆ = −10)∧(A = 101)∧(B = 100) within a fraction of a second using Mjollnir (see Sec. 6).
In this case, there is a single acceptable inductive invariant of the suggested shape.In general, there may be parameters to optimise, as explained in Sec.3.4.The result of this analysis is therefore a map from parameters defining sets of states to parameters defining sets of pairs of states (the abstraction of a transition relation).This abstract transition relation (a subset of X × Y where X and Y are the input and output state sets) can be transformed into an abstract transformer in X ♯ → Y ♯ as explained in Sec.3.2.Such an interprocedural analysis may also be used to enhance the analysis of loops [72].

Implementations and Experiments
We have implemented the techniques of Sec. 3 in quantifier elimination packages, including Mathematica 13 and Reduce 3.8 14 + Redlog 15 in addition to our own package, Mjollnir [81]. 16We ignore which exact techniques are implemented within Mathematica. 17Redlog appears to implement some virtual substitution method [44,107].
As test cases, we took a library of operators for synchronous programming, having streams of floating-point values as input and outputs.These operators are written in a restricted subset of C and take as much as 20 lines.A front-end based on CIL [85] converts them into formulas, then these formulas are processed and the corresponding abstract transfer functions are pretty-printed.Since for our application, it is important to bound numerical quantities, we chose the interval domain.
For instance, the rate limiter presented in Sec.3.4.3was extracted from that library.Since this operator includes a memory (a variable whose value is retained from a call to the operator to the next one), its data-flow semantics is expressed using a fixed-point.When considered with real variables, the resulting expanded formula was approximately 13 Mathematica is a commercial computer algebra package available under an unfree license from Wolfram Research [108]. 14Reduce is a computer algebra package from Anthony C. Hearn, now available under a modified BSD licence. 15Redlog is an extension to Reduce for working over quantified formulas. 16Mjollnir is available under a free license from the author's home page.In addition to the author's own quantifier elimination techniques, it implements Ferrante and Rackoff and Loos and Weispfenning's. 17Loos and Weispfenning's quantifier elimination procedure is used by Mathematica to perform simplifications over linear inequalities [108, §A.9.5], but we are unsure whether this is the algorithm called by the Reduce function.
1000 characters long, and with floating point variables approximately 8000 characters long.Despite the length of these formulas, they can be processed by Mjollnir in a matter of seconds.The result can then be saved once and for all.Analysers such as Astrée [14,15,37] must have special knowledge about such operators, otherwise the analysis results are too coarse (for instance, the intervals do not get stabilized at all).The Astrée development team therefore had to provide specialized, hand-written analyses for certain operators.In contrast, all linear floating-point operators in the library were analysed within a fraction of a second using the method in the present article, assuming that floating-point values in the source code were real numbers.If one considered instead the abstraction of floating-point computations using real numbers from Sec. 4.5, computation times did not exceed 17 seconds per operator; the formulas produced are considerably more complex than in the real case.Note that this computation is done once and for all for each operator; a static analyser can therefore cache this information for further use and need not recompute abstractions for library functions or operators unless these functions are updated.
Our analyser front-end currently cannot deal with non-numerical operations and data structures (pointers, records, and arrays).It is therefore not yet capable of directly dealing with the real control programs that e.g.Astrée accepts, which do not consist purely of numerical operators.We plan to integrate our analysis method into a more generic analyser.Alternatively, we plan to adapt a front-end for synchronous programming languages such as Simulink, a tool widely used by control/command engineers.
The correctness of the methods described in this article does not rely on any particularity of the quantifier elimination procedure used, provided one also has symbolic computation procedures for e.g.putting formulas in disjunctive normal form and simplifying them.The difference between the various quantifier elimination and simplification procedures is efficiency; experiments showed that ours was vastly more efficient than the others tested for this kind of application.For instance, our implementation of our quantifier elimination algorithm [81] was able to complete the analysis of the rate limiter of Sec.3.4.3,implemented over the reals, in 1.4 s, and in 17 s with the same example over floating-point numbers, while Redlog took 182 s for the former and could not finish the latter, and Mathematica could analyse neither (out-of-memory).On other examples, our quantifier elimination procedure is faster than the other ones, or can complete eliminations that the others cannot.

Extensions to other numerical domains
We have so far concerned ourselves solely with real (and possibly Boolean) variables appearing in linear arithmetic formulas.In §4, we have seen how to reason over certain other data types (integers, floating-point values), but again modeling them as real numbers.Yet, in §2.2, we pointed out that linear real arithmetic is not the only arithmetic theory with quantifier elimination algorithms; other well-known examples include Presburger arithmetic (linear integer arithmetic) and the theory of real closed fields (that is, polynomial real arithmetic).
In this section, we report on using quantifier elimination in both these theories for computing optimal transformers or fixed points.Unfortunately, experiments have shown that the high complexity of the quantifier elimination algorithms for these theories and the lack of simplifications for the formulas they produce preclude their use in practice.The results in this section are thus mainly of theoretical interest.7.1.Presburger arithmetic.The approach from §4.6 relaxes integers to reals.It therefore yields correct results, but might lead to overapproximations that could be avoided if integers were used instead.What if we used quantifier elimination on Presburger arithmetic instead?Consider the following simple example: i n t a , m; . . .i = 0 ; while ( i < m) { i = i +a ; } Let us abstract the loop variable a using an interval [l, h].These bounds must satisfy I △ = l ≤ 0 ≤ h (initialisation of the loop) and The condition for the existence of a finite range of values for i is therefore ∃l∃u I ∧ S. 18Intuitively, this condition is equivalent to m < 0 ∨ a ≥ 0. Yet, the formula produced by the quantifier elimination for Presburger arithmetic from Redlog yields a very large formula, more than one page long, which we have therefore omitted. 19Redlog cannot even perform simplifications such as replacing i = 0 ∨ i = 1 ∨ i = 2 by 0 ≤ i ≤ 2.
In comparison, the real relaxation gives exact and fast results: Eliminating quantifiers from ∃l∃u I ∧ S R yields immediately the answer a ≥ 0 ∨ m ≤ −1.

7.2.
Real polynomial constraint domains.We now consider the abstraction of program states (in R V ) using domains defined by polynomial constraints, a natural extension of those seen previously ( §3.1); the orthogonal extensions from § 4 also apply.Instead of quantifier elimination in linear real arithmetic, we shall use quantifier elimination in the theory of real closed fields.One difference, though, is that we will not be able to produce nice, closed form formulas, at least not in general.
7.2.1.Method.We generalise the constructs of §3, except those of § 3.3, to formulas over polynomial inequalities.The same results hold: • For any loop-free program code, and any template polynomial abstract domain with parameters p 1 , . . ., p n , there is a family of formulas F 1 , . . ., F n that uniquely defines the optimal parameters p ′ 1 , . . ., p ′ n ′ of the postcondition with respect those p 1 , . . ., p n in the precondition (the free variables of F i are among p 1 , . . ., p n , p ′ i ).
• For any loop, and any template polynomial abstract domain with parameters p 1 , . . ., p n , there is a family of formulas F 1 , . . ., F n that uniquely defines the optimal parameters p ′ 1 , . . ., p ′ n ′ of the least inductive invariant for that loop, with respect those p 1 , . . ., p n in the precondition (the free variables of F i are among p 1 , . . ., p n , p ′ i ).The main obstacle is the high cost of quantifier elimination in the theory of real closed fields.The other crucial difference is that it is in general impossible to move from such a formula to a formula computing p ′ i from p 1 , . . ., p n , as we did in § 3.3.By performing the cylindrical algebraic decomposition with variables p 1 , . . ., p n first, we could obtain the tree structure with case disjunctions, as the output of Algorithm 1.But at the leaves, we would obtain formulas defining p ′ i as a specific root of a polynomial in the variable p ′ i , with coefficients themselves polynomials in p 1 , . . ., p n .
In Sec.3.3, we explained how to turn a formula over linear arithmetic defining a partial function from p 1 , . . ., p n to p ′ i -that is, a relation between (p 1 , . . ., p n , p ′ i ) such that for any choice of p 1 , . . ., p n there is at most one suitable p ′ i -into an algorithmic function, expressed using linear assignments and if-then-else over linear inequalities.Can we do the same here?That is, can we turn a formula over nonlinear arithmetic defining a partial function from p 1 , . . ., p n to p ′ i into a simple algorithm written using, say, if-then-else tests and normal arithmetic operators as well as the n-th root operations n √ ?For instance, if we have a formula p ′ ≥ 0 ∧ p ′ 2 = 2p, we would like to obtain p ′ = √ 2p.
Unfortunately, this is impossible in the general case.The Abel-Ruffini theorem, from Galois theory, states that for polynomials in one variable of degrees higher or equal to 5, there is in general no way to express the value of the roots using only arithmetic operations (+, −, ×, /) and radicals ( n √ ).Thus, we cannot hope to obtain in general a simple algorithm expressing p ′ i as a function of p 1 , . . ., p n using tests, arithmetic operations (+, −, ×, /) and radicals ( n √ ).Let us now assume that there are no precondition parameters p 1 , . . ., p n or, equivalently, that we know exactly their value.Would it be at least possible to compute the values of the p ′ i ?p ′ i is defined as the only solution of a logical formula with a single free variable, built using polynomial arithmetic.By putting this formula into DNF, we can reduce the problem to computing the only solution, if any, of a conjunction of polynomial inequalities and equalities.Such a solution can be computed to arbitrary precision ; in fact, for any ǫ, one can obtain bounds [l, h] such that l ≤ p ′ i ≤ h and h − l ≤ ǫ.Unfortunately, the cost of such computations is high.[10] 7.2.2.Experiments.Consider l x ≤ x ≤ h x , l y ≤ y ≤ h y and the problem of generating the optimal abstract transfer function for the multiplication operation, z := x * y.We wish to obtain l z and h z such that l z ≤ z ≤ h z , and l z and h z are optimal (l z is maximal, h z is minimal).We first define the set of admissible (not necessarily minimal) h z : Now we define the least value for h z : The free variables of this formula are l x , h x , l y , h y and h z .Mathematica 7 performs quantified elimination by cylindrical algebraic decomposition on this formula in 4.3 s and yields a large formula (Fig. 7) with many case disjunctions, ) .(7.5) most notably on the sign of l x , h x , l y , h y .This is quite natural: the monotonicity of the function y → xy changes according to the sign of x.This function is equivalent to the much more terse h z = max(l x l y , h x l y , l x h y , h x h y ) (7.6)This illustrates the limit of our approach on nonlinear problems: even on simple program constructions and with simple invariants, quantifier elimination takes nonneglible time and outputs complicated formulas.We therefore did not pursue this direction further.

Related work
Since the first numerical abstract analysis techniques were proposed in the 1970s, there has been considerable work on improving precision, efficiency, or both.Without attempting to be exhaustive, we shall now describe a few of the approaches and how they differ from ours.
8.1.Relational abstract domains and modular analysis.There is a sizeable amount of literature concerning relational numerical abstract domains; that is, domains that express constraints between numerical variables.Convex polyhedra were proposed in the 1970s [36,58], and there have been since then many improvements to the technique; a bibliography was gathered by Bagnara et al. [4].Algorithms on polyhedra are costly and thus a variety of domains intermediate between simple interval analysis and convex polyhedra were proposed [25,74,99].
It is possible to use relational abstract domains such as polyhedra to model the input/output relationship of a program, function or block [58, §7.2, p. 112].Instead of considering only the current values of the program variables (v 1 , . . ., v n ) at the various program points, one also considers the initial values (v 0 1 , . . ., v 0 n ) of these variables at the beginning of the program, function or block; thus the computed polyhedra, for instance, relate (v 0 1 , . . ., v 0 n , v 1 , . . ., v n ).We employed this approach when dealing with recursive procedures ( §5.2).Such an approach is modular: one can for instance analyse a procedure in such a way, and plug the result of the analysis at each point of call; though of course one loses optimality.It also provides for modular analysis of loop nests: one first analyses the innermost loop, and then replace this innermost loop by the result of the analysis, considered as a nondeterministic program; one then proceeds to the next innermost loop.
One limit of this approach is that the relationship between input and output is constrained by the abstract domain.Most numerical abstract domains concern convex relations: difference bound constraints, octagons, polyhedra etc. are all geometrically convex (given two points a, b in the concretisation, the segment [a, b] is also in the concretisation).Note that the result of the analysis of the absolute value function ( §3.2), as expressed by Rel.3.5, or that of the rate limiter ( §3.4.3), are piecewise linear but not convex.
The idea of producing procedure summaries [103] as formulas mapping input bounds to output bounds is not new.Rugina and Rinard [97], in the context of pointer analysis (with pointers considered as a base plus an integer offset), proposed a reduction to linear programming.This reduction step, while sound, introduces an imprecision that is difficult to measure in advance; our method, in contrast, is guaranteed to be "optimal" in a certain sense.Rugina and Rinard's method, however, allows some nonlinear constructs in the program to be analysed.Martin et al. [72] proposed applying interprocedural analysis to loops.
Seidl et al. [102] also produce procedure summaries as numerical constraints.Our procedure summaries are implementations of the corresponding abstract transformer over some abstract domain, while theirs outputs a relationship between input and output concrete values.Their analysis considers a convex set of concrete input-output relationships, expressed as simplices, a restricted class of convex polyhedra.This restriction trades precision for speed: the generator and constraint representations of simplices have approximately the same size, while in general polyhedra exponential blowup can occur.Tests by arbitrary linear constraints cannot be adequately represented within this framework.Seidl et al. [102,Sec. 4] propose deferring those constraints using auxiliary variables; this, however, loses some precision.Their analysis and ours are therefore incomparable, since they make different choices between precision and efficiency.
Lal et al. [67] proposed an interprocedural analysis of numerical properties of functions using weighted pushdown automata.The "weights" are taken in a finite height abstract domain, while the domains we consider have infinite height.8.2.Computations of exact fixed points.The limitations of the widening approach explained in §2.1 have been recognized for long.There has therefore been extensive research about computing precise inductive invariants, if possible the least inductive invariant inside the abstract domain considered.
Several methods have been proposed to synthesize invariants without using widening operators [29,32,98].In common with us, they express as constraints the conditions under which some parametric invariant shape truly is an invariant, then they use some resolution or simplification technique over those constraints.Again, these methods are designed for solving the problem for one given set of constraints on the inputs, as opposed to finding a relation between the output or fixed-point constraints and the input constraints.In some cases, the invariant may also not be minimal.
Bagnara et al. [5,6] proposed improvements over the "classical" widenings on linear constraint domains [58].Gopan and Reps [54] introduced "lookahead widenings": standard widening-based analysis is applied to a sequence of syntactic restrictions of the original program, which ultimately converges to the whole program; the idea is to distinguish phases or modes of operation in order to make the widening more precise.Gonnord [52], Gonnord and Halbwachs [53], Leroux and Sutre [68] have proposed acceleration techniques: when the transition relation τ is of certain particular forms, it is possible to compute its transitive closure τ + exactly or with small imprecision.Typically, acceleration techniques have difficulties dealing with programs where the control flow is not flat, for instance when there are paths through a loop body that affect the iteration variables in different ways, such as the circular buffer example (Listing 6).
Adjé et al. [1], Costan et al. [31], Gaubert et al. [49] proposed a "policy iteration" or "strategy iteration" approach, 20 by downwards iterations providing successive over-approximations of the least fixed point.Their approach can fail to converge to the least fixed point, for instance with expansive semantics such as those of the "circular buffer" example (Listing 6), though for some classes of programs it converges to the least fixed point [1].
Gawlitza and Seidl [51] proposed another policy iteration approach, which is guaranteed to provide the least fixed point of the system of abstract equations.In contrast to the above method, they use upwards iterations, so each value computed is an under-approximation of the abstract least fixed point.They extended that approach to template linear constraint domains [50].The differences with our approach are twofold: • Their approach computes the least fixed point of a system of min/max abstract equations, derived from the source code of the program.In intuitive terms, min's correspond to conditions (and closures operators in relational domains), and max's to "merge points" in the control flow graph (end of if-then-else).This approach thus incurs the same problem of "undistinguished paths" as the example from §2.1.
Even then, the policy iteration algorithm may iterate across a number of iterations exponential in the number of merge points.An alternative would be to consider a cut-set for the control flow graph and distinguish each path between two points in the cut-set. in other words, for a single loop, one would consider each individual control path inside the loop body.The number of such paths is exponential in the number of tests, and thus the 20 The terminology comes from game theory.In broad terms, their consider equations with max/min operators, which are similar to the "minimax" operators appearing in the definition of the value of games in game theory.Choosing which argument of a "min" is used corresponds, in game theory terms, of choosing a strategy or policy for a player that tries to minimise the value of the game.
"max" operation in the abstract equations would also have an exponential number of arguments.Such an approach would compute the same result as ours; however, it would be unworkable without further work to get rid of the explicitly exponential number of arguments in the "max" operation.
• Their approach needs all preconditions exactly known.In contrast, we compute it as an explicit function of the precondition.In short, our invariants are parametric in the precondition, while theirs are not.Gulwani et al. [56] have also proposed a method for generating linear invariants over integer variables, using a class of templates.The methods described in the present article can be applied to linear invariants over integer variables in two ways: either by abstracting them using rationals (as in examples in Sec.3.4.2,5.1), either by replacing quantifier elimination over rational linear arithmetic by quantifier elimination over linear integer arithmetic, also known as Presburger arithmetic ( §4.6).Gulwani et al. instead chose to first consider integer variables as rationals, so as to be able to compute over rational convex polyhedra, then bound variables and constraint parameters so as to model them as finite bit vectors, finally obtaining a problem amenable to SAT solving.Program variables are finite bit vectors in most industrial programming languages, and parameters to useful invariants over integer variables are often small, thus their approach seems justified.We do not see, however, how their method could be applied to programs operating over real or floating-point variables, which are the main motivation for the present article.8.3.Limitations of template-based approaches.Much, if not all, of the published results on computing least inductive invariants in abstract domains, or at least abstract fixed points, deal with template domains [49,50,56], including intervals [31,51].The fundamental reasons for this are: • The domain of convex polyhedra is not closed under infinite intersection.Thus, in general, there is no best abstraction of a set of states.In general, there is no least inductive invariant inside the domain.• These methods replace the problem of dealing with arbitrarily complex shapes such as convex polyhedra by dealing with a vector of real numbers in finite, fixed dimension.The coefficients of these vectors are then amenable to a variety of constraint solving techniques.A common criticism of template-based approaches, is that they suppose that one knows the interesting templates beforehands -in the case of linear constraint domains, interesting directions in space.In contrast, methods based on general convex polyhedra infer these directions themselves, through the convex hull and widening operations [36,58].For instance, an analysis using standard polyhedra of the following program will infer that 2x = y, while a template based approach would succeed in doing so only if a template of the form 2x − y has been provided: x = y = 0 ; while ( t r u e ) { x = x + 1 ; y = y + 2 ; } We understand and share that criticism.In some cases there exist "natural" templates, such as intervals, or when dealing with timing or scheduling constraints, difference bounds v i −v j ≤ C, but in general, finding the correct templates seems a hard problem.A suggestion is to run iterations using general convex polyhedra and look at the stable directions of the faces of these polyhedra.
We think however that criticism of template domains should be part of wider considerations on how to choose the proper abstract domain.Abstract interpretation essentially replaces the unsolvable problem of computing the least inductive invariant of the concrete problem by the solvable problem of computing an inductive invariant in an abstract domain -in some lucky cases such as the one dealt with in this article, actually obtaining the least inductive invariant in the abstract domain.The choice of the abstract domain is at present somewhat arbitrary, typically hardwired into the analysis tool, at best chosen by some command-line flags.
Convex polyhedra [36,58,59] are popular because many programming idioms naturally exhibit convexity -for instance, the set of loop indices (i, j, k, . . . ) of nested loops occurring in numerical analysis programs is often convex.Yet, one can easily think of programs where some interesting properties are not convex (a test |x| ≥ 1, for instance).There are some cases where it is important not to enforce convexity and instead implement disjunctive domains, capable of representing properties such as x ≤ −1 ∨ x ≥ 1; an example is trace partitioning [93].Thus, a static analysis tool based on general convex polyhedra also enforces an a priori convex shape that might not be representative of the useful program invariants.
While convex polyhedra are not enough, they are often "too much": their complexity is too high for many applications.Indeed, even for less costly constraint domains such as octagons [74,77], it is simply too expensive to compute constraints between all visible program variables, so some analysers choose a priori to consider relations only between certain variables -an approach knowing as packing in the Astrée tool [37,39].Again, the packing choice is a form of a priori template guess made by the analyser, using heuristics that look at the way the program is organized.
Further research is obviously needed in how to choose and adapt abstract domains so that they can represent interesting inductive properties at reasonable costs.This problem is similar to the problem of finding the correct predicates in predicate abstraction.For this, various methods for finding predicates in addition to those syntactically present in the program have been proposed, especially those based on the analysis of spurious counterexamples (counterexample-guided abstraction refinement or CEGAR). 21Perhaps similar ideas could be employed for suggesting suitable additional numerical templates, finer-grained packing or disjunctions.8.4.Relational domains beyond polyhedra.In earlier works, we have proposed a method for obtaining input-output relationships of digital linear filters with memories, taking into account the effects of floating-point computations [79].This method computes an exact relationship between bounds on the input and bounds on the output, without the need for an abstract domain for expressing the local invariant; as such, for this class of problems, it is more precise than the method from this article.This technique, however, cannot be easily generalized to cases where the operator block contains tests and other nonlinear constructs; the semantics of nonlinear constructs must be approximated by e.g.interval analysis.
There have been several published approaches to finding nonlinear relationships between program variables.One approach obtains polynomial equalities through computations on ideals using Gröbner bases [94,95].This work only deals with equalities (not inequalities), uses a classical approach of computing output constraints from a set of input constraints (instead of finding relationships between the two sets of constraints), and deals with loops using a widening operator.In comparison, our approach abstracts whole program fragments, and is modular -it is possible to "plug" the result of the analysis of a procedure at the location of a procedure call, though of course this is less precise than inlining the procedure.
Since nonlinear relations are notoriously costly to compute upon, Bagnara et al. [7] have proposed using further abstraction to be able to reduce the problem to computations over convex polyhedra.
Kapur [64] also proposed to use quantifier elimination to obtain invariants: he considers program invariants with parameters, and derives constraints over those parameters from the program.Our work improves on his by noting that least invariants of the chosen shape can be obtained, not just any invariant; that the abstraction can be done modularly and compositionally (a program fragment can be analysed, and the result of its analysis can be plugged into the analysis of a larger program), or combined into a "conventional" abstract interpretation framework (by using invariants of a shape compatible with that framework), and that the resulting invariants can be "projected" to obtain numerical quantities.

Conclusion and future prospects
Writing static analysers by hand has long been found tedious and error-prone.One may of course prove an existing analyser correct through assisted proof techniques, which removes the possibility of soundness mistakes, at the expense of much increased tediousness.In this article, we proposed instead effective methods to synthesize abstract domains by automatic techniques.The advantages are twofold: new domains can be created much more easily, since no programming is involved; a single procedure, testable on independent examples, needs be written and possibly formally proved correct.To our knowledge, this is the first effective proposal for generating numerical abstract domains automatically, and one of the few methods for generating numerical summaries.Also, it is also the only method so far for computing summaries of floating-point functions.
We have shown that floating-point computations could be safely abstracted using our method.The formulas produced are however fairly complex in this case, and we suspect that further over-approximation could dramatically reduce their size.There is also nowadays significant interest in automatizing, at least partially, the tedious proofs that computer arithmetic experts do and we think that the kind of methods described in this article could help in that respect.
We have so far experimented with small examples, because the original goal of this work was the automatic, on-the-fly, synthesis of abstract transfer functions for small sequences of code that could be more precise than the usual composition of abstract of individual instructions, and less tedious for the analysis designer than the method of pattern-matching the code for "known" operators with known mathematical properties.A further goal is the precise analysis of longer sequences, including integer and Boolean computations.We have shown in Sec.4.4 how it was possible to partition the state space and abstract each region of the state-space separately; but naive partitioning according to n Booleans leads to 2 n regions, which can be unbearably costly and is unneeded in most cases.We think that automatic refinement and partitioning techniques [62] could be developed in that respect.
The main practical application that we envision is to be able to analyse numerical operator blocks from synchronous programming languages such as Simulink, 22 Scicos, 23  Lustre, 24 Scade 25 or Sao, 26 which are widely used for programming control systems [3], particularly in the automative and avionic industries.In order to obtain good analysis precision, such blocks often have to be analysed as a whole instead of decomposing them into individual components and applying individual transfer functions, as in our rate limiter example.The static analysis tool Astrée [15,37,38,39,43,104] outputs few, if any, false alarms on some classes of control programs because it has specific specialized transfer functions for certain operator blocks or coding patterns.Such transfer functions had to be implemented by hand; the techniques described in the present article could have been used to implement some of them automatically and even on-the-fly.
There are two important drawbacks to our method, which make it currently only useful for very precise analysis of small parts of programs.The first is that we need to "see" the whole of the loop or function that we are analysing, the instructions of which must belong to the class of constructs that we are capable of dealing with, or at least can be abstracted by them.In contrast, iterative techniques are more tolerant: they see the program state locally, at each program point, and the numerical analysis may easily interact with other analyses, such as pointers [14,15].The second issue is the high cost of quantifier elimination.Despite our work on new algorithms [81], in which we are still making progress, scalability remains an issue.

Listing 1 :
Original program i n t x , y ; bool a ; i f ( x < y ) a = f a l s e ; Listing 2: Boolean program bool a ; i f ( nondet ( ) ) a = f a l s e ;

Figure 1 :
Figure 1: Transformation of a program into a Boolean program by erasing the numeric part and replacing tests over numerical quantities by nondeterministic choice (nondet() nondeterministically returns true or false).
[a, b] ⊓ [c, d] = [min(a, c), max(b, d)].Again, for a loop-free program, one can obtain a suitable effective forward abstract transfer function by a program transformation of the source code.

Listing 3 :Figure 2 :
Figure2: The transformer for the interval domain obtained by composition of locally optimal abstract transformers is imprecise.For each statement (on the left) we use a corresponding optimal transformer (on the right), but the composition of these transformers is not optimal.For the sake of simplicity, all variables are considered to be real numbers.

Figure 5 :
Figure 5: The least fixed point representable in the domain (lfp(α • f • γ)) is not necessarily the least approximation of the least fixed point (α(lfp f )) inside the abstract domain.For instance, if we take a program initialised by x ∈ [−1, 1] and y = 0, and at each iteration, we rotate the point by 45 • , the least invariant is an 8point star, and its best approximation inside the abstract domain of intervals is the square [−1, 1] 2 .However, this square is not an inductive invariant: no rectangle (product of intervals) is stable under the iterations, thus there is no abstract inductive invariant within the interval domain.Using the domain of convex polyhedra, one would obtain a regular octagon.

3. 4 . 3 .
Synchronous Data Flow Example: Rate Limiter.To go back to the original problem of floating-point data in data-flow languages, let us consider the following library block: a rate limiter.As seen in Listing 8, such a block in inserted in a reactive loop, as shown below, where assume(c) stands for if (c) {} else { fail ();} and fail () aborts execution.

r n n− 10 ;
} e l s e { r e t u r n M(M( n + 1 1 ) ) ; } }

Figure 7 :
Figure 7: This formula is the result of quantifier elimination.It defines h z to be the least upper bound of xy for x ∈ [l x , h x ] and y ∈ [l y , h y ].