Automatic Modular Abstractions for Linear Constraints

arXiv:0811.0166v1 [cs.PL] 2 Nov 2008

David Monniaux
VERIMAG∗
June 27, 2008

Abstract
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.
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. Our algorithms are based on new quantifier elimination and symbolic manipulation techniques.
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.
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.

1

Introduction

In program analysis, it is often necessary to prove or infer numerical properties of programs, for
instance, in order to prove certain relationships between array indices, or to prove the absence
of overflows. Static program analysis by abstract interpretation obtains properties of variables,
or of relationships between variables, representable in an abstract domain. Examples of “classical” numerical abstract domains for numerical properties include intervals Cousot and Cousot
[1976] — to each variable x one attaches an interval [xmin , xmax ] — and convex polyhedra
Cousot and Halbwachs [1978] — conjunctions of inequalities a1 x1 + · · · + an xn ≤ c are inferred.
For each implemented numerical domain and each program instruction, the static analyzer
must provide an abstract transfer function, which maps the property before the instruction
to a safe property after the instruction (for forward analysis; the reverse is true of backward
analysis). For instance, over the intervals, z=x+y is optimally abstracted as zmax = xmax +ymax
and zmin = xmin + ymin ; the transfer functions for polyhedra are more complex. While the
designers of abstract interpreters generally strive so that the output property is “optimal”
(the interval [zmin , zmax ] defined above is the least possible one for the inclusion ordering),
∗

VERIMAG is a joint laboratory of CNRS, Université Joseph Fourier and Grenoble-INP.

1

optimality is not preserved by composition. Consider, for instance, y=x; z=x-y; with the
precondition that x ∈ [0, 1]. The interval for z, obtained from those for x and y by applying
the rules of interval arithmetics, is [−1, 1]; yet, the optimal interval is {0}. The reason for
this loss of precision is that while the computation of the interval for z from those for x and
y is locally optimal, it does not take into account the relationship between x and y.
Our initial target application was programs written in synchronous data-flow languages
such as Lustre Caspi et al. [1987], Simulink or Scade Caspi et al. [2003]. In these languages, operators are built out of elementary operators, introducing many intermediate variables. Successions of small elementary operations may also occur when analyzing low-level
code, e.g. assembly Gopan and Reps [2007], Balakrishnan and Reps [2004] or Java bytecode,
and they hamper certain static analysis methods due to the reduced size of the code window
used for transfer functions Logozzo and Fähndrich [2008]. Analyzing floating-point code at
the assembly level may actually be easier than analyzing higher-level programs, since the semantics of elementary floating-point operations are usually fairly well-defined while the definition and compiling processes of higher-level languages may leave significant leeway Monniaux
[2008b]. It is therefore important, for such applications, to be able to analyze program blocks
as a whole and not as a succession of independent operations.
In the above simple example, we could obtain better precision by using a relational abstract
domain linking the inputs and the outputs of the procedure. In general, though, the code
fragment may contain tests and loops (or, more generally, semantic fixed points), which
complicates the matter (see Sec. 3.4.3 for a short example whose semantics involves a fixed
point).
Ideally, for better precision, the analyzer should provide a (hopefully optimal) abstract
transfer function for each possible program block (fragment of code without loops). However,
the designers of the analyzer cannot include a hand-coded function for each possible program
block to be analyzed, if only because the number of possible program blocks is infinite. Also,
the user might want to use abstract domains not pre-programmed in the analyzer. We would
like that abstract transfer functions be obtained automatically from the definition of the
abstract domain and the source code (or semantics) of the program block.
In this article, we show how to automatically transform program blocks without loops
into an effective implementation of their optimal abstract transfer function. This optimal
transformer maps constraints on the block inputs to the tightest possible constraints on the
block output. This transformation is parametric in the abstract domain used: it takes as
inputs both the program block and a specification of the abstract domain, and outputs the
corresponding transfer function. The same method applies for both forward and backward
analysis by abstract interpretation, though, for the sake of simplicity, the article focuses on
forward analysis.
For short, our analysis considers the exact transition relation of loop-free program fragments as an existentially quantified formula. From that formula, it is able to compute the
optimal abstract transformer for the fragment with respect to a user-specified abstract domain, or even for the least invariant of the fragment in that abstract domain. The user may
specify any abstract domain in the wide class of template linear abstract domains Colon et al.
[2003].
Our method is based upon quantifier elimination in the theory of rational linear arithmetic.
It has long been known that this theory admitted quantifier elimination, but algorithms
remained mostly impractical. Recent improvements in SAT/SMT solving techniques have
2

made it possible to perform quantifier elimination on larger formulas Monniaux [2008a].
We also show how to obtain transfer functions for loops, which are also optimal in a certain
sense (they compute the least inductive invariant representable in the abstract domain).
In the beginning of the article, we focus on simple forward analysis of loop-free blocks,
then single loops (or single fixed points), for programs dealing with real or rational variables.
The same methods apply to integer variables, at the expense of some added abstraction.
We show in later sections how to deal with various constructions, including nested loops and
arbitrary control-flow graphs, recursive procedures and floating-point computations. Our focus
was indeed, originally, synchronous data-flow programs operating over real (for modeling) or
floating-point (for execution) variables, but we realized that the same technique could apply
to a wider spectrum of languages.
Our analysis goes further than most constraint-based static analysis Sankaranarayanan et al.
[2005, 2004] in that it computes the general form of the optimal postcondition or least inductive invariant as a function of the precondition parameters, not just for specific values of those
parameters. For a simple example, if the procedure is invoked on the interval domain and the
z := x + y operation, our transformation outputs zmin := xmin + ymin and zmax := xmax + ymax .
This is especially important since the function mapping the input parameters to the output
parameters may be non convex (a simple example is the abstraction of the absolute value
with respect to intervals from Sec. 3.2).
In the above case, the abstract transfer function is linear, but in general it is only piecewise
linear. It can be expressed as a simple executable program, consisting only of tests and
assignments (see an example at the end of Sec. 3.2). The analysis thus amounts to a program
transformation from the concrete to the abstract program. An advantage of obtaining the
abstract transfer functions in such a form is that it can be compiled as an ordinary program
and loaded back into the analyzer for maximal efficiency. The abstract transfer function
obtained by the analysis of a block may be retained for future use, since it is valid in any
context. An application of our transformation is therefore modular interprocedural analysis.
We have so far considered analyzes where the constraints apply to program variables
at a given control point. It is also possible to consider relationships between variables at
two different control points, especially the entry and exit of procedures. This way, we can
also analyze programs with recursive procedures, including the famous McCarthy 91 function Manna and McCarthy [1969], Manna and Pnueli [1970].
Contrary to most analyzes of numerical properties based on abstract interpretation, our
analysis for loops does not use widening operators for finding over-approximations of least
fixed points. For instance, the set of reachable states at the start of a loop (a loop invariant)
is expressed as the least fixed point of the transition relation that contains the input precondition. In widening-based analyzes, over-approximations of the set of reachable states after
1, 2, 3, etc. loop iterations are computed, and the analyzer tries to extrapolate these results
in order to obtain some “candidate” for being a loop invariant. For instance, an abstract
analyzer based on intervals may obtain [1, 2], [1, 3], [1, 5], and, because the lower bound of the
interval stays stable and the upper bound is unstable, may try [1, +∞[. If [1, +∞[ is stable
under the transition relation, then it is a safe invariant, otherwise further widening is needed.
Widenings are a major source of imprecision in many static analyzers and their design is
somewhat of a “black art”. While the soundness of the transition relation and the stability
test ensure that the analysis results are correct, and the correct construction of the widening

3

operator ensures termination, the quality of the over-approximation obtained (whether it is
close to the actual least invariant or far from it) depends on various factors. In contrast, our
method is guaranteed to yield least inductive invariants.
In Sec. 2, we recall facts of formulas built out of linear inequalities. In Sec. 3.1 we define
the class of abstract domains that we consider. In Sec. 3.2, we show how we obtain optimal
abstract transformers as logical formulas, and in Sec. 3.3 how to compile these formulas into
executable functions. In Sec. 3.4 we show how the same process applies to least inductive
invariants. In Sec. 4 we show how to deal with various extensions to the admissible domains
and operations: how to allow infinite values for constraint parameters, how to allow some
class of non-convex domains, how to partition the state space, and how to model floatingpoint computations using real numbers. In Sec. 5 we shall see how to deal with recursive
procedures and arbitrary control-flow graphs.

2

Linear formulas

We consider logical formulas built out of linear inequalities. A linear expression is a sum
a1 v1 + · · · + an vn where the ai ∈ Q and the vi are variables. Q denotes the field of rational
numbers, R the field of real numbers. A linear inequality is of the form l > 0 or l ≥ 0, where
l is a linear expression. Linear inequalities can always be scaled so that they use only integer
coefficients, as opposed to rationals. a ≤ b ≤ c is shorthand for a ≤ b ∧ b ≤ c. Unquantified
formulas are built out of atomic formulas (linear inequalities) using logical connectives ∧
and ∨. l = 0 means l ≥ 0 ∧ l ≤ 0. A formula is said to be in disjunctive normal form (DNF) if
it is written as a disjunction C1 ∨· · ·∨Cn , where each of the Ci is a conjunction Ai,1 ∧· · ·∧Ai,nj
where the Ai,j are atomic formulas or negations thereof. Quantified formulas are built out of
the same, plus the universal and existential quantifiers ∀ and ∃.
The Q-models (respectively, R-models) of a formula F are mappings m from the free
variables of F to Q (respectively, R) such that m verifies the formula; we then note m |= F .
F is said to be true if every assignment is a model (a model is a mapping from the set of
variables to Q or R), satisfiable if it has a model, and false or unsatisfiable otherwise. Truth
and satisfiability are equivalent if F has no free variables.
We say that two formulas F and G with the same free variables are equivalent, noted
F ≡ G, if they have the same models. Any formula is equivalent to a formula in disjunctive
normal form, which can be obtained by repeated application of distributivity: a ∧ (b ∨ c) ≡
(a ∧ b) ∨ (a ∧ c). F is said to imply G, noted F ⇛ G, if all models of F are models of G. We
say that F and G are equivalent modulo assumptions T , noted F ≡T G, if F ∧ T ≡ G ∧ T ;
we define similarly F ⇛T G as F ∧ T ⇛ G ∧ T . Equivalences modulo assumptions are often
used when simplifying formulas. For instance, if we know that a certain program is always
△
used in a context where T = a < b holds, and program analysis, at some point, generates the
△
△
formula F = ∃x a ≤ x ≤ b, then this formula can be simplified to G = true.
The theory of linear inequalities admits quantifier elimination: for any formula F with
quantifiers, there exists a formula G without quantifiers such that G ≡ F . There exist several
algorithms that compute such a G from F . Ferrante and Rackoff [1975] proposed a doubly
exponential method [Bradley and Manna, 2007, Sec. 7.3], which is too slow in practice; we
have since proposed another algorithm that takes advantage of the recent improvements in
satisfiability testing technology. Monniaux [2008a] Our algorithm also allows conversion to
4

disjunctive normal form, and formula simplification modulo assumptions.

3

Optimal Abstraction over Template Linear Constraint Domains

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 p1 , . . . , pn and m state variables s1 , . . . , sm . We note
γF : Qn → P(Qm ) defined by γF (~
p) = {~s ∈ Qm | (~
p, ~s) |= F }. As an example, the interval
abstract domain for 3 program variables s1 , s2 , s3 uses 6 parameters m1 , M1 , m2 , M2 , m3 , M3 ;
the formula is m1 ≤ s1 ≤ M1 ∧ m2 ≤ s2 ≤ M2 ∧ m3 ≤ s3 ≤ M3 .
In this section, we focus on the case where F is a conjunction L1 (s1 , . . . , sm ) ≤ p1 ∧
· · · ∧ Ln (s1 , . . . , sm ) ≤ pn of linear inequalities whose left-hand side is fixed and the righthand sides are parameters. Such conjunctions define the class of template linear constraint
domains Colon et al. [2003]. Particular examples of abstract domains in this class are:
• the intervals (for any variable s, consider the linear forms s and −s);
• the difference bound matrices (for any variables s1 and s2 , consider the linear form
s1 − s2 );
• the octagon abstract domain (for any variables s1 and s2 , distinct or not, consider the
linear forms ±s1 ± s2 ) Miné [2001]
• the octahedra (for any tuple of variables s1 , . . . , sn , consider the linear forms ±s1 · · · ±
sn ). Clarisó and Cortadella [2004]
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 s1 ≤ p1 ∧s2 ≤ p2 ∧s1 +s2 ≤ p3 , then both (p1 , p2 , p3 ) = (1, 1, 2)
and (1, 1, 3) define the same set for state variables s1 and s2 . 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 ⊆ Qm . Take for all i = 1, . . . , m: pi = sup~s∈W Li (~s).
Clearly, W ⊆ γF (p1 , . . . , pm ), and in fact p~ is such that γF (~
p) is the least solution to this
inclusion. pi belongs in general to R ∪ {+∞}, not√necessarily to Q ∪ {+∞}. (for instance,
if W = {s1 | s21 ≤ 2} and L1 = s1 , then p1 = 2). We have therefore defined an αF :
P(Rm ) → {⊥} ∪ (R ∪ {+∞})n , and (αF , γF ) form a Galois connection: αF maps any set to
its best upper-approximation. The fixed points of αF ◦γF are the normal forms. For instance,
s1 ≤ 1 ∧ s2 ≤ 1 ∧ s1 + s2 ≤ 2 is in normal form, while s1 ≤ 1 ∧ s2 ≤ 1 ∧ s1 + s2 ≤ 3 is not.

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 consisting
in linear arithmetic assignments, linear tests, and sequences. 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 JP K such that (a, b, . . . , a′ , b′ , . . . ) |= P if
5

and only if the memory state (a′ , b′ , . . . ) can be reached at the end of an execution starting
in memory state (a, b, . . . ):
△

Arithmetic Ja := L(a, b, . . . ) + KKF = a′ = L(a, b, . . . ) + K ∧ b′ = b ∧ c′ = c ∧ . . . where K
is a real constant and L is a linear form, and b, c, d . . . are all the variables except a;
△

Tests Jif c then p1 else p2 K = (c ∧ Jp1 KF ) ∨ (¬c ∧ Jp2 KF );
△

Non deterministic choice Ja := randomK = b′ = b ∧ c′ = c ∧ . . . , for all variables except a;
△

Failure JfailK = false;
△

Skip JskipK = a′ = a ∧ b′ = b ∧ c′ = c ∧ . . .
△

Sequence JP1 ; P2 KF = ∃a′′ , b′′ , . . . f1 ∧ f2 where f1 is JP1 KF where a′ has been replaced by
a′′ , b′ by b′′ etc., f2 is JP2 KF 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. This is why we defer quantifier elimination to the point where we
compute the abstract transfer relation.
△
Consider now a domain definition formula F = L1 (s1 , s2 , . . . ) ≤ p1 ∧ · · · ∧ Ln (s1 , s2 , . . . ) ≤
△
pn on the program inputs, with parameters p~ and free variables ~s, and another F ′ = L′1 (s′1 , s′2 , . . . ) ≤
p′1 ∧ · · · ∧ L′n (s′1 , s′2 , . . . ) ≤ p′n 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
∀~s, s~′ F ∧ JP K =⇒ F ′

(1)

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: Ci′ = ∀~s, s~′ F ∧ JP K =⇒ L′i (s~′ ) ≤ p′i .
△
Let us take a simple example: if P is the program instruction z := x + y, F = x ≤
△
△
p1 ∧ y ≤ p2 , F ′ = z ≤ p′1 , then JP K = z ′ = x + y, and the soundness condition is ∀x, y, z (x ≤
p1 ∧ y ≤ p2 ∧ z = x + y =⇒ z ≤ p′1 ). Remark that this soundness condition is equivalent
to a formula without quantifiers p′1 ≥ p1 + p2 , 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 p1 = p2 = 1), only one value is optimal (p′1 = 2 for p1 = p2 = 1). An
△
optimal value for the output parameter p′i is defined by Oi′ = Ci′ ∧ ∀qi′ (Ci′ [qi′ /p′i ] =⇒ p′i ≤ qi′ ).
Again, quantifier elimination can be applied; on our simple example, it yields p′1 = p1 + p2 .
6

If there are n input constraint parameters p1 , . . . , pn , then the optimal value for each output
constraint parameter p′i is defined by a formula Oi′ with n + 1 free variables p1 , . . . , pn , p′i . This
formula defines a partial function from Qn to Q, in the mathematical sense: for each choice
of p1 , . . . , pn , there exist at most a single p′i . The values of p1 , . . . , pn 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:
if (x >= 10) { y = random; } else { y = 0; }
If F = x ≤ p1 and F ′ = y ≤ p′1 , then O1′ ≡ p1 < 10 ∧ p′1 = 0, and the function is defined only
for p1 < 10.
At this point, we have a characterization 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) Oi′ 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 ∈ [xmin , xmax ] and the post-condition is y ∈ [ymin , ymax ]. Acceptable values for (ymin , ymax )
are characterized by formula
△

G = ∀x xmin ≤ x ≤ xmax =⇒ ymin ≤ |x| ≤ ymax

(2)

′
′
′
. Quantifier
/ymax ] =⇒ ymax ≤ ymax
G[ymax
The optimal value for ymax is defined by G∧∀ymax
elimination over this last formula gives as characterization for the least, optimal, value for
ymax :

(xmin + xmax ≥ 0 ∧ ymax = xmax )∨

(xmin + xmax < 0 ∧ ymax = −xmin ). (3)

We shall see in the next sub-section that such a formula can be automatically compiled into
code such as:
if (xmin + xmax >= 0) {
ymax = xmax;
} else {
ymax = -xmin;
}

3.3

Generation of the Implementation of the Abstract Domain

Consider formula 3, 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, xmin + xmax < 0 and xmin + xmax ≥ 0), and
on each of these regions, the output parameter is a linear function of the input parameters.
Given a disjunct (such as ymax = −xmin ∧ xmin + xmax < 0), the domain of validity of the
disjunct can be obtained by existential quantifier elimination over the result variable (here
∃ymax (ymax = −xmin ∧xmin +xmax < 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
7

overlap, but in this case, since Oi′ defines a function in the mathematical sense, the functions
defined by such disjuncts coincide on their overlapping domains of validity.
This suggests a first algorithm for conversion to an executable form:
1. Put Oi′ into quantifier-free, disjunctive normal form C1 ∧ · · · ∧ Cn .
2. For each disjunct Ci , obtain the validity domain Vi as a conjunction of linear inequalities
and solve for p′i (obtain p′i as a linear function vi of the p1 , . . . , pn ).
3. Output the result as a cascade of if-then-else and assignments, as in the example at the
end of Sec. 3.2.
Algorithm 1 ToITEtree(F, z, T ): turn a formula defining z as a function of the other free
variables of F into a tree of if-then-else constructs, assuming that T holds.
D(= C1 ∧ · · · ∧ Cn ) ← QElimDNFModulo({}, F, T )
for all Ci ∈ D do
Pi ← QElimDNFModulo(z, F, T )
end for
P ← Predicates(P1 , . . . , Pn )
if P = ∅ then
Ensure: ∃z F is always true
O ← Solve(D, z)
else
K ← Choose(P )
O ← IfThenElse(K, ToITEtree(F, z, T ∧ K), ToITEtree(F, z, T ∧ ¬K))
end if
An if-then-else cascade may be inefficient, since identical conditions may have to be tested
several times. We could of course factor out all conditions and assign them to Boolean
variables, but then, some of the tests performed may actually not be needed. We therefore
propose an algorithm for building an if-then-else tree. 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 Oi′ 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. One
then solves this conjunction to obtain the computed variable (here, ymax ) as a function.
Our algorithm ToITEtree(F, z, T ) (Alg. 1) uses a function QElimDNFModulo(~v , F, T )
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. Predicates(F ) returns the set of atomic predicates of F . Solve(D, z)
solves a minimal disjunction D of inequalities for variable z, assuming that there is at most
one solution for z for each choice of the other variables; one simple way to do that is to look
8

for any constraint of the form z ≥ L or z ≤ L and output z = L. Choose(P ) chooses any
predicate in P (one good heuristic seems to be to choose the most frequent in P1 , . . . , Pn ).
Let us take, as a simple example, formula 3. We wish to obtain ymax as a function of
△
xmin and xmax , so in the algorithm ToITEtree we set z = ymax . C1 is the first disjunct
xmin + xmax ≥ 0 ∧ ymax = xmax , C2 is the second disjunct xmin + xmax < 0 ∧ ymax = −xmin .
We project C1 and C2 parallel to ymax , obtaining respectively P1 = (xmin + xmax ≥ 0) and
P2 = (xmin + xmax < 0). We choose K to be the predicate xmin + xmax ≥ 0 (in this case, the
choice does not matter, since P1 and P2 are the negation of each other).
△

• The first recursive call to ToITEtree is made in the context of T = (xmin + xmax ≥ 0).
Obviously, F ∧ T ≡ (ymax = xmax ) ∧ T and thus (∃ymax F ) ∧ T ≡ T .

QElimDNFModulo(ymax , F, T ) will then simply output the formula “true”. It then
suffices to solve for ymax in ymax = xmax . This yields the formula for computing the
correct value of ymax in the cases where xmin + xmax ≥ 0.
△

• The second recursive call is made in the context of T = (xmin + xmax < 0. The result is
ymax = −xmin , the formula for computing the correct value of ymax in the cases where
xmin + xmax < 0.
These two results are then reassembled into an 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.

3.4

Least Inductive Invariants

We have so far considered programs without loops. We shall now see that not only can we
compute the optimal abstract post-condition of a block as a simple, executable function of the
parameters of the precondition, but we can also compute the parameters of the least inductive
invariant of a program block that is of the form specified by the abstract domain.1 Beware
that this least inductive invariant found in the abstract domain is in general different from the
least element of the abstract domain that includes the least inductive invariant of the system
(Fig. 1).
3.4.1

Stability Inequalities

Consider a program fragment: while (c) { p; }. We have domain definition formulas
△
F = L1 (s1 , . . . , sm ) ≤ p1 ∧ · · · ∧ Ln (s1 , . . . , sm ) ≤ pn for the precondition of the program
△
fragment , and F ′ = L′1 (s1 , . . . , sm ) ≤ p′1 ∧ · · · ∧ L′n (s1 , . . . , sm ) ≤ p′n for the invariant.
1

In order to specify the least invariant, we would have to quantify over all sets of states, then filter those
which are inductive invariants. This is second-order quantification, which we cannot handle. By restricting
ourselves to invariants of a certain shape, we replace it by first order quantification.

9

Figure 1: 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 initialized by x ∈ [−1, 1] and y = 0, and at each iteration,
we rotate the point by 45◦ , the least invariant is an 8-point star, and the 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.

10

Define G = JcK ∧ JpK. G is a formula whose free variables are s1 , . . . , sm , s′1 , . . . , s′m such
that (s1 , . . . , sm , s′1 , . . . , s′m ) |= G if and only if the state (s′1 , . . . , s′m ) can be reached from
the state (s1 , . . . , sm ) in exactly one iteration of the loop. A set W ⊆ Qm 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:
∀~s, s~′ F ′ ∧ G =⇒ F ′ [s~′ /~s].

(4)

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
△
H = (∀~s, F =⇒ F ′ ) ∧ (∀~s, s~′ F ′ ∧ G =⇒ F ′ [s~′ /~s])

(5)

This formula links the values of the input constraint parameters p1 , . . . , pn 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
(∃p′1 , . . . , p′i−1 , p′i+1 , . . . , p′n H) ∧ (∀q~′ H[q~′ /p~′ ] =⇒ p′i ≤ qi′ )

(6)

The free variables of this formula are p1 , . . . , pn , p′i . This formula defines a function (in the
mathematical sense) defining p′i from p1 , . . . , pn . As before, this function can be compiled to
an executable version using cascades or trees of tests.
3.4.2

Simple Loop Example

To show how the method operates in practice, let us consider first a very simple example
(something happens is a nondeterministic choice):
int i=0;
while (i <= n) {
if (something_happens) {
i=i+1;
if (i == n) {
i=0;
}
}
}

11

Let us abstract i at the head of the loop using an interval [imin , imax ]. 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 6= 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:
imin ≤ 0 ∧ 0 ≤ imax ∧ ∀i∀i′ ((imin ≤ i ∧ i ≤ imax ∧

(((i + 1 ≤ n − 1 ∨ i + 1 ≥ n + 1) ∧ i′ = i + 1)∨

(i + 1 = n + 1 ∧ i′ = 0) ∨ i′ = i)) =⇒ (imin ≤ i′ ∧ i′ ≤ imax )) (7)

Quantifier elimination produces:
(imin ≤ 0 ∧ imax ≥ 0 ∧ imax < n ∧ −imin + n − 2 < 0)∨

(imin ≤ 0 ∧ imax ≥ 0 ∧ imax − n + 1 ≥ 0 ∧ imax < n) (8)

The formulas defining optimal imin and imax are:
imin ≥ 0 ∧ imin ≤ 0 ∧ n > 0

(imax = 0 ∧ ∧n > 0 ∧ n < 2) ∨ (imax = n − 1 ∧ imax ≥ 1)

(9)
(10)

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:
if (n <= 0) {
fail();
} else {
iMin = 0;
if (n < 2) {
iMax = 0;
} else /* n >= 2 */
iMax = n-1;
}
}
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
Cousot and Halbwachs [1978] or octagons Miné [2001].
Let us now consider the same program, simply replacing n by the constant 20. All implementations of intervals (and thus of octagons and polyhedra, since we only have one variable),
will overshoot the imax = 19 target when using the traditional widening and narrowing strategies: they will compute i ∈ [0, 0], then ∈ [0, 1], ∈ [0, 2] and widen to [0, +∞[, and narrowing
will not reduce the interval. Even if we replaced i == 20 by i >= 20, narrowing would still
fail to reduce the interval due to the nondeterministic choice since the concrete transfer function f , mapping sets of states at the head of the loop to sets of states at the next iteration, is
12

expansive: for all set of states W , W ⊆ f (W ). This is a well-known weakness of the widening/narrowing approach, and the workaround is a syntactic trick known as widening up to or
widening with thresholds: for all variables, the constants to which it is compared are gathered
and used as widening steps [Blanchet et al., 2003, Sec. 7.1.2]. This syntactic approach fails if
tests are more indirect, whereas our semantic approach is not affected.
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. When compiled into C, such a block in inserted in
a reactive loop, as shown below, where assume(c) stands for if (c) {} else {fail();}:
while (true) {
...
e1 = random(); assume(e1 >= e1min && e1 <= e1max);
e2 = random(); assume(e2 >= e2min && e2 <= e2max);
e3 = random(); assume(e3 >= e3min && e3 <= e3max);
olds1 = s1;
if (random) {
s1 = e3;
} else {
if (e1 - olds1 < -e2) {
s1 = olds1 - e2;
}
if (e1 - olds1 > e2) {
s1 = olds1 + e2;
}
}
...
}
We are interested in the input-output behavior of that block: obtain bounds on the output
s1 of the system as functions of bounds on the inputs (e1, e2, e3). Note that in this case, s1,
e1, e2, e3 are streams, not single scalars. One difficulty is that the s1 output is memorized,
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 s1min ≤ s1 ≤ s1max under the
assumption that e1min ≤ e1max ∧e2min ≤ e2max ∧e3min ≤ e3max . The stability condition yields,
after quantifier elimination and projection on s1max the condition s1max ≥ e1max ∧ s1max ≥
e3max . Minimization then yields an expression that can be compiled to an if-then-else tree:
if (e1max > e3max) {
s1max = e1max;
} else {
s1max = e3max;
}

13

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 behavior if all variables and
operations are IEEE-754 floating-point, since rounding errors introduce slight differences of
regimes between ranges of inputs (Sec. 4.4, 6). Rounding errors in the program to be analyzed
introduce difficulties for analyzes 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 [Blanchet et al., 2003, Sec. 7.1.4].

4

Extensions to the Admissible Domains and Operations

The class of domains and program constructs of the preceding section may seem too limited.
We shall see here a few extensions.

4.1

Infinities

Consider the interval abstract domain, defined by x ≤ p2 ∧−x ≤ p1 . The techniques explained
in Sec. 3.1 allow only finite bounds. Yet, it makes sense that p1 and p2 could be equal to
+∞ so as to represent infinite intervals. This can be easily achieved by a minor alteration
∞
to our definitions. Each parameter pi is replaced by two parameters pbi and p∞
i . pi 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); pi = 0 means that pi is finite and
b
∞
∞
equal to pi , pi = 1 means pi = +∞. Li ≤ pi becomes (pi > 0) ∨ (Li ≤ pbi ), Li < pi becomes
b
(p∞
i > 0) ∨ (Li < pi ). After this rewriting, all formulas are formulas of the theory of linear
inequalities without infinities and are amenable to the appropriate algorithms.

4.2

Non-Convex Domains

Section 3.1 constrains formulas to be conjunctions of inequalities of the form Li ≤ pi . 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
T (infinite) intersection; that is, if we have a family p~ ∈ W , then there
p) = γF (~
p0 ) (besides, p0 = inf{~
p | p~ ∈ W }). This is what enables
exist p0 such that p~∈W γF (~
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 ((−x ≤ p1 ∧ x ≤ p2 ) ∨ (−x ≤ p3 ∧ x ≤ p4 )) ∧ p2 ≤ p3 . There are two,
incomparable, minimal elements of the form γF (~
p): p1 = p2 = 0 ∧ p3 = −1 ∧ p4 = 2 and
p1 = 0 ∧ p2 = 1 ∧ p3 = −2 ∧ p4 = 2.
We consider formulas F built out of linear inequalities Li (s1 , . . . , sn ) ≤ pi as atoms,
conjunctions, and disjunctions. By induction on the structure of F , we can show that γF :
(R ∪ {−∞})n → P(Rn ) is inf-continuous; that
pi )i∈I such that
T is, for any descending chain (~
limi p~i = p~∞ , then γF (~
pi ) is decreasing and i∈I γF (~
pi ) = γF (~
p∞ ). The property is trivial for
atomic formulas, and is conserved by greatest lower bounds (∧) as well as binary least upper
bounds (∨).
14

Let us consider a set S ⊆ P(Rn ), 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 . γF−1 (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 . γF−1 (S) is stable under greatest lower
bounds of descending chains: take a descending chain (~
pi )i∈I , then γF (limi ~pi ) = ∩i γF (~
pi ) ∈ S
−1
by inf-continuity and stability of S. By Zorn’s lemma, γF (S) has at least one minimal
element.
Let P [~
p] be a formula representing γF−1 (S) (Sec. 3.1 proposes formulas defining safe post△
conditions and inductive invariants). The formula G[~
p] = P [~
p] ∧ ∀p~′ P [p~′ ] ∧ p~′ ≤ p~ =⇒ p~ ≤ p~′
−1
defines the minimal elements of γ (S).
△
For instance, consider ~
p = (a, b, c, d), F = (−x ≤ a ∧ x ≤ b) ∨ (−x ≤ c ∧ x ≤ d),
representing unions of two intervals [−a, b] ∪ [−c, d]. We want upper-approximations of the
△
set {0, 1, 3}; that is P [~
p] = ∀x (x = 0 ∨ x = 1 ∨ x = 3 =⇒ F [~
p, x]). We add the constraint
that −a ≤ b ∧ b ≤ −c ∧ −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 ∧ b = 1 ∧ c = −3 ∧ d = 3) ∨ (a = 0 ∧ b = 0 ∧ c = −1 ∧ d = 3), that is, either
[0, 0] ∪ [1, 3] or [0, 1] ∪ [3, 3].

4.3

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.
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. Consider pairwise disjoint subsets (Ci )i∈I of the
state space Qm , and abstract domains stable under
S intersection (Si )i∈I , S
Si ⊆ P(Ci ). Elements
of the partitioned abstract domain are unions i∈I si where si ∈ Si . If ( i si,j ])j∈J is a family
 S T
T
S
of elements of the domain, then j∈J
i∈I
j∈J si,j ; that is, intersections are
i∈I si,j ] =
taken separately in each Ci .
Take a family (Fi [~
p])i∈I of formulas defining template linear constraint domains (conjunctions of linear inequalities Li (s1 , . . . , sn ) ≤ pi ) and a family (Ci )i∈I of formulas such
that for all i and i′ , Ci ∧ Ci′ is equivalent to false and C1 ∨ · · · ∨ Cl is equivalent to true.
F = (C1 ∧ F1 ) ∨ · · · ∨ (Cl ∧ Fl ) then defines an an abstract domain such that γF is a infmorphism. All the techniques of Sec. 3.1 then apply.

4.4

Floating-Point Computations

Real-life programs do not operate on real numbers; they operate on fixed-point or floatingpoint 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.
√
for
In IEEE floating-point IEE [1985], each atomic operation (noting ⊕, ⊖, ⊗, ⊘,
f
√
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

15

function.2 This rounding function, depending on user choice, maps each real x to the nearest
floating-point value rn (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 +∞), r0 (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
The semantics of the rounding operation cannot be exactly represented inside the theory of
linear inequalities.3 As a consequence, we are forced to use an axiomatic over-approximation
of that semantics: a formula linking a real number x to its rounded version r(x).
Miné [2004] 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.
Because our language for axioms is richer than the interval linear forms used by Miné, 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 = ±s.m where 1 ≤ m < 2 is the mantissa or significand, which has a fixed number
p of bits, and s = 2e the scaling factor (Emin ≤ e ≤ Emax 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 normalized representation.
Except in the case of numbers of very small magnitude, IEEE-754 always works with normalized representations. There exists a least positive normalized number mnormal and a least
positive denormalized number mdenormal , and the denormals are the multiples of mdenormal
less than mnormal . All representable numbers are multiples of mdenormal .
Consider for instance floating-point addition or subtraction x = ±a ± b. Suppose that
0 ≤ x ≤ mnormal . a and b are multiples of mdenormal and thus a − b is exactly represented as a
denormalized number; therefore r(x) = x. If x > mnormal , then |r(x) − x| ≤ εrel .x. The cases
for x ≤ 0 are symmetrical. We can therefore characterize r(x) − x using linear inequalities
through case analysis over x: Round+ (a ⊕ b, a + b) (respectively, Round+ (a ⊖ b, a − b)) holds,
where
△

Round+ (r, x) = (x ≤ mnormal ∧ r = x)

∨ (x > mnormal ∧ −εrel .x ≤ r − x ≤ εrel .x (11)

2

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. Monniaux [2008b]
3
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.

16

△

Round(r, x) = (x = 0 ∧ r = 0)∨

(x > 0 ∧ r ≥ 0 ∧ Round+ (r, x))∨

(x < 0 ∧ r ≤ 0 ∧ Round+ (−r, −x)) (12)

To each floating-point expression e, we associated a “rounded-off” variable re , the value of
which we constrain using Round(re , e) or Round+ (re , e). For instance, a expression e = a ⊕ b
is replaced by a variable re , and the constraint Round+ (re , a + b) is added to the semantics.
In the case of a compound expression e = ab + c, we introduce e1 = ab, and we obtain
Round+ (re , re1 + c) ∧ Round(re1 , ab). If we know that the compiler uses a fused multiply-add
operator, we can use Round(re , ab + c) instead.

5

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

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:
if (e1max > e3max) {
s1max = e1max;
} else {
s1max = e3max;
}
assume(s1 <= s1max);
/* and similar for 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. [2008, §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 [Gulwani et al.,
2008, §2.3]; Bourdoncle [1992, 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, parameterized by a number of variables. The values of these variables for all points in the cut-set
defines 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
17

a stability condition (Formula 4). The conjunction of all these stability conditions defines
acceptable inductive invariants. As above, the least inductive invariant is obtained by writing
a minimization formula (Sec. 3.4).
Let us take a simple example:
i=0;
while(true) { /* A */
if (choice()) {
j=0;
while(j < i) { /* B */
/* something */
j=j+1;
}
i=i+1;
if (i==20) {
i=0;
}
} else {
/* something */
}
}
We choose program points A and B as cut-set. At program point A, we look for an
△
invariant of the form IA (i, j) = imin,A ≤ i ≤ imax,A , and at program point B, for an invariant
△
of the form IB (i, j) = imin,B ≤ i ≤ imax,B ∧ jmin ≤ j ≤ jmax ∧ δmin ≤ i − j ≤ δmax (a
difference-bound invariant). The (somewhat edited for brevity) stability formula is written:

∀j IA (0, j) ∧ ∀i∀j ((IB (i, j) ∧ j ≥ i ∧ (i + 1 ≤ 19∨

i + 1 = 20 ∨ i + 1 ≥ 21)) ⇒ If[i + 1 = 20, IA (0, j), IA (i + 1, j)])∧
∀i∀j (IA (i, j) ⇒ IB (i, 0)) ∧ ∀i∀j ((IB (i, j) ∧ j < i)

⇒ IB (i, j + 1)) (13)

Replacing IA and IB into this formula, then applying quantifier elimination, we obtain a
formula defining all acceptable tuples (imin,A , imax,A , imin,B , imax,B , jmin , jmax , δmin , δmax ). Optimal values are then obtained by further quantifier elimination: imin,A = imin,B = jmin = 0,
imax,A = imax,B = 19, jmax = 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.
18

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. Astrée Blanchet et al.
[2002, 2003], Cousot et al. [2005]. 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
(mz ≤ z ≤ Mz ); suppose we obtain Mz := 1000; mz := 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. This method is however insufficient for dealing with recursive procedures.
In order to analyze 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” Manna and McCarthy
[1969], Manna and Pnueli [1970], which, non-obviously, returns 91 for all inputs less than 101:
int M(int n) {
if (n > 100) {
return n-10;
} else {
return M(M(n+11));
}
}
The concrete semantics of that function is a relationship R between its input n and its
output r. It is the least solution of
R ⊇ {(n, r) ∈ Z2 | (n > 100 ∧ r = n − 10)∨

(n ≤ 100 ∧ ∃n2 ∈ Z(n + 11, n2 ) ∈ R ∧ (n2 , r) ∈ R)} (14)
△

We look for a inductive invariant of the form I = ((n ≥ A) ∧ (r − n ≥ δ) ∧ (r − n ≤
∆)) ∨ ((n ≤ B) ∧ (r = C)), a non-convex domain (Sec. 4.2). By replacing R by I into
inclusion 14, and by universal quantification over n, r, n2 , 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 optimize, 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 Martin et al. [1998].

19

6

Implementations and Experiments

We have implemented the techniques of Sec. 3 in quantifier elimination packages, including
Mathematica4 and Reduce 3.85 + Redlog6 in addition to our own package, Mjollnir Monniaux [2008a].7
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 Necula et al. [2002] 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.3 was 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 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.
Analyzers such as Astrée Blanchet et al. [2002, 2003], Cousot et al. [2005] 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 Astrée development team therefore had to provide specialized, hand-written analyzes. In contrast, all linear floating-point
operators in the library were analyzed 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.4, 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 analyzer can therefore cache this information for further use and need not recompute abstractions for library functions or operators unless these
functions are updated.
Our analyzer 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-command programs that e.g. Astrée accepts, which do not consist
purely of numerical operators. We plan to integrate our analysis method into a more generic
analyzer. 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
4

http://www.wolfram.com/
http://www.uni-koeln.de/REDUCE/
6
http://www.algebra.fim.uni-passau.de/~ redlog/
7
Source code and GNU/Linux/IA32 binaries of this implementation
http://www-verimag.imag.fr/~ monniaux/download/automatic_abstraction.zip.
5

20

are

available

from

application. For instance, our implementation 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 analyze 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 Monniaux [2008a].

7

Related Works

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 Halbwachs [1979], Cousot and Halbwachs [1978], and there have been since
then many improvements to the technique; a bibliography was gathered by Bagnara et al.
[2006]. Algorithms on polyhedra are costly and thus a variety of domains intermediate between
simple interval analysis and convex polyhedra were proposed Miné [2001], Clarisó and Cortadella
[2004], Sankaranarayanan et al. [2005]. All these domains compute invariants using a widening
operator Cousot and Cousot [1976], Cousot and Halbwachs [1978], Cousot and Cousot [1992].
There is, however, no guarantee that the resulting invariant is the best representable in the
abstract domain, even with the use of narrowing iterations; this is one difference with our
proposal, which computes the best representable inductive invariant.
Another difference is that these domains are designed to work with numerical values for the
input constraints, thus the computation must be done for every value of the input constraints
parameters. Using simple program transformations, they may also apply to symbolic input
constraints (constraint parameters being taken as extra variables), but in general this will lead
to bad results; for instance, the input-output relationship for the rate limiter of Sec. 3.4.3 is
not convex, while numerical abstract domains in the literature are convex. In comparison the
algorithm in this article can be run once to obtain a formula that gives the best invariant
depending on the input constraints, allowing modular analysis.
Several methods have been proposed to synthesize invariants without using widening operators Colon et al. [2003], Cousot [2005], Sankaranarayanan et al. [2004]. 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. [2005a,b] proposed improvements over the “classical” widenings on linear constraint domains Halbwachs [1979]. Gopan and Reps [2006] 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 programs; the idea
is to distinguish phases or modes of operation in order to make the widening more precise. Gonnord and Halbwachs [2006] have proposed acceleration techniques for linear constraints. These do not replace widenings altogether, but they alleviate the need for some of
the costly workarounds to the imprecision introduced by widenings, such as delayed widening [Blanchet et al., 2003, Sec. 7.1.3]. These address a different problem from ours. On the
one hand, neither improved widenings nor acceleration guarantee that the inductive invariant

21

obtained at the end is the least one (indeed, they can yield the top element ⊤). 8 Furthermore, the invariant that these methods obtain is not parametric in the precondition, contrary
to the one that our method obtains. On the other hand, improved widenings work regardless
of the form of the transition relation, which our method constrains to be piecewise linear.
Some of the cited methods operate on general polyhedra, while our method constrains the
shape of the polyhedra that are found to a certain template.
Gaubert et al. [2007], Gawlitza and Seidl [2007] proposed replacing the usual widening/narrowing
iteration techniques by a policy iteration (or strategy iteration) approach. Their approach converges on a fixed point, but not necessarily the least one. Their idea is to replace computing
the least fixed point of a complex abstract operator (the point-wise minimum of a family of
simpler operators) by a sequence of least fixed point computations for these simple operators.
Their technique anyway needs to compute these latter least fixed points, and it is possible
that our method can help in that respect.
Techniques using quantifier elimination for generating nonlinear invariants for programs
using nonlinear arithmetic have also been proposed Kapur [2004] and shown capable of producing optimal invariants parameterized by input constraints Monniaux [2007]. Quantifier
elimination in the theory of real closed fields is, however, a very costly technique. Experimentally, the formulas generated by common implementations tend to grow huge (due to difficult
simplifications) and both time and space requirements grow very fast with the number of
variables. This is why we considered the linear case in the present article.
Gulwani et al. [2008] 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. Quantifier elimination over Presburger arithmetic is however
very expensive Fischer and Rabin [1974]. 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.
The idea of producing procedure summaries Sharir and Pnueli [1981] as formulas mapping
input bounds to output bounds is not new. Rugina and Rinard [2005], 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 analyzed. Martin et al. [1998] proposed applying interprocedural analysis to
loops.
Seidl et al. [2007] also produce procedure summaries as numerical constraints. Our procedure summaries are implementations of the corresponding abstract transformer over some
8

There exist exact acceleration techniques but these rather apply to discrete automata.

22

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 a 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. [2007, 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. [2005] 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.
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 Monniaux [2005]. 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.

8

Conclusion and Future Work

Writing static analyzers by hand has long been found tedious and error-prone. One may of
course prove an existing analyzer 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.3 how it was possible to partition the state space and abstract each region of

23

the state-space separately; but naive partitioning according to n Booleans leads to 2n regions,
which can be unbearably costly and is unneeded in most cases. We think that automatic
refinement and partitioning techniques Jeannet [2003] could be developed in that respect.

References
Roberto Bagnara,
Patricia M. Hill,
Elena Mazzi,
and Enea Zaffanella.
Widening operators for weakly-relational numeric abstractions.
In Static Analysis
(SAS), volume 3672 of LNCS, pages 3–18. Springer, 2005a. DOI: 10.1007/11547662 3.
Roberto Bagnara,
Patricia M. Hill,
Elisa Ricci,
and Enea Zaffanella.
Precise widening operators for convex polyhedra.
Sci. Comput. Program., 58(1-2):
28–56, 2005b. DOI: 10.1016/j.scico.2005.02.003.
Roberto Bagnara, Patricia M. Hill, and Enea Zaffanella. The Parma Polyhedra Library,
version 0.9, 2006. URL http://www.cs.unipr.it/ppl.
Gogul Balakrishnan and Thomas Reps. Analyzing memory accesses in x86 executables. In
Compiler Construction (CC), volume 2985 of LNCS, pages 5–23. Springer, 2004. DOI:
10.1007/b95956.

Bruno Blanchet,
Patrick
Cousot,
Radhia Cousot,
Jérôme
Feret,
Laurent Mauborgne,
Antoine Miné,
David Monniaux,
and Xavier Rival.
Design and implementation of a special-purpose static program analyzer for safety-critical real-time embedd
In The Essence of Computation: Complexity, Analysis, Transformation, number 2566 in
LNCS, pages 85–108. Springer, 2002. DOI: 10.1007/3-540-36377-7 5.
Bruno Blanchet,
Patrick
Cousot,
Radhia Cousot,
Jérôme
Feret,
Laurent Mauborgne,
Antoine Miné,
David Monniaux,
and Xavier Rival.
A static analyzer for large safety-critical software.
In Programming Language Design
and Implementation (PLDI), pages 196–207. ACM, 2003. DOI: 10.1145/781131.781153.
François Bourdoncle. Sémantique des langages impératifs d’ordre supérieur et interprétation abstraite.
PhD thesis, École polytechnique, Palaiseau, 1992.
Aaron R. Bradley and Zohar Manna. The Calculus of Computation: Decision Procedures with
Applications to Verification. Springer, October 2007.
Paul Caspi, Daniel Pilaud, Nicolas Halbwachs, and John A. Plaice. LUSTRE: a declarative
language for real-time programming. In Principles of Programming Languages (POPL),
pages 178–188. ACM, 1987. DOI: 10.1145/41625.41641.
Paul Caspi, Adrian Curic, Aude Maignan, Christos Sofronis, Stavros Tripakis, and Peter
Niebert. From Simulink to Scade/Lustre to TTA: a layered approach for distributed embedded applications. SIGPLAN notices, 38(7):153–162, 2003. DOI: 10.1145/780731.780754.
CAV05. Computer Aided Verification (CAV), number 4590 in LNCS, 2005. Springer. DOI:
10.1007/b138445.

24

Robert Clarisó and Jordi Cortadella. The octahedron abstract domain. In Static Analysis
(SAS), number 3148 in LNCS, pages 312–327. Springer, 2004.
Michael Colon, Sriram Sankaranarayanan, and Henny Sipma. Linear invariant generation
using non-linear constraint solving. In Computer Aided Verification (CAV), number 2725
in LNCS, pages 420–433. Springer, 2003.
Patrick Cousot. Proving program invariance and termination by parametric abstraction, lagrangian relaxation
In VMCAI05, pages 1–24. DOI: 10.1007/b105073.
Patrick Cousot and Radhia Cousot. Abstract interpretation and application to logic programs.
J. of Logic Programming, 13(2–3):103–179, 1992.
Patrick Cousot and Radhia Cousot. Static determination of dynamic properties of programs.
In Proceedings of the Second International Symposium on Programming, pages 106–130.
Dunod, Paris, France, 1976.
Patrick Cousot and Nicolas Halbwachs. Automatic discovery of linear restraints among variables of a program.
In Principles of Programming Languages (POPL), pages 84–96. ACM, 1978. DOI:
10.1145/512760.512770.
Patrick Cousot, Radhia Cousot, Jérôme Feret, Laurent Mauborgne, Antoine Miné, David
Monniaux, and Xavier Rival. The ASTRÉE analyzer. In Programming Languages and
Systems (ESOP), number 3444 in LNCS, pages 21–30, 2005.
ESOP07. Programming Languages and Systems (ESOP), volume 4421 of LNCS, 2007.
Springer. DOI: 10.1007/978-3-540-71316-6.
Jeanne Ferrante and Charles Rackoff. A decision procedure for the first order theory of real
addition with order. SIAM Journal of Computation, 4(1):69–76, March 1975.
Michael J. Fischer and Michael O. Rabin. Super-exponential complexity of Presburger arithmetic.
In Complexity of Computation, number 7 in SIAM–AMS proceedings, pages 27–42. American Mathematical Society, 1974.
Stéphane Gaubert, Éric Goubault, Ankur Taly, and Sarah Zennou. Static analysis by policy iteration on relational domains. In ESOP07, pages 237–252. DOI: 10.1007/978-3-540-71316-6.
Thomas Gawlitza and Helmut Seidl. Precise fixpoint computation through strategy iteration.
In ESOP07, pages 300–315. DOI: 10.1007/978-3-540-71316-6 21.
Laure Gonnord and Nicolas Halbwachs. Combining widening and acceleration in linear relation analysis.
In Static Analysis (SAS), volume 4134 of LNCS, pages 144–160. Springer, 2006. DOI:
10.1007/11823230 10.
Denis Gopan and Thomas W. Reps. Lookahead widening. In Computer Aided Verification
(CAV), volume 4144 of LNCS, pages 452–466. Springer, 2006. DOI: 10.1007/11817963 41.
Denis Gopan and Thomas W. Reps. Low-level library analysis and summarization. In Computer Aided Verification (CAV), volume 4590 of LNCS, pages 68–81. Springer, 2007. DOI:
10.1007/978-3-540-73368-3 10.
25

Sumit
Gulwani,
Saurabh
Srivastava,
and
Ramarathnam
Venkatesan.
Program analysis as constraint solving.
In Programming Language Design and Implementation (PLDI). ACM, 2008. DOI: 10.1145/1375581.1375616.
Nicolas Halbwachs. Détermination automatique de relations linéaires vérifiées par les variables
d’un programme. PhD thesis, Université scientifique et médicale de Grenoble, 1979.
IEEE standard for Binary floating-point arithmetic for microprocessor systems. IEEE, 1985.
ANSI/IEEE Std 754-1985.
Bertrand Jeannet. Dynamic partitioning in linear relation analysis. application to the verification of reactive systems. Formal Methods in System Design, 23(1):5–37, July 2003.
Deepak Kapur. Automatically generating loop invariants using quantifier elimination.
ACA (Applications of Computer Algebra), 2004.

In

Akash Lal, Gogul Balakrishnan, and Thomas Reps. Extended weighted pushdown systems.
In CAV05, pages 343–357. DOI: 10.1007/11817963 32.

Francesco Logozzo and Manuel Fähndrich. On the relative completeness of bytecode analysis versus source code
In Compiler Construction (CC), volume 4959 of LNCS, pages 197–212. Springer, 2008.
DOI: 10.1007/978-3-540-78791-4 14.
Zohar Manna and John McCarthy. Properties of programs and partial function logic. In
Machine Intelligence, 5, pages 27–38. Edinburgh University Press, 1969.
Zohar Manna and Amir Pnueli. Formalization of properties of functional programs. J. ACM,
17(3):555–569, 1970. DOI: 10.1145/321592.321606.
Florian Martin, Martin Alt, Reinhard Wilhelm, and Christian Ferdinand. Analysis of loops.
In Compiler Construction (CC), volume 1383 of LNCS, pages 80–94. Springer, 1998.
Antoine Miné. The octagon abstract domain. In Reverse Engineering (WCRE), pages 310–
319. IEEE, 2001. DOI: 10.1109/WCRE.2001.957836.
Antoine Miné. Relational abstract domains for the detection of floating-point run-time errors.
In Programming Languages and Systems (ESOP), volume 2986 of LNCS, pages 3–17.
Springer, 2004.
David Monniaux. Compositional analysis of floating-point linear numerical filters. In CAV05,
pages 199–212. DOI: 10.1007/b138445.
David Monniaux. A quantifier elimination algorithm for linear real arithmetic. In LPAR
(Logic for Programming, Artificial Intelligence, and Reasoning), LNCS. Springer, 2008a.
David Monniaux. Optimal abstraction on real-valued programs. In Static analysis (SAS),
number 4634 in LNCS, pages 104–120. Springer, 2007. DOI: 10.1007/978-3-540-74061-2 7.
David Monniaux. The pitfalls of verifying floating-point computations. ACM Transactions
on programming languages and systems, 30(3):12, 2008b. DOI: 10.1145/1353445.1353446.

26

George C. Necula, Scott McPeak, Shree P. Rahul, and Westley Weimer.
CIL: Intermediate language and tools for analysis and transformation of C programs.
In Compiler Construction (CC), volume 2304 of LNCS, pages 209–265. Springer, 2002.
DOI: 10.1007/3-540-45937-5 16.

Radu Rugina and Martin Rinard. Symbolic bounds analysis for pointers, array indices, and accessed memory re
ACM Trans. on Programming Languages and Systems (TOPLAS), 27(2):185–235, 2005.
DOI: 10.1145/349299.349325.
Sriram
Sankaranarayanan,
Henny
Sipma,
and
Zohar
Manna.
Constraint-based linear-relations analysis. In SAS, number 3148 in LNCS, pages 53–
68. Springer, 2004.
Sriram
Sankaranarayanan,
Henny
Sipma,
and
Zohar
Scalable analysis of linear systems using mathematical programming.
In
pages 21–47. DOI: 10.1007/b105073.

Manna.
VMCAI05,

Helmut Seidl, Andrea Flexeder, and Michael Petter. Interprocedurally analysing linear inequality relations. In ESOP07, pages 284–299. DOI: 10.1007/978-3-540-71316-6 20.
Micha Sharir and Amir Pnueli. Two approaches to inter-procedural data-flow analysis. In
Program Flow Analysis: Theory and Application. Prentice-Hall, 1981.
VMCAI05. Verification, Model Checking and Abstract Interpretation (VMCAI), number 3385
in LNCS, 2005. Springer. DOI: 10.1007/b105073.

27

