University of Leeds SCHOOL OF COMPUTING RESEARCH REPORT SERIES Report 2005.06

A Linear Domain for Analyzing the Distribution of Numerical Values^{1}

by

Roberto Bagnara^{2} & Katy Dobson^{3} & Patricia M. Hill^{3} & Matthew Mundell^{3} & Enea Zaffanella^{2}

December 2005

1 This work has been partly supported by EPSRC project EP/C520726/1 “Numerical Domains for Software Analysis”, by MIUR project “AIDA — Abstract Interpretation: Design and Applications,” and by a Royal Society (UK) International Joint Project (ESEP) award.

2 Department of Mathematics, University of Parma, Italy

{bagnara,zaffanella}@cs.unipr.it

3 School of Computing, University of Leeds, UK

{katyd,hill,mattm}@comp.leeds.ac.uk

Abstract. This paper explores the abstract domain of grids, a domain that is able to represent sets of equally spaced points and hyperplanes over an n-dimensional vector space. Such a domain is useful for the static analysis of the patterns of distribution of the values program variables can take. Besides the bare abstract domain, we present a complete set of operations on grids that includes all that is necessary to define the abstract semantics and the widening operators required to compute it in a finite number of steps. The definition of the domain and its operations exploit well-known techniques from linear algebra as well as a dual representation that allows, among other things, for a concise and efficient implementation.

1 Introduction

The static analysis of numerical information about the values program variables can take is a challenging problem that has led to considerable research in both the theory and its practical realization. We distinguish between two kinds of numerical information: outer limits (or bounds within which the values must lie) and the pattern of distribution of these values. Both kinds of information have important applications. For example, in the field of automatic program verification, limit information is crucial to ensure that array accesses are within bounds, while distribution information is what is required to ensure that the external memory accesses obey the alignment restriction imposed by the host architecture. In the field of program optimization, limit information can be used to compile out various kinds of run-time tests, whereas distribution information enables several transformations for efficient parallel execution as well as optimizations that enhance cache behavior.

Both limit and distribution information often come in a relational form; for instance, the outer limits or the pattern of possible values of one variable may depend on the values of one or more other variables. Domains that can capture relational information are generally much more complex than domains that do not have this capability; in exchange they usually offer significantly more precision, often important for the overall performance of the client application. Relational limit information can be captured, among other possibilities, by means of polyhedral domains, that is, domains that represent regions of some n-dimensional vector space bounded by a finite set of hyperplanes [11]. While several polyhedral domains have been proposed and that of convex polyhedra, the most popular relational polyhedral domain, has been thoroughly researched and widely used, relational domains for representing the (linear) distribution of numerical values have been less well researched. For instance, previous investigations have hardly touched on the issue of widening operators; ignored the domain difference operation; and have failed to provide any support for a non-relational approximation. In addition, proposed algorithms are unnecessarily complex and standard mathematical procedures have not been utilized to their best advantage. Moreover, as far as we know and at the time of writing, there is no available implementation providing all the basic operations needed

-4 -2 4 2 -2 -4  6 4 2

(a) The grid L

-4 -2 4 2 -2  4 2

(b) The grid L′

Fig. 1. Two grids in ℝ^{2}

by a relational abstract domain for distribution information. This is in spite of the fact that previous research has shown that a knowledge about the (discrete) distribution of numerical information, especially when combined with that of the limit information, can significantly improve the quality of the analysis results [1]. The purpose of this paper is to fill this gap: we will present a complete relational domain of grids that captures numerical distribution information and briefly describe its implementation within the Parma Polyhedra Library (PPL) [5, 6].

1.1 Grids in a Nutshell

The diagrams in Figure 1 illustrate alternative approaches to the description of any grid; either by means of a finite set of congruence relations that all the grid points must satisfy (given by the dashed lines) or by means of a finite set of generating vectors used for constructing the grid points and lines (given by the filled squares and the arrow).

Consider first Figure 1(a). The set of points marked by squares at the intersection of the dashed lines represent the grid L. This denotes the distribution of possible values of integer variables x and y resulting from the execution of the following program fragment (based on an example in [11]) for any value of m:

    x := 2; y := 0;
    for i := 1 to m
      if ... then
        x := x + 4
      else
        x := x + 2; y := y + 1
      endif
    endfor

Note that the vertical dashed lines represent the set of points satisfying the congruence relation x = 0 (mod 2) while the sloping lines represent the set of points satisfying x + 2y = 2 (mod 4). The set

C = {x = 0 (mod 2), x + 2y = 2 (mod 4)}

is called a congruence system and said to describe L. The filled squares in Figure 1(a) represent the points

p_{1} = (\begin{matrix}2 \\ 0 \end{matrix}) , p_{2} = (\begin{matrix}6 \\ 0 \end{matrix}) and p_{3} = (\begin{matrix}4 \\ 1 \end{matrix})

while all the squares (both filled and unfilled) in the diagram mark position vectors v = π_{1}p_{1} + π_{2}p_{2} + π_{3}p_{3}, where π_{1}, π_{2}, π_{3} ∈ ℤ and π_{1} + π_{2} + π_{3} = 1. The set of points P = {p_{1}, p_{2}, p_{3}} is said to generate L. Some of these generating points can be replaced by parameters that give the direction and spacing for the neighbouring points. Specifically, by subtracting the point p_{1} from each of the other two generating points p_{2}, p_{3} we obtain

q_{2} = (\begin{matrix}4 \\ 0 \end{matrix}) and q_{3} = (\begin{matrix}2 \\ 1 \end{matrix}) .

It follows that each point v ∈ L can be written as v = p_{1} + π_{2}q_{2} + π_{3}q_{3} for some π_{2}, π_{3} ∈ ℤ; the set Q = {q_{2}, q_{3}} is called a parameter set for L.

Consider next the grid L′ illustrated by the dashed line in Figure 1(b) which is all the points that satisfy the equality x = y + 1. This denotes the distribution of values of the real variables x and y after an assignment x := y + 1, assuming that nothing is known about the value of y. In this case, as we regard equalities as congruences modulo 0, we say that the congruence system

C′ = {x − y = 1}

describes L′. Observe also that the grid L′ consists of all points that can be obtained as λℓ + p′, for any λ ∈ ℝ, where

ℓ = (\begin{matrix}1 \\ 1 \end{matrix}) and p′ = (\begin{matrix}1 \\ 0 \end{matrix});

the vector ℓ defines the direction of the line and the vector p′ its position (illustrated in Figure 1(b) by the arrow and the filled square, respectively).

From what we have just seen, any grid can be represented both by a congruence system and by a generator system. The latter may consist of three components: a set of lines, a set of parameters and a set of points. For instance, the triples G_{1} = (∅, ∅, P ) and G_{2} = (∅, Q, {p_{1}}) are both generator systems for L while the triple G′ = ({ℓ}, ∅, {p′}) is a generator system for L′.

1.2 Implementing Grids

The congruence and generator systems for representing elements of the domain of rational grids can be seen to have close parallels with the constraint and generator systems for describing convex polyhedra. In this paper, we develop and exploit this analogy. First, as for the systems representing convex polyhedra, by adding an extra dimension, the representations can be made homogeneous [27, 29]. Secondly, as a direct result of [35, Section 4.4], the homogenized congruence and generator systems are dual and, hence, equivalent representations for the domain of rational grids. Thus the congruence and generator systems for representing grids form the two components of a double description method for the grid domain that is very similar to that of the double description method for convex polyhedra [25].

Just as for the domain of convex polyhedra, for a double description method for grids to be viable, we require algorithms that can convert between a rational grid’s generator and congruence systems. In fact, for this domain, the standard techniques for matrix inversion can be used as the basis for the conversion algorithms. Moreover, since a congruence or generator system may contain some redundancies, or, even when there are no explicit redundancies, there may exists smaller a set of congruences or generators that generates or describes the same grid, we require algorithms for minimizing the representations. As before, there are standard matrix algorithms (for instance, the Hermite normal form algorithm [28, 35]) that can be used as a basis for minimization algorithms. Note that the minimal representations are needed not just to economize on memory space and subsequent computation time, but also for the conversion algorithms mentioned above.

1.3 Related Work

As far as we know, Granger provided the first account of an analysis for congruence information; in [14] it is shown how a simple integral and non-relational grid domain could be used in a static analyzer. In this domain, an element is characterized as the set of all integers that are congruent to c modulo m where c and m are integers; that is, the set of all integers x where, for some k ∈ ℤ, x = c + km; which does, in fact, immediately provide an equivalent one-dimensional generator system consisting of a point (c) and parameter (m). One particularly useful aspect of this work is that it indicates some practical applications of the domain that carry over to the more complex domain considered here. For instance, it is shown how an analysis using the grid domain can obtain more precise information than an analysis using a constant propagation domain or a sign domain. It is then indicated how the results may be applied for automatic vectorization. Independently, Larsen et al. [22] have also developed a static analyzer over a non-relational grid domain; this analyzer is specifically designed to detect when the dynamic memory addresses are congruent with respect to a given modulus (the actual modulus being dependent on the application). It is shown how this information may be used to construct a comprehensive set of program transformations that can lead to energy savings on low-power architectures and performance improvements on multimedia processors.

Note that, prior to the current stream of work on congruence relations, Karr [21] had pioneered the consideration of a linear relational equality domain for program analysis. More recently, Granger, building on his previous work for non-relational congruences [14], generalized Karr’s work to linear congruence relations [16] (see also [15]). This paper focuses on an integral domain, but as indicated in the subsequent paper by the same author [17], this work easily generalizes to the domain of rational grids. Two representations for the domain are described, sets of congruences and sets of generators. In [16], it is shown how standard matrix algorithms [36] can be used to convert from a generator system to a congruence system; a more complex algorithm is also provided for conversion from a congruence system to a generator system. Operators for comparing grids and computing the greatest lower and least upper bounds are also described in [16].

Most previous work (such as that in [16, 17]) as well as the work in this paper assumes fully relational congruence systems for representing grids. In contrast, by applying a generic technique, Miné [24] shows how to construct a weakly relational zone-congruence domain (that is, a domain that only allows congruences that have the form x − y = a (mod b) where a and b are rationals) from the non-relational congruence domain described in [14].

Müller-Olm and Seidl [27] have recently improved on the work of Granger [16, 17]; one aspect of the improvement is obtained by recognizing, as we do here, that the inhomogeneous representations may be transformed to homogeneous ones by adding an extra dimension. However, we note that [27] is focused, not on the definition of the domain itself, but on the definition of an abstract interpretation of a constraint system that can characterize the concrete program semantics, particularly in the case when there are procedural calls. Thus the main contributions of [27] are largely orthogonal to the contributions here. Note also that as stated in [27], the framework described there subsumes analysis described in previous work by the same authors; in particular, it subsumes the analysis described in [26] designed specifically for analyzing certain integer arithmetic computations that are modulo 2^{w} where, for example, w is 32 or 64.

Following a completely independent stream of research and publications, Ancourt [1] considered a domain of grids (i.e., lattices) combined with the domain of convex polyhedra; that is the domain of Z-polyhedra defined as “the intersections of polyhedra and integral lattices”. We are primarily interested here in the “integral lattices” component which may be seen as a subdomain of the domain of rational grids where the grid points are all integral vectors and is fulldimensional. In all the papers concerning the Z-polyhedral domain [1, 29, 31, 32], the representation of these integral lattices is similar to that of the generator representation given in Section 3. All the operations on Z-polyhedra (and therefore the lattices) require canonic representations; hence Quinton et al. [31, 32] defines a canonical form for these lattices with a method for its computation. We note that the algorithm for computing the canonic form has complexity O(n^{4}), where n is the number of dimensions of the vector space. Other operations provided for lattices are those of intersection, affine image and affine preimage. As there is no congruence representation, the intersection of two lattices is computed directly from the generator representations [1]; a refined version of this method is provided in [31] which we note that, as for computing the canonic form, it has complexity O(n^{4}). On the other hand, the operations of grid join ‘⊕’ and grid difference ‘⊖’ (as defined in Section 4) are not considered. Instead the union operator takes two lattices L_{1} and L_{2} and returns the set {L_{1}, L_{2}} of the two unless one (say L_{1}) is contained in the other, in which case they will return the larger L_{2}. Moreover, the difference operation is also exact and, in general, returns a set of lattices.

The domain of integral lattices and operations on them has been implemented in PolyLib[23] following the methodology described in [31, 32]. Thus several of the main operations on lattices, such as union, return sets of lattices while others are concerned with manipulating and simplifying these sets. We note in particular that, as explained by Nookala and Risset [29], all the operations are implemented using a homogeneous representation. This is obtained, as we do here, by adding an extra dimension to the inhomogeneous representation that is input to the system via the user interface.

1.4 Contributions

The paper provides a full account of the relational domain of grids including its alternative representations, operations required for static analysis, and key algorithms needed for its implementation. Although some of this is already in the literature, it is widely dispersed and has never before been assembled into a single description. There are also several specific contributions.

– Although the Gaussian elimination and Hermite normal form algorithms form the basis of the proposed conversion and minimization algorithms for grids, we have had to extend them so as to apply to any grid and its representations as well as adapt them so as to be suitable for a real, efficient implementation. Therefore we outline the actual algorithms we use in our implementation; some of which are shown to have complexities strictly better than that of previous proposals [16, 27].

– By taking the union of the generator systems (resp., congruence systems) representing two grids we can compute their least upper bound (resp., the greatest lower bound); therefore, since we have reduced the complexity of the algorithms for minimizing and converting between representations, our implementations for these operations also have smaller complexity than that of previous proposals [16].

– We provide an algorithm for computing the smallest grid containing the difference of two grids. Observe that there are no previous proposals for this operator. The only papers that considers the computation of the difference of two grids are those concerning the Z-polyhedral domain [1, 29, 31, 32]; here, as mentioned above, the difference operation returns the exact difference of the grids which has to be represented by a set of grids.

– One major difference between the domains of rational and integral grids is that the rational grid domain 𝔾_{n} does not satisfy the ascending chain condition. In [17], a widening is proposed for non-relational rational grids and it is briefly stated that it might be generalized to the case of rational grids but without any further details; how this might be done is unclear. Here we provide a practical widening operator for rational grids that, when restricted to the non-relational grid, does correspond to that proposed in [17].

– In order to dynamically manage the complexity/precision tradeoff of program analysis, the analyzer should be able to dynamically switch (in response to timeouts or driven by heuristics) to a simpler analysis domain. We provide a non-relational approximation for an arbitrary grid that can be implemented by means of an interval domain together with algorithms for conversion from a grid to the approximation and vice-versa.

1.5 Plan of the Paper

The rest of the paper is structured as follows. The required notations and preliminary concepts are given in Section 2. Section 3 introduces a grid together with its congruence and generator representations; Section 4 describes the basic operators on grids needed for program analysis; Section 5 provides the main algorithms needed to support the double description; Section 6 introduces the grid widening operation and its implementation; Section 7 gives some examples of grids applications. We conclude the main body of the paper in Section 8. Appendix A contains the proofs of all the stated results.

2 Preliminaries

The cardinality of a set S is denoted by # S. The set of integers is denoted by ℤ, rationals by ℚ and reals by ℝ. For any a, b ∈ ℝ where a ≠ 0, we say a divides b, denoted by a|b, if, for some m ∈ ℤ, am = b. The complexities we give for the different algorithms assume a unit cost for every arithmetic operation; we take the computation of the greatest common divisor of a pair of numbers a, b ∈ ℝ to be a single operation. Given sets X, Y and any relation R ⊆ X × Y, the image for R on a subset A of X is { y ∈ Y | ∃x ∈ A . (x, y) ∈ R }, and the preimage for R on a subset B of Y is { x ∈ X | ∃y ∈ B . (x, y) ∈ R }.

2.1 Vectors and Matrices

For each i ∈ {1, …, n}, v_{i} denotes the i-th component of the (column) vector v ∈ ℝ^{n}. The empty vector (in ℝ^{0}) is denoted by ϵ; a vector that has all its elements equal to zero is denoted by 0; and, for 1 ≤ i ≤ n, the notation e_{i} denotes the vector in ℝ^{n} with 1 in the i-th position and zeroes in every other position. For v ∈ ℝ^{n}, piv_{<}(v) denotes the maximum index i such that v_{i} ≠ 0; if v = 0, we define piv_{<}(v) := 0. Similarly, piv_{>}(v) denotes the minimum index i such that v_{i} ≠ 0; if v = 0, we define piv_{>}(v) := n+1. Any vector v ∈ ℝ^{n} is also a matrix in ℝ^{n×1} so that it can be manipulated with the usual matrix operations of addition and multiplication, both by a scalar and by another matrix. On the other hand, it is often convenient to consider a matrix H = (h_{1}, …, h_{m}) ∈ ℝ^{n×m} as a finite set of vectors {h_{1}, …, h_{m}} ⊆ ℝ^{n}. For each i ∈ {1, …, n} and j ∈ {1, …, m}, the ij-th component of a matrix H ∈ ℝ^{n×m} is denoted by H_{ij} and the i-th row by H_{i}. The transposition of a matrix H is denoted by H^{T}. The scalar product of v, w ∈ ℝ^{n}, denoted ⟨v, w⟩, is the real number v^{T}w = ∑_{i=1}^{n} v_{i}w_{i}.

Vectors v_{1}, …, v_{m} ∈ ℝ^{n} are said to be affinely independent if, for all λ ∈ ℝ^{m}, λ = 0 is the only solution of the set of equations {∑_{i=1}^{m} λ_{i}v_{i} = 0, ∑_{i=1}^{m} λ_{i} = 0}.

2.2 Integer Combinations

Let S = {v_{1}, …, v_{k}} ⊆ ℝ^{n} be a set of k vectors. For all scalars λ_{1}, …, λ_{k} ∈ ℝ, the vector v = ∑_{j=1}^{k} λ_{j}v_{j} is said to be a linear combination of the vectors in S. Such a combination is said to be

– an affine combination, if ∑_{j=1}^{k} λ_{j} = 1;

– an integral combination, if λ_{1}, …, λ_{k} ∈ ℤ;

– an integral affine combination, if it is both integral and affine.

We denote by affine.hull(S) (resp., int.hull(S), int.affine.hull(S)) the set of all the affine (resp., integer, integer affine) combinations of the vectors in S.

2.3 Congruences and Congruence Relations

For any a, b, f ∈ ℝ, a ≡_{f} b denotes the congruence ∃μ ∈ ℤ . a − b = μf. In the case that f = 0, the congruence denotes the equality a = b. Let 𝕊 be either ℚ or ℝ. For each vector a ∈ 𝕊^{n} and scalars b, f ∈ 𝕊, the notation ⟨a, x⟩ ≡_{f} b stands for the linear congruence relation in 𝕊^{n} defined by the set of vectors { v ∈ ℝ^{n} | ∃μ ∈ ℤ . ⟨a, v⟩ = b + μf }; f is called the frequency and b the base value of the relation; when f = 0, the congruence relation denotes the equality ⟨a, x⟩ = b; when f ≠ 0, the congruence relation is said to be proper. Thus, provided a ≠ 0, the congruence relation ⟨a, x⟩ ≡_{f} b defines the set of affine hyperplanes { (⟨a, x⟩ = b+μf) | μ ∈ ℤ }. The congruence ⟨0, x⟩ ≡_{f} b defines the universe ℝ^{n} if b ≡_{f} 0, and the emptyset, otherwise. We will assume that in such a congruence (when a = 0) we have b ≠ 0. Any vector that satisfies one of the equalities ⟨a, x⟩ = b + μf for any μ ∈ ℤ is said to satisfy the congruence relation ⟨a, x⟩ ≡_{f} b. We do not distinguish between syntactically different congruences defining the same set of vectors in 𝕊^{n} so that, e.g., x ≡_{1} 2 and 2x ≡_{2} 4 are considered to be the same congruence.

We extend the pivot notation to congruences as follows. If β = (⟨a, x⟩ ≡_{f} a_{0}) and piv_{<}(a) = k, then piv_{<}(β) := k. We say that congruence β is pivot equivalent to γ = (⟨c, x⟩ ≡_{g} c_{0}), and write β ⇑ γ, if piv_{<}(γ) = piv_{<}(β) = k and ga_{k} = fc_{k}. Observe that, as a_{k}, c_{k} ≠ 0, this means that β and γ are either both equalities or both proper congruences.

3 Rational Grids

In this section we introduce rational grids and their representation. Note that previously the terminology ‘lattice’ instead of ‘grid’ has been used for elements of domains similar to the grids domain described here. However, to avoid any confusion with the use of ‘lattice’ in its set-theoretic context (which is particularly relevant when working in the framework of abstract interpretation), we prefer the terminology ‘grid’ for sets of vectors that can represent discrete and linearly repetitive information. Moreover, as G is used to denote a generator system, we use L to denote a rational grid.

3.1 The Grid Domain and the Congruence Representation

A congruence system in ℚ^{n} is a finite set of congruence relations C in ℚ^{n}. As we do not distinguish between syntactically different congruences defining the same set of vectors, we can assume that all proper congruences in C have modulus 1.

Definition 1. The set of vectors L is a rational grid in ℝ^{n} described by a congruence system C in ℚ^{n} if and only if L is the set of points in ℝ^{n} that satisfy all the congruences in C. We also say that C is a congruence system for L and write L = gcon(C).

If gcon(C) = ∅, then we say that C is inconsistent. For example, the congruence systems {⟨0, x⟩ ≡_{0} 1} and {⟨a, x⟩ ≡_{2} 0, ⟨a, x⟩ ≡_{2} 1}, for any a ∈ ℚ^{n}, both describe the empty grid in ℝ^{n}. In fact, the first congruence system requires that 0 = 1, while the second one requires that the value of an expression is both even and odd, so that they are both inconsistent.

The grid domain 𝔾_{n} is the set of all rational grids in ℝ^{n}. When ordering grids by the set inclusion relation, the empty set ∅ and the vector space ℝ^{n} (which is described by the empty set of congruence relations) are, respectively, the smallest and the biggest elements of 𝔾_{n}. The vector space ℝ^{n} is also called the universe grid. In set theoretical terms, 𝔾_{n} is a lattice under set inclusion.

The space dimension of a grid L ∈ 𝔾_{n} is the dimension n ∈ ℕ of the corresponding vector space ℝ^{n}. If the maximum number of affinely independent points in L is k + 1, then dim(L) = k denotes the affine dimension of L. The affine dimension of an empty grid is defined to be 0. Thus we have 0 ≤ dim(L) ≤ n.

Let C be a congruence system and L = gcon(C). Suppose also that the congruence relation β = (⟨a, x⟩ ≡_{f} b) is such that L_{β} = gcon({β}). We say that

– L is disjoint from β if L ∩ L_{β} = ∅; that is, adding β to C gives us the empty grid.

– L strictly intersects β if L ∩ L_{β} ≠ ∅ and L ∩ L_{β} ⊂ L; that is, adding β to C gives us a non-empty grid strictly smaller than L.

– L is included in β if L ⊆ L_{β}; that is, adding β to C leaves L unchanged.

Many algorithms used in the implementation require the congruence systems to have a minimal number of elements.

Definition 2. Suppose C is a congruence system in ℚ^{n}. Then we say that C is in minimal form if either C = {⟨0, x⟩ ≡_{0} 1} or C is consistent and, for each congruence β = (⟨a, x⟩ ≡_{f} b) ∈ C, the following hold:

1. if piv_{<}(β) = k, then k > 0 and a_{k} > 0;

2. for all β′ ∈ C \ {β}, piv_{<}(β′) ≠ piv_{<}(β).

Proposition 1. Let C be a congruence system in ℚ^{n}. Then there exists an algorithm for finding a congruence system C′ in minimal form such that gcon(C) = gcon(C′). Letting m = # C, the complexity of the algorithm is O(n^{2}m).

The pivot equivalence relation between congruences can be extended to congruence systems. That is, we say that congruence systems C_{1} and C_{2} are pivot equivalent if, for each β ∈ C_{1}, there exists γ ∈ C_{2} such that β ⇑ γ and, for each γ ∈ C_{2}, there exists β ∈ C_{1} such that γ ⇑ β. This equivalence relation provides an easy syntactic check to establish if two grids are equal, given that it is known that one is a subset of the other.

Proposition 2. Let L_{1} = gcon(C_{1}) and L_{2} = gcon(C_{2}) be non-empty grids in 𝔾_{n} such that L_{1} ⊆ L_{2}. If C_{1} and C_{2} are pivot equivalent congruence systems in minimal form, then L_{1} = L_{2}.

3.2 The Generator Representation

Let L be a grid in 𝔾_{n}. Then

– a vector p ∈ L is called a point of L;

– a vector q ∈ ℝ^{n}, where q ≠ 0, is called a parameter of L if L ≠ ∅ and p + μq ∈ L, for all points p ∈ L and all μ ∈ ℤ;

– a vector l ∈ ℝ^{n}, where l ≠ 0, is called a line of L if L ≠ ∅ and p + λl ∈ L, for all points p ∈ L and all λ ∈ ℝ.

We can generate any rational grid in 𝔾_{n} from a finite subset of its points, parameters and lines; each point in a grid is obtained by adding a linear combination of its generating lines to an integral combination of its parameters and an integral affine combination of its generating points.

If L, Q and P are each finite subsets of ℚ^{n} and

L := linear.hull(L) + int.hull(Q) + int.affine.hull(P)  (1)

where the symbol ‘+’ denotes the Minkowski’s sum,^{2} then L ∈ 𝔾_{n} is a rational grid [15, Theorem 48]. The 3-tuple (L, Q, P) is said to be a generator system for L and we write L = ggen(L, Q, P). L is said to subsume a generator g if adding g to any generator system representing L does not change L.

Note that the grid L = ggen(L, Q, P) = ∅ if and only if the set of points P = ∅. If P ≠ ∅, then L = ggen(L, ∅, Q_{p} ∪ P) where, for some p ∈ P, Q_{p} = { p + q ∈ ℝ^{n} | q ∈ Q }.

As for congruence systems, for many procedures in the implementation, it is useful if the generator systems have a minimal number of elements.

2 This is defined, for each S, T ⊆ ℝ^{n}, by S + T := { s + t ∈ ℝ^{n} | s ∈ S, t ∈ T }.

Definition 3. Suppose G = (L, Q, P) is a generator system in ℚ^{n}. Then we say that G is in minimal form if either L = Q = P = ∅ or # P = 1 and, for each generator v ∈ L ∪ Q, the following hold:

1. if piv_{>}(v) = k, then v_{k} > 0;

2. for all v′ ∈ (L ∪ Q) \ {v}, piv_{>}(v′) ≠ piv_{>}(v).

Proposition 3. Let G = (L, Q, P) be a generator system in ℚ^{n}. Then there exists an algorithm for finding a generator system G′ in minimal form such that ggen(G′) = ggen(G). Letting m = # L + # Q, the complexity of the algorithm is O(n^{2}m).

3.3 Double Description

We have shown that any grid L can be described by using a congruence system C and also generated by a generator system G. For the same reasons as for the polyhedral domain, it is useful to represent the grid L by the double description (C, G). In order to maintain and exploit such a view of a grid, we will use a double description method for grids which is a collection of results very similar to those for convex polyhedra showing that, given one kind of representation, there are algorithms for computing a representation of the other kind and for minimizing both representations. As for convex polyhedra, having easy access to both the representations is assumed in the implementation of many operators such as, for example, intersection and grid join described in Section 4.

Suppose we have a double description (C, G) of a grid L ∈ 𝔾_{n}, where both C and G are in minimal form. Then, it follows from the definition of minimal form that # C ≤ n + 1 and # L + # Q ≤ n. In fact, we have a stronger result.

Proposition 4. Let (C, G) be a double description where both C and G are in minimal form. Letting C = E ∪F where E and F are sets of equalities and proper congruences, respectively, and G = (L, Q, P), then

# F = # Q = n − # L − # E.

Example 1. Consider the grids L and L′ in Figure 1. The congruence systems C and C′ are in minimal form and the generator systems G_{2} and G′ are also in minimal form; however, G_{1} is not in minimal form as it contains more than one point. Furthermore, for i = 1, 2, the pairs (C, G_{i}) are double descriptions for L while (C′, G′) is a double description for L′.

4 Operations on Rational Grids

In this section we describe some operations on rational grids that are analogous to ones that are already provided for the domain of convex polyhedra and used in program analyzers.

\begin{aligned}8 & 10 & 12 & 4 & 2 & 6 & 12 & 10 & 8 & 6 & 4 & 2 & L, & L′, & L ∩ L′ \\ & & (a) Grid L_{1} ∩ L_{2} \end{aligned}  \begin{aligned}8 & 10 & 12 & 4 & 2 & 6 & 12 & 10 & 8 & 6 & 4 & 2 & L, & L′, & L ⊕ L′ \\ & & (b) Grid L_{1} ⊕ L_{2} \end{aligned}

Fig. 2. Grid intersection and join

4.1 Grid Intersection and Grid Join

For any pair of grids L_{1}, L_{2} ∈ 𝔾_{n}, the intersection of L_{1} and L_{2}, defined as the set intersection L_{1} ∩L_{2}, is the largest grid included in both L_{1} and L_{2}; similarly, the grid join of L_{1} and L_{2}, denoted by L_{1} ⊕ L_{2}, is the smallest grid that includes both L_{1} and L_{2}.

In theoretical terms, the intersection and grid join operators defined above are the binary meet and the binary join operators on the lattice 𝔾_{n}.

Example 2. Consider the 𝔾_{2} grids L = gcon(C) and L′ = gcon(C′) where

C := {x ≡_{2} 0, −x + y ≡_{3} 0} and C′ := {x ≡_{4} 0, −x + 2y ≡_{6} 0}.

The grids L and L′ are illustrated by the squares and circles in Figure 2(a), respectively. The grid-intersection L ∩ L′ = gcon(C ∪ C′) is illustrated by the filled circles in the same diagram. Observe that, if

C′′ := {x ≡_{12} 0, y ≡_{3} 0}

then C′′ is a reduced form of C ∪ C′ so that we also have L ∩ L′ = gcon(C′′).

Example 3. Consider the grids L = ggen((∅, ∅, P)) and L′ = ggen((∅, ∅, P′)) in 𝔾_{2} such that

P := (\begin{matrix}2 & 0 & 0 \\ 2 & 3 & 0 \end{matrix}) and P′ := (\begin{matrix}4 & 0 & 0 \\ 2 & 3 & 0 \end{matrix}) .

The grids L and L′ are illustrated by the squares and filled circles in Figure 2(b), respectively. The grid-join L ⊕ L′ = ggen((∅, ∅, P ∪ P′)) is illustrated by all (open and filled) circles in the same diagram. Observe that, if

Q′′ := (\begin{matrix}2 & 0 \\ 0 & 1 \end{matrix}) and P′′ := (\begin{matrix}0 \\ 0 \end{matrix}) ,

\begin{aligned}4 & 5 & 6 & 2 & 1 & 3 & 6 & 5 & 4 & 3 & 2 & 1 & L_{1}, & L_{2}, & L_{1} ⊖ L_{2} \\ & & (a) Grid L_{1} ⊖ L_{2} \end{aligned}  \begin{aligned}4 & 5 & 6 & 2 & 1 & 3 & 6 & 5 & 4 & 3 & 2 & 1 & L_{1}, & L_{3}, & L_{1} ⊖ L_{3} \\ & & (b) Grid L_{1} ⊖ L_{3} \end{aligned}

Fig. 3. Grid difference

then (∅, Q′′, P′′) is minimal form of (∅, ∅, P ∪ P′) so that we also have L ⊕ L′ = ggen((∅, Q′′, P′′)). Note that, in this example, L ⊕ L′ ≠ L ∪ L′.

4.2 Grid Difference

For any pair of grids L_{1}, L_{2} ∈ 𝔾_{n}, the grid difference of L_{1} and L_{2}, denoted by L_{1} ⊖ L_{2}, is defined as the smallest grid containing the set-theoretic difference of L_{1} and L_{2}.

Example 4. Consider grids

L_{1} := gcon({x ≡_{1} 0, y ≡_{1} 0})
L_{2} := gcon({x ≡_{1} 0, x + y ≡_{2} 0})
L_{3} := gcon({x ≡_{1} 0, x + y ≡_{4} 0}).

The grids L_{1} and L_{2} are illustrated by all the circles (open and solid) and squares, respectively, in Figure 3(a); the grid difference

L_{1} ⊖ L_{2} = gcon({x ≡_{1} 0, x + y ≡_{2} 1}).

is illustrated by the filled circles. The grids L_{1} and L_{3} are illustrated by the circles and squares, respectively, in Figure 3(b); in this case, the grid difference is L_{1} ⊖ L_{3} = L_{1}.

Algorithm 1: The grid difference algorithm

Input: Two grids L_{1} = gcon(C_{1}) and L_{2} = gcon(C_{1}) in 𝔾_{n}.

Output: A grid in 𝔾_{n}.

(1) if L_{1} = ∅ ∨ L_{2} = ∅

(2) return L_{1}

(3) if L_{1} ⊆ L_{2}

(4) return ∅

(5) L′ := the empty grid in 𝔾_{n}

(6) while ∃β = (e ≡_{f} 0) ∈ C_{2}

(7) C_{2} := C_{2} \ {β}

(8) if L_{1} ⊈ gcon({β})

(9) if L_{1} ⊆ gcon({2e ≡_{f} 0})

(10) L_{β} := gcon(C_{1} ∪ {2e − f ≡_{2f} 0})

(11) L′ := L′ ⊕ L_{β}

(12) else

(13) return L_{1}

(14) return L′

Algorithm 1 provides an implementation for grid difference.

Proposition 5. Let L_{1}, L_{2} ∈ 𝔾_{n} and suppose that L is the grid returned by Algorithm 1. Then L = L_{1} ⊖ L_{2}.

4.3 Affine Images and Preimages

Given a grid L ∈ 𝔾_{n}, a variable x_{k} and linear expression e = ⟨a, x⟩ + b with coefficients in ℚ, the affine image operator ϕ(L, x_{k}, e) maps the grid L to

{ (p_{1}, …, p_{k−1}, ⟨a, p⟩ + b, p_{k+1}, …, p_{n})^{T} ∈ ℝ^{n} | p ∈ L }.

Conversely, the affine preimage operator ϕ^{−1}(L, x_{k}, e) maps the grid L to

{ p ∈ ℝ^{n} | (p_{1}, …, p_{k−1}, ⟨a, p⟩ + b, p_{k+1}, …, p_{n})^{T} ∈ L }.

Observe that the affine image ϕ(L, x_{k}, e) and preimage ϕ^{−1}(L, x_{k}, e) are invertible if and only if the coefficient a_{k} in the vector a is non-zero.

Example 5. Suppose L = ggen((∅, ∅, P)) ∈ 𝔾_{2}, where P = ( \begin{matrix}0 & 0 & 3 \\ 0 & 3 & 0 \end{matrix} ). Then

L_{1} = ϕ(L, x, 3x + 2y + 1) = ggen((∅, ∅, P_{1})) = gcon(C_{1})

where P_{1} = ( \begin{matrix}1 & 7 & 10 \\ 0 & 3 & 0 \end{matrix} ) and C_{1} = \{(x ≡_{3} 1), (x + y ≡_{9} 1)\}. On the other hand,

L_{2} = ϕ(L, x, y) = ggen((∅, ∅, P_{2})) = gcon(C_{2})

where P_{2} = ( \begin{matrix}0 & 3 \\ 0 & 3 \end{matrix} ) and C_{2} = \{(x ≡_{3} 0), (y ≡_{3} 0), (x − y = 0)\}; that is, the grid containing all the points whose coordinates are integral multiples of 3 and lie on line x = y.

Note that the affine preimage of L_{1} using the same variable x and linear expression 3x + 2y + 1 is the original grid L. Moreover, the affine preimage of L_{2} using variable x and linear expression y is the grid ggen((L_{3}, Q_{3}, P_{3})) where

L_{3} = (\begin{matrix}0 \\ 1 \end{matrix}) , Q_{3} = (\begin{matrix}0 \\ 3 \end{matrix}) and P_{3} = (\begin{matrix}0 \\ 0 \end{matrix}) ,

that is, all the points in ℝ^{2} where the y coordinate is an integral multiple of 3.

The generalized affine image and generalized affine preimage are extensions of the affine image and affine preimage operators defined above. Given a grid L ∈ 𝔾_{n}, linear expressions e′ = ⟨a′, x⟩ + b′, e = ⟨a, x⟩ + b with coefficients in ℚ and f ∈ ℚ, the generalized affine image operator ψ(L, e′, f, e) will transform the grid L to the grid

\{ w ∈ ℝ^{n} | \begin{matrix}v ∈ L \\ (i ∈ \{1, …, n\} ∧ a′_{i} = 0) ⟹ w_{i} = v_{i} \\ ⟨a′, w⟩ + b′ ≡_{f} ⟨a, v⟩ + b \end{matrix} \}.

Note that, when e′ = x_{k} and f = 0, then the transformation is equivalent to the standard affine transformation on L with respect to the variable x_{k} and the affine expression e; that is

ψ(L, x_{k}, 0, e) = ϕ(L, x_{k}, e).

Example 6. Let L and L_{1} be as defined in Example 5. Then

ψ(L, x, 0, 3x + 2y + 1) = ϕ(L, x, 3x + 2y + 1) = L_{1}.

On the other hand,

ψ(L, x, 3, 3x + 2y + 1) = ggen((∅, ∅, P′_{1})) = gcon(C′_{1})

where P′_{1} = ( \begin{matrix}1 & 4 & 1 \\ 0 & 0 & 3 \end{matrix} ) and C′_{1} = \{(x ≡_{3} 1), (y ≡_{3} 0)\}.

4.4 Rectilinear Grids and Covering Boxes

We show how we can reuse the standard interval domain [7] to represent a rectilinear grid, which is a grid that can be defined by a set of non-relational congruences (i.e., congruences where at most one coefficient is non-zero). We also show how such grids can provide safe approximations for any rational grid.

A non-empty n-dimensional rational box B is a sequence (I_{1}, ….I_{n}) of rational intervals. A rational interval is a pair (μ_{i}, ν_{i}) ∈ (ℚ ∪ {∞})^{2} of bounds, called lower and upper, respectively. The lower bound is always less than or equal to the upper bound, where q < ∞ for each q ∈ ℚ. If both the bounds are in ℚ, the interval is said to be bounded.

4 5 6 2 1 3  6 5 4 3 2 1  L_{1} (a) Grid L_{1} 4 5 6 2 1 3  6 5 4 3 2 1  L_{2}, L′_{2}

(b) Grids L_{2} and L′_{2}

Fig. 4. Covering boxes for a grid

The essential idea we use here exploits the fact that the shape and size of an n-dimensional box B defines a rectangular “tile” that may be used to cover an n-dimensional vector space. It therefore follows that the box B determines a covering (i.e., tiling) of the n-dimensional vector space, where the given box B provides the position for one of the tiles. By defining a grid to be the vertices of the tiles in such a tiling, we obtain a rational rectilinear grid L and call B a covering box for L.

Let B = (I_{1}, ….I_{n}) be a non-empty box. For each i = 1, …, n, let I_{i} = (μ_{i}, ν_{i}); then, if μ_{i} ≠ ν_{i}, let

β_{i} := \begin{cases}(⟨e_{i}, x⟩ ≡_{ν_{i}−μ_{i}} μ_{i}), & if I_{i} is bounded; \\ (⟨e_{i}, x⟩ = μ_{i}), & if I_{i} is not bounded. \end{cases}

Then we say that the box B represents the grid L := gcon(C), where

C := { β_{i} | 1 ≤ i ≤ n, μ_{i} ≠ ν_{i} }.  (2)

Note that the congruence system C is in minimal form. Observe also that, when μ_{i} = ν_{i} for some 1 ≤ i ≤ n, there is no corresponding congruence in C for (μ_{i}, ν_{i}); this is because, in this case, the tiling will cover every value in this dimension and hence there will be a line e_{i} in the generator representation of L.

Consider now any non-empty rational grid L. A covering box for L is a rational box representing the smallest rectilinear grid that contains L.

Example 7. Consider grid L_{1} = gcon({x ≡_{3} 0, y ≡_{2} 1}) illustrated by all the circles in Figure 4(a); then L_{1} is rectilinear and and box B_{1} = {(0, 3), (1, 3)} is a covering box representing L_{1}. Consider the grid L_{2} = gcon({x ≡_{1} 0, x+y ≡_{3} 2}) illustrated by all the filled circles in Figure 4(b); then L_{2} is not rectilinear and box B_{2} = {(1, 2), (1, 2)}, is a covering box for L_{2}. Thus B_{2} represents the grid L′_{2} = gcon({x ≡_{1} 1, y ≡_{1} 1}) which is illustrated by all the circles in Figure 4(b).

We now provide a procedure for computing the covering box of a grid.

Proposition 6. Let L = ggen(G) where G = (L, Q, {p}). Let q = # Q and Q = {q_{1}, …, q_{q}}. Let B = (I_{1}, …, I_{n}) such that, for each i ∈ {1, …, n}:

1. if, for some ℓ ∈ L, ℓ_{i} ≠ 0, then let I_{i} := (0, 0);

2. if, for all ℓ ∈ L, ℓ_{i} = 0, and q_{1i} = ⋯ = q_{qi} = 0, let I_{i} := (p_{i}, ∞);

3. otherwise, let I_{i} := (p_{i}, p_{i} + |g|) where g = gcd({q_{1i}, …, q_{qi}}).

Then B is a covering box for L.

The complexity of computing the covering box for a grid in 𝔾_{n} using the procedure given in Proposition 6, is O(n^{2}).

It should be stressed that a grid L ∈ 𝔾_{n} does not have, in general, a unique covering box. For instance, if B = (I_{1}, …, I_{i}, …, I_{n}) is a covering box for L, and the interval I_{i} = (μ_{i}, ν_{i}) is bounded, then the box B′ = (I_{1}, …, I′_{i}, …, I_{n}) is also a covering box for L if, for some m ∈ ℤ, I′_{i} = (μ_{i} + m(ν_{i} − μ_{i}), ν_{i} + m(ν_{i} − μ_{i})).

4.5 Other Grid Operators

There are many operations that a practical domain of grids must provide for applications in program analysis and verification. For instance, the concatenation of two grids L_{1} ∈ 𝔾_{n} and L_{2} ∈ 𝔾_{m} (taken in this order) is the grid in 𝔾_{n+m} defined as

\{ (x_{1}, …, x_{n}, y_{1}, …, y_{m})^{T} ∈ ℝ^{n+m} | \begin{matrix}(x_{1}, …, x_{n})^{T} ∈ L_{1} \\ (y_{1}, …, y_{m})^{T} ∈ L_{2} \end{matrix} \}.

Other operators are required that add, remove, rename and map the space dimensions, or expand and fold them along the lines of [13]. Even though all these operations have been specified and implemented, space considerations make their discussion here inadvisable: the interested reader is referred to the web site of the Parma Polyhedra Library, http://www.cs.unipr.it/ppl/, where all the code and documentation is publicly available.

5 Implementation

In this section, we first describe convenient internal representations of the congruence and generator systems in terms of arrays (i.e., matrices). Then we show how matrix inversion provides a basis for converting between the (possibly singular) matrix representations of the congruence and generator systems in minimal form.

5.1 Homogeneous Representations

A congruence system C is homogeneous if, for all (⟨a, x⟩ ≡_{f} b) ∈ C, we have b = 0. Similarly, a generator system (L, Q, P) is homogeneous if 0 ∈ P. For the implementation, it is convenient to work with a homogeneous system. Thus we first convert any congruence or generator systems in ℚ^{n} to a homogeneous representation in ℚ^{n+1}. The extra dimension is denoted with a 0 subscript and x̂ denotes the vector (x_{0}, …, x_{n})^{T}; we also let e_{0} denote the vector (1, 0^{T})^{T}.

Consider the congruence system C = F ∪E in ℚ^{n} where E is a set of equalities and F is a set of proper congruences. Then the homogeneous form for C is the congruence system Ĉ = F̂ ∪ Ê in ℚ^{n+1} defined by:

F̂ := { ⟨f^{−1}(−b, a^{T})^{T}, x̂⟩ ≡_{1} 0 | (⟨a, x⟩ ≡_{f} b) ∈ F } ∪ {⟨e_{0}, x̂⟩ ≡_{1} 0}, (3)
Ê:= { ⟨(−b, a^{T})^{T}, x̂⟩ = 0 | (⟨a, x⟩ = b) ∈ E }.  (4)

The congruence ⟨e_{0}, x̂⟩ ≡_{1} 0, expressing the fact that 1 ≡_{1} 0, is called the integrality congruence. By writing F̂ = (F^{T}x ≡_{1} 0) and Ê = (E^{T}x = 0) where F, E ⊆ ℚ^{n+1}, it can be seen that the pair (F, E), called the matrix form of the homogeneous system Ĉ, is sufficient to determine C. Observe also that we can recover the original grid gcon(C) from gcon(Ĉ) since:

gcon(C) = { v ∈ ℝ^{n} | (1, v^{T})^{T} ∈ gcon(Ĉ) }.

Consider next a generator system G = (L, Q, P) in ℚ^{n}. Then the homogeneous form for G is the generator system Ĝ:= (L̂, Q̂ ∪ P̂, {0}) ∈ ℚ^{n+1} where

L̂ := {(0, ℓ^{T})^{T} | ℓ ∈ L}, Q̂ := {(0, q^{T})^{T} | q ∈ Q}, P̂ := {(1, p^{T})^{T} | p ∈ P }.

(5)

We can recover the original grid ggen(G) from ggen(Ĝ) since:

ggen(G) = { v ∈ ℝ^{n} | (1, v^{T})^{T} ∈ ggen(Ĝ) }.

Suppose that C is a congruence (resp., G is a generator) system in ℚ^{n+1}; then we say that C (resp., G) is in homogeneous form if there exists a congruence system Č (resp., generator system Ǧ) in ℚ^{n} such that C (resp., G) is the homogeneous form of Č (resp., Ǧ). Note that, if (C, G) is a double description for a grid and Ĉ and Ĝ are homogeneous forms for C and G, then (Ĉ, Ĝ) is also a double description.

5.2 Converting Representations

In this subsection, we outline the conversion algorithms (that is, algorithms for converting from the congruence to the generator minimized representations of a nonempty grid and vice versa).

By considering the matrix forms of the (homogeneous forms of the) representations, we can build the conversion algorithms on top of standard techniques for matrix inversion. For an informal explanation why this is appropriate, suppose that the generator system G = (∅, Q, {0}) in ℚ^{n} is in minimal form and that # Q = n. Then Q is a non-singular square matrix. Letting L = ggen(G) = { Qπ | π ∈ ℤ^{n} }, then we also have L = { v ∈ ℝ^{n} | Q^{−1}v ≡_{1} 0 }, so that (Q^{−1}, ∅) is the matrix form of a congruence system for the same grid L. Similarly we can use matrix inversion to convert the matrix form of a homogeneous congruence system in minimal form consisting of n proper congruences for a grid L to a generator system for L.

Proposition 7. Let C be a congruence system in ℚ^{n} in minimal form. Let (F, E) be the matrix form of the homogeneous form for C. Let N be a matrix in ℤ^{n+1} whose vectors are of the form e_{i}, i ∈ {0, …, n}, and such that (N, F̂, Ê) is square and nonsingular. Let

(L̂, Q̂, M) := ((N, F̂, Ê)^{−1})^{T}

where # L̂ = # N, # Q̂ = # F̂ and # M = # Ê. Then Ĝ = (L̂, Q̂, {0}) is the homogeneous form for a generator system G in minimal form and ggen(G) = gcon(C).

As matrix inversion and transposition are both computable and invertible, the conversion from generator to congruence system is similar.

Proposition 8. Let G be a generator system in ℚ^{n} in minimal form and Ĝ = (L̂, Q̂, {0}) be the homogeneous form for G. Let M be a matrix in ℤ^{n+1} whose vectors are of the form e_{i}, i ∈ {0, …, n}, and such that (L̂, Q̂, M) is square and nonsingular. Let

(N, F̂, Ê) := ((L̂, Q̂, M)^{−1})^{T}

where # N = # L̂, # F̂ = # Q̂ and # Ê = # M. Then (F̂, Ê) is the matrix form of the homogeneous form for a congruence system C in minimal form and gcon(C) = ggen(G).

Both conversion algorithms just perform matrix inversion so that the complexity is that of matrix inversion, namely O(n^{3}).

6 Grid Widening

In this section we describe a widening for the domain of grids, and compare our proposal with an operator based on the standard widening for the domain of convex polyhedra.

A simple and general characterization of a widening for enforcing and accelerating convergence of an upward iteration sequence is given in [7–10]. We assume here a minor variation of this classical definition (see footnote 6 in [10, p. 275]).

Definition 4. (Widening.) Let ⟨D, ⊢, 0, ⊕⟩ be a join-semilattice. The partial operator ∇: D × D ↣ D is a widening if

1. for each d_{1}, d_{2} ∈ D, d_{1} ⊢ d_{2} implies that d_{1} ∇ d_{2} is defined and d_{2} ⊢ d_{1} ∇ d_{2};

2. for each increasing chain d_{0} ⊢ d_{1} ⊢ ⋯, the increasing chain defined by d′_{0} := d_{0} and d′_{i+1} := d′_{i} ∇ (d′_{i} ⊕ d_{i+1}), for i ∈ ℕ, is not strictly increasing.

In addition to the formal requirements in Definition 4, it is also important to have a widening that has an efficient implementation, preferably, one that depends on a simple syntactic mapping of the representations. At the same time, so that the widening is well-defined, the result of this operation should be independent of the actual representation used. For this reason, the widening we propose requires the concept of strong minimal form for the congruence systems.

Definition 5. A congruence system C in ℚ^{n} is in strong minimal form if C is in minimal form and, for each pair of distinct proper congruences

β = (⟨a, x⟩ ≡_{1} b), γ = (⟨c, x⟩ ≡_{1} d) ∈ C,

if piv_{<}(γ) = k > 0, then −c_{k} < 2a_{k} ≤ c_{k}.

A congruence system in minimal form can always be reduced to a congruence system in strong minimal form that describes the same grid.

Proposition 9. Let C be a congruence system in ℚ^{n} in minimal form. Then there exists an algorithm with complexity O(n^{3}) for converting C to a congruence system C′ in strong minimal form such that C is pivot equivalent to C′ and gcon(C) = gcon(C′).

We also have that, relative to the affine hull of a grid, the strong minimal form of any congruence system that describes that grid is canonical.

Proposition 10. Let C_{1} and C_{2} be congruence systems in strong minimal form in ℚ^{n} such that L = gcon(C_{1}) = gcon(C_{2}). Then, for all β ∈ C_{1}, there exists γ ∈ C_{2} such that

gcon({β}) ∩ affine.hull(L) = gcon({γ}) ∩ affine.hull(L).

We now define a widening for grids.

Definition 6. Let L_{1} = gcon(C_{1}) and L_{2} = gcon(C_{2}) be two grids in 𝔾_{n} such that L_{1} ⊆ L_{2}, C_{1} is in minimal form and C_{2} is in strong minimal form. Then the grid widening L_{1} ∇ L_{2} is defined by

L_{1} ∇ L_{2} := \begin{cases}L_{2}, & if L_{1} = ∅ or dim(L_{1}) < dim(L_{2}); \\ gcon(C_{S}), & otherwise, \end{cases}

where

C_{S} := { γ ∈ C_{2} | ∃β ∈ C_{1} . β ⇑ γ }.  (6)

Proposition 11. The operator ‘∇’ is a widening on 𝔾_{n}.

Assuming that the two congruence systems are already available and in strong minimal form, the complexity of this widening is O(n^{2}).

By Proposition 10, for a fixed variable ordering, the operator ‘∇’ is welldefined. The next example illustrates that ‘∇’ depends on the variable ordering.

Example 8. Consider the grids L_{1} := gcon(C_{1}) and L_{2} := gcon(C_{2}) in 𝔾_{2}, where

C_{1} := {5x + y ≡_{1} 0, 22x ≡_{1} 0},  C_{2} := {5x + y ≡_{1} 0, 44x ≡_{1} 0}.

Assume that variables are ordered so that x precedes y, as in the vector (x, y)^{T}; then, the congruence systems C_{1} and C_{2} are in strong normal form and, according to Definition 6, we obtain L_{1} ∇ L_{2} = gcon({5x + y ≡_{1} 0}). On the other hand, the same grids as before are also described by the congruence systems

C′_{1} := {9y + x ≡_{1} 0, 22y ≡_{1} 0},  C′_{2} := {9y + x ≡_{1} 0, 44y ≡_{1} 0},

respectively, which are in strong normal form when considering the variable order where y precedes x. In this case, by Definition 6, L_{1} ∇L_{2} = gcon({9y+x ≡_{1} 0}).

In Definition 6, we require that the second congruence systems C_{2} is in strong minimal form. The following example shows that this is necessary for the operator ‘∇’ to be well-defined,

Example 9. Consider the congruence systems

C_{1} := {x ≡_{2} 0, y ≡_{2} 0},
C_{2} := {x ≡_{1} 0, x + y ≡_{2} 0},
C_{3} := {x ≡_{1} 0, 3x + y ≡_{2} 0},

and let L_{1} := gcon(C_{1}), L_{2} := gcon(C_{2}) and L_{3} := gcon(C_{3}); then L_{2} = L_{3}. Note that only C_{1} and C_{2} are in strong minimal form. Therefore, assuming C_{S} (resp., C_{S}′) is defined as in (6) using C_{1} and C_{2} (resp., C_{1} and C_{3}), we have

C_{S} = {x + y ≡_{2} 0},  C_{S}′ = {3x + y ≡_{2} 0}.

Thus L_{1} ∇ L_{2} = gcon(C_{S}) ≠ gcon(C_{S}′).

The grid widening ‘∇’ defined above is a fully functional widening operator and, as such, can be used for the development of more refined widenings (or extrapolation operators) by following well-known techniques and, in particular, the frameworks proposed in [2, 4].

For instance, when instantiating the finite powerset construction [4] with the domain of grids, a widening for the powerset of grids can be obtained by designing a finite convergence certificate for the grid widening. In formal terms, this certificate is a triple (O, ≻, μ) where (O, ≻) is a well-founded ordered set and μ: 𝔾_{n} → O, which is called level mapping, is such that, for all L_{1} ⊂ L_{2} ∈ 𝔾_{n}, μ(L_{1}) ≻ μ(L_{1} ∇ L_{2}). One such certificate can be easily defined by taking O equal to {0, …, n} × {0, …, n}, ≻ the lexicographic ordering on O and, for all L ∈ 𝔾_{n}, letting μ(L) := (# E, # C) where L = gcon(C), C is in minimal form, and E ⊆ C is the set of equalities in C. By Definition 6 and Proposition 2, if follows that L_{1} ≠ L_{2} implies μ(L_{1}) ≻ μ(L_{1} ∇ L_{2}); hence we have a finite convergence certificate for ‘∇’.

The grid widening ‘∇’ can be enhanced with additional delay constructs such as the with token option previously described for widenings on convex polyhedra [3]. Often in analysis or verification, the convergence guarantee that comes with a widening operator is not essential and in such cases, all that is required are extrapolation operators. These differ from widenings in that their use along an upper iteration sequence does not ensure convergence in a finite number of steps. In particular, following the widening “up to” technique as described in [20] for convex polyhedra, we define a limited extrapolation operator to be an operator that takes a congruence system as an additional parameter and uses it to improve the approximation yielded by the grid widening. We can also exploit the covering boxes, defined in Subsection 4.4, to provide a covered extrapolation operator that improves the approximation yielded by the widening operator by ensuring that the result cannot be worse than the covering box for the grid in the second argument of the widening.

6.1 Discussion

To motivate further our choice of ‘∇’ for the grid widening, consider the following grid counterpart of the standard widening for convex polyhedra as specified in the PhD thesis of N. Halbwachs [18], also described in [20].

Let L_{1} = gcon(C_{1}) and L_{2} = gcon(C_{2}) ∈ 𝔾_{n} where L_{1} ⊆ L_{2} and C_{1} and C_{2} are congruence systems in ℝ^{n} in strong minimal form. Then h(L_{1}, L_{2}) ∈ 𝔾_{n} is defined to be L_{2}, if L_{1} = ∅, and gcon(C′_{1} ∪ C′_{2}), if L_{1} ≠ ∅, where

C′_{1} := { β ∈ C_{1} | L_{2} ⊆ gcon({β}) },  (7)
C′_{2} := { γ ∈ C_{2} | ∃β ∈ C_{1} . L_{1} = gcon(C_{1}[γ/β]) }.  (8)

One difference here with the standard widening for convex polyhedra is in the handling of equalities. This is because, whereas an equality is equivalent to two inequalities, there is no finite set of proper congruences that is semantically equivalent to an equality, so the above definition has ignored any distinction between equalities and proper congruences.

A second difference is that to ensure that the operator ‘h’ is well-defined, we require that both the congruence systems C_{1} and C_{2} are in strong minimal form whereas the standard widening just requires that the constraint system describing the smaller polyhedron is minimized.

Example 10. Consider again the grids and congruence systems in Example 9. Then, assuming definition (7) for C′_{1} and (8) for C′_{2} and C′_{3}, we have C′_{1} = ∅, C′_{2} = {x + y ≡_{2} 0} and C′_{3} = {3x + y ≡_{2} 0} so that h(L_{1}, L_{2}) = C′_{1} ∪ C′_{2} = gcon({x + y ≡_{2} 0}) although C′_{1} ∪ C′_{3} = {3x + y ≡_{2} 0}.

Observe that the computation of the congruence system C′_{2} should be carefully tuned, since any naive implementation is going to be rather expensive. On the other hand, Example 10 illustrates that, simply ignoring the C′_{2} component (as done in some implementations of the standard widening for convex polyhedra) can lose precision, even when the affine hull of both the grids is the universe.

7 Applications

In this section we discuss applications for the domain of rational grids in program analysis and verification.

Many program properties are quantitative or depend on quantitative information. While such information may depend directly on the values of numerical data objects, it could instead reflect some numerical measures of the structure of the program and its data. We first discuss applications where the values of numeric variables are abstracted. In this case, any collection of numerical data objects (such as the numerical variables of some procedure) is a candidate for a non-trivial approximation via the grid domain.

Example 11. Consider again the example program given in Section 1 but annotated with program points Pj, for j = 1, …, 5:

    x := 2; y := 0;              (P1)
    for i := 1 to m              (P2)
      if ... then
        x := x + 4               (P3)
      else
        x := x + 2; y := y + 1   (P4)
      endif                      (P5)
    endfor

We show here that with the grid domain, we can derive the invariant set of congruences C = {x + 2y ≡_{4} 2, x ≡_{2} 0} that hold at point P2 and hence at end of the program, when the loop terminates. Let G_{j}^{i} ∈ 𝔾_{2} denote the grid computed at the i-th iteration executed by the point Pj. Initially, G_{j}^{0} = ∅ = gcon({1 = 0}), for j = 1, …, 5. After the first iteration of the analysis, we have:

G_{1}^{1} = gcon({x = 2, y = 0}),
G_{2}^{1} = gcon({x = 2, y = 0}),
G_{3}^{1} = gcon({x = 6, y = 0}),
G_{4}^{1} = gcon({x = 4, y = 1}),
G_{5}^{1} = gcon({x = 4, y = 1}) ⊕ gcon({x = 6, y = 0})
= gcon({x + 2y = 6, x ≡_{2} 0}).

8 6 4 2  2  4  6  8  10 x y

(a) Example 12: L_{r} and L_{w}

8 6 4 2  2  4  6  8  10 j  i

(b) Example 13: L, P and P′

Fig. 5. Examples 12 and 13

Then we have

G_{2}^{2} = gcon({x = 2, y = 0}) ⊕ gcon({x + 2y = 6, x ≡_{2} 0})
= gcon({x + 2y ≡_{4} 2, x ≡_{2} 0})

and subsequent steps in the computation show that we have already computed an invariant for P2 since we have G_{3}^{2} = G_{3}^{1}, G_{4}^{2} = G_{4}^{1}, G_{5}^{2} = G_{5}^{1} so that G_{3}^{2} = G_{2}^{2}.

Example 11 shows that, even for a very simple piece of code, using the grids domain, we can find non-trivial relational congruence properties that cannot be found using the polyhedra domain [11], constraint-based analysis [34] or polynomial invariants [33].

Data dependence analysis for arrays —deciding if two elements of an array can refer to the same element and, if so, under what conditions— is required for advanced optimizing compilers [30].

Example 12. Consider the following program (adapted from a simple example given in [30]):

    for i := 0 to 100
      for j := 2i to 100
        A[i, 2j + 1] := A[i, 2j]
      endfor
    endfor

Then, the program reads from array elements (0, 0), (0, 2), (1, 4), … and writes to array elements (0, 1), (0, 3), (1, 5), …. The two sets of points generate, respectively, the two grids L_{r} and L_{w} in ℝ^{2}: L_{r} = ggen{(0, 0), (0, 2), (1, 4)} includes all the array elements that are read from, while L_{w} = ggen{(0, 1), (0, 3), (1, 5)} includes all the array elements that are written to. Notice that, by computing the intersection L_{r} ∩L_{w} we can deduce that no location is both read and written.

Figure 5(a) illustrates the grids L_{r} and L_{w} where squares denote the points of the grid L_{r} and circles denote the points of the grid L_{w}.

As mentioned in Section 1, most of the previous work on domains for numerical information concerns the representation of bounds on the values of the data objects. When one of these bounds is strict, an implementation normally has to use a different and, usually, less efficient representation compared with that for domain elements that are known to be topologically closed (and therefore defined entirely by non-strict constraints). However, if the data objects can only take discrete values and we have information about the frequency of their distribution, we can convert all the inequalities into non-strict ones, allowing for a more efficient implementation as well as possibly improving the precision.

Example 13. Consider the following program fragment:

    for i := 1 to 10
      j := 0
      while j < i
        A[j] := A[j + 2];
        j := j + 2
      endwhile
    endfor

Then the linear constraints representing the values of i and j at the start of each iteration of the while loop will lead to a convex polyhedron bounded by constraints {1 ≤ i ≤ 10, 0 ≤ j < i}. If we add the information that i, j must be integers, we obtain a polyhedron P bounded by constraints {1 ≤ i ≤ 10, 0 ≤ j ≤ i − 1}; so that it may be deduced that at most 9 elements of the vector A are changed. However, using a grid domain we can represent the distribution of the possible values of i and j by a grid L described by the congruences {i ≡_{1} 0, j ≡_{2} 0}; so that, by combining this information with that represented by P, we can see that the values of i and j must lie in the polyhedron P′ bounded by constraints {1 ≤ i ≤ 10, 0 ≤ j ≤ i − 1, j ≤ 8} and deduce that at most 8 elements of the vector A are changed since the value of A[9] is unaltered.

Figure 5(b) illustrates L, P, and P′ indicated, respectively, by the squares, the shaded area, and the region enclosed by the dashed lines.

Example 14. This example is taken from [19] and concerns a water tank monitor. Let t be the time in seconds and w the water level in the tank in meters. Let p denote the pump switch position (on or off). When the pump is running the water fills the tank at 1 meter a second. When the pump is off, the water empties from the tank at 2 meters a second. At the start when t = 0, the water level w = 1 and the pump position p = on. By changing the position of the pump switch, the monitor ensures that the water level never reaches 12 meters and

12 10 8 6 4 2  1  2  3  4  5  6  7  8  9  10 t w

Fig. 6. Water tank monitor example

that there is always 2 meters of water in the tank. It checks the water level every half second. The procedure is described by the following program where T is a constant denoting the maximum time:

    t = 0; w = 1; p = on;
    while t < T
      t = t + 0.5;
      if p = on then
        w := w + 0.5;
        if w >= 11.5
          p := off
      else
        w := w - 1;
        if w <= 2
          p := on
        endif
      endif
    endwhile

This will result in a grid of points relating the time t elapsed and the level of the water w. As the monitor checks every half second, we just have the discrete points at half second intervals. The grid is therefore generated by the set of points {(0, 1), (1, 2), (10.5, 11.5), (11, 10.5)}. This is also described by the set of congruences {2t ≡_{1} 0, 2t − 2w ≡_{3} 1}; the grid gcon({2t ≡_{1} 0, 2t − 2w ≡_{3} 1}) is illustrated in Figure 6.

8 Conclusion

We have presented the abstract domain of grids, a domain that allows to represent and manipulate sets of regular spaced points and hyperplanes over ℝ^{n}. Two alternative grid representations in ℚ^{n} are described: one, a set of congruences, where the congruences can be either equalities or proper congruences; and the other, a set of generators, which can be lines, parameters or points.

We have provided algorithms for minimizing and converting the grid representations. Assuming the systems in ℚ^{n} consist of m congruences or generators, the minimization algorithms have complexity O(n^{2}m). Note that other, previously proposed algorithms for minimization, which only consider the addition of one generator at a time, have complexity, at best, O(n^{3}) [27]. This is comparable with ours since, in this case, m = O(n). We have shown that existing algorithms for matrix inversion can be used for converting between generators and congruence systems so that these conversion procedures have complexity O(n^{3}) (for congruence to generator conversion this is an improvement over the previous proposals, which have complexity no better than O(n^{4}) [17]).

We have observed that the grid meet (resp., join) operations can be easily implemented if we have the congruence (resp. generator) representations for the grids as we then just need to take the union of the systems; hence both the operators have complexity O(n^{3}). In contrast, in [16] the grid join operation has complexity O(n^{4} log_{2} n) since the generators of one of the grids are added, one at a time, to the generators of the other grid, reducing the new set at each stage to a linearly independent set. Moreover, because of the complexity of the conversion algorithm which is needed to convert the representation to a congruence system, the complexity of the grid meet operation in [16] is also O(n^{4}).

An algorithm for computing the grid difference, the best approximation by a grid for the difference of two grids, has been proposed and proved to be correct (see Section 4). We have proposed a grid widening operator for rational grids that is straightforward to implement. To ensure it is well-defined, the congruence system representing the grid has to be in a strong minimal form; that is, one in minimal form for which the non-diagonal elements of the matrix of coefficients of the proper congruences have to have a minimal absolute value. We also discuss how the proposed widening compares with the grid counterpart of the standard widening for convex polyhedra [20]. We have based the widening on the congruence representations of the grids. However, we anticipate that a dual form of the widening could also be defined that uses the generator representations instead; further investigation of this is future work.

Apart from a widening, we also show how the polyhedral domain of intervals may be used to represent a non-relational grid. We have also provided algorithms for finding the best approximation of an arbitrary relational grid as a non-relational grid. We observe that this idea of using a polyhedral domain to represent the “holes” in a grid could easily be generalized from that of representing a rectangular grid using covering boxes to the representation of any relational grid using n-dimensional parallelograms.

The grid domain is being implemented in the PPL [5, 6] following the approach and algorithms described in this paper. The current code is already publicly available in the grids branch of the PPL’s CVS repository (see the PPL web site at http://www.cs.unipr.it/ppl/ for more information). As the grid domain is implemented in a library providing several polyhedral domains, it can be combined with any of these domains to give not only the Z-polyhedra domain but also many variations such as grid-polyhedra, grid-octagons, grid-bounded-differences and grid-intervals. Moreover, as the PPL provides full support for lifting any domain to the powerset of that domain, a user of the PPL will be able to experiment with powersets of grids (and of grid-polyhedra combinations) and the extra precision this provides.

References

1. C. Ancourt. Génération Automatique de Codes de Transfert pour Multiprocesseurs à Mémoires Locales. PhD thesis, Université de Paris VI, March 1991.

2. R. Bagnara, P. M. Hill, E. Ricci, and E. Zaffanella. Precise widening operators for convex polyhedra. In R. Cousot, editor, Static Analysis: Proceedings of the 10th International Symposium, volume 2694 of Lecture Notes in Computer Science, pages 337–354, San Diego, California, USA, 2003. Springer-Verlag, Berlin.

3. R. Bagnara, P. M. Hill, E. Ricci, and E. Zaffanella. Precise widening operators for convex polyhedra. Science of Computer Programming, 58(1–2):28–56, 2005.

4. R. Bagnara, P. M. Hill, and E. Zaffanella. Widening operators for powerset domains. In B. Steffen and G. Levi, editors, Verification, Model Checking and Abstract Interpretation: Proceedings of the 5th International Conference (VMCAI 2004), volume 2937 of Lecture Notes in Computer Science, pages 135–148, Venice, Italy, 2003. Springer-Verlag, Berlin.

5. R. Bagnara, P. M. Hill, and E. Zaffanella. The Parma Polyhedra Library User’s Manual. Department of Mathematics, University of Parma, Parma, Italy, release 0.8 edition, January 2006. Available at http://www.cs.unipr.it/ppl/.

6. R. Bagnara, E. Ricci, E. Zaffanella, and P. M. Hill. Possibly not closed convex polyhedra and the Parma Polyhedra Library. In M. V. Hermenegildo and G. Puebla, editors, Static Analysis: Proceedings of the 9th International Symposium, volume 2477 of Lecture Notes in Computer Science, pages 213–229, Madrid, Spain, 2002. Springer-Verlag, Berlin.

7. P. Cousot and R. Cousot. Static determination of dynamic properties of programs. In B. Robinet, editor, Proceedings of the Second International Symposium on Programming, pages 106–130, Paris, France, 1976. Dunod, Paris, France.

8. P. Cousot and R. Cousot. Abstract interpretation: A unified lattice model for static analysis of programs by construction or approximation of fixpoints. In Proceedings of the Fourth Annual ACM Symposium on Principles of Programming Languages, pages 238–252, New York, 1977. ACM Press.

9. P. Cousot and R. Cousot. Abstract interpretation frameworks. Journal of Logic and Computation, 2(4):511–547, 1992.

10. P. Cousot and R. Cousot. Comparing the Galois connection and widening/narrowing approaches to abstract interpretation. In M. Bruynooghe and M. Wirsing, editors, Proceedings of the 4th International Symposium on Programming Language Implementation and Logic Programming, volume 631 of Lecture Notes in Computer Science, pages 269–295, Leuven, Belgium, 1992. Springer-Verlag, Berlin.

11. P. Cousot and N. Halbwachs. Automatic discovery of linear restraints among variables of a program. In Conference Record of the Fifth Annual ACM Symposium on Principles of Programming Languages, pages 84–96, Tucson, Arizona, 1978. ACM Press.

12. R. Giacobazzi, editor. Static Analysis: Proceedings of the 11th International Symposium, volume 3148 of Lecture Notes in Computer Science, Verona, Italy, 2004. Springer-Verlag, Berlin.

13. D. Gopan, F. DiMaio, N. Dor, T. Reps, and M. Sagiv. Numeric domains with summarized dimensions. In K. Jensen and A. Podelski, editors, Tools and Algorithms for the Construction and Analysis of Systems, 10th International Conference, TACAS 2004, volume 2988 of Lecture Notes in Computer Science, pages 512–529, Barcelona, Spain, 2004. Springer-Verlag, Berlin.

14. P. Granger. Static analysis of arithmetical congruences. International Journal of Computer Mathematics, 30:165–190, 1989.

15. P. Granger. Analyses Sémantiques de Congruence. PhD thesis, École Polytechnique, 921128 Palaiseau, France, July 1991.

16. P. Granger. Static analysis of linear congruence equalities among variables of a program. In Samson Abramsky and T. S. E. Maibaum, editors, TAPSOFT’91: Proceedings of the International Joint Conference on Theory and Practice of Software Development, Volume 1: Colloquium on Trees in Algebra and Programming (CAAP’91), volume 493 of Lecture Notes in Computer Science, pages 169–192, Brighton, UK, 1991. Springer-Verlag, Berlin.

17. P. Granger. Static analyses of congruence properties on rational numbers (extended abstract). In P. Van Hentenryck, editor, Static Analysis: Proceedings of the 4th International Symposium, volume 1302 of Lecture Notes in Computer Science, pages 278–292, Paris, France, 1997. Springer-Verlag, Berlin.

18. N. Halbwachs. Détermination Automatique de Relations Linéaires Vérifiées par les Variables d’un Programme. Thèse de 3^{ème} cycle d’informatique, Université scientifique et médicale de Grenoble, Grenoble, France, March 1979.

19. N. Halbwachs, Y.-E. Proy, and P. Raymond. Verification of linear hybrid systems by means of convex approximations. In B. Le Charlier, editor, Static Analysis: Proceedings of the 1st International Symposium, volume 864 of Lecture Notes in Computer Science, pages 223–237, Namur, Belgium, 1994. Springer-Verlag, Berlin.

20. N. Halbwachs, Y.-E. Proy, and P. Roumanoff. Verification of real-time systems using linear relation analysis. Formal Methods in System Design, 11(2):157–185, 1997.

21. M. Karr. Affine relationships among variables of a program. Acta Informatica, 6:133–151, 1976.

22. S. Larsen, E. Witchel, and S. P. Amarasinghe. Increasing and detecting memory address congruence. In Proceedings of the 2002 International Conference on Parallel Architectures and Compilation Techniques (PACT’02), pages 18–29, Charlottesville, VA, USA, 2002. IEEE Computer Society Press.

23. V. Loechner. PolyLib: A library for manipulating parameterized polyhedra. Avail-

able at http://icps.u-strasbg.fr/^{∼}loechner/polylib/, March 1999. Declares

itself to be a continuation of [37].

24. A. Miné. A few graph-based relational numerical abstract domains. In M. V. Hermenegildo and G. Puebla, editors, Static Analysis: Proceedings of the 9th International Symposium, volume 2477 of Lecture Notes in Computer Science, pages 117–132, Madrid, Spain, 2002. Springer-Verlag, Berlin.

25. T. S. Motzkin, H. Raiffa, G. L. Thompson, and R. M. Thrall. The double description method. In H. W. Kuhn and A. W. Tucker, editors, Contributions to the Theory of Games – Volume II, number 28 in Annals of Mathematics Studies, pages 51–73. Princeton University Press, Princeton, New Jersey, 1953.

26. M. Müller-Olm and H. Seidl. Analysis of modular arithmetic. In M. Sagiv, editor, Programming Languages and Systems, Proceedings of the 14th European Symposium on Programming, volume 3444 of Lecture Notes in Computer Science, pages 46–60, Edinburgh, UK, 2005. Springer-Verlag, Berlin.

27. M. Müller-Olm and H. Seidl. A generic framework for interprocedural analysis of numerical properties. In C. Hankin and I. Siveroni, editors, Static Analysis: Proceedings of the 12th International Symposium, volume 3672 of Lecture Notes in Computer Science, pages 235–250, London, UK, 2005. Springer-Verlag, Berlin.

28. G. L. Nemhauser and L. A. Wolsey. Integer and Combinatorial Optimization. Wiley Interscience Series in Discrete Mathematics and Optimization. John Wiley & Sons, 1988.

29. S. P. K. Nookala and T. Risset. A library for Z-polyhedral operations. Publication interne 1330, IRISA, Campus de Beaulieu, Rennes, France, 2000.

30. W. Pugh. A practical algorithm for exact array dependence analysis. Communications of the ACM, 35(8):102–114, 1992.

31. P. Quinton, S. Rajopadhye, and T. Risset. On manipulating Z-polyhedra. Technical Report 1016, IRISA, Campus Universitaire de Bealieu, Rennes, France, July 1996.

32. P. Quinton, S. Rajopadhye, and T. Risset. On manipulating Z-polyhedra using a canonic representation. Parallel Processing Letters, 7(2):181–194, 1997.

33. E. Rodríguez-Carbonell and D. Kapur. An abstract interpretation approach for automatic generation of polynomial invariants. In Giacobazzi [12], pages 280–295.

34. S. Sankaranarayanan, H. Sipma, and Z. Manna. Constraint-based linear-relations analysis. In Giacobazzi [12], pages 53–68.

35. A. Schrijver. Theory of Linear and Integer Programming. Wiley Interscience Series in Discrete Mathematics and Optimization. John Wiley & Sons, 1999.

36. G. Stewart. Introduction to Matrix Computations. Academic Press, New York, 1973.

37. D. K. Wilde. A library for doing polyhedral operations. Master’s thesis, Oregon State University, Corvallis, Oregon, December 1993. Also published as IRISA Publication interne 785, Rennes, France, 1993.

Appendix

A Proofs

Here we prove those results stated but not proved in the main part of the paper.

In the following two proofs, we require an extended version of the Euclidean algorithm. Suppose v, v′ ∈ ℤ. Let r ∈ ℤ be the greatest common divisor of v, v′ and s, s′ ∈ ℤ be relatively prime such that sv + s′v′ = r. Then we write

gcd(v, v′) := r, gcdext(v, v′) := (r, (s, s′)).

Proof (of Proposition 1). To prove the result, we first define the key transformation step in the algorithm and show that the resulting congruence system describes the same grid. Suppose there exist distinct congruences

β_{1} = (⟨a_{1}, x⟩ ≡_{f_{1}} b_{1}),  β_{2} = (⟨a_{2}, x⟩ ≡_{f_{2}} b_{2})  (9)

in C such that piv_{<}(β_{1}) = piv_{<}(β_{2}) = i > 0. We will define a congruence system C′′ such that either C′′ = C \ {β_{1}, β_{2}} ∪ {β′′_{1}, β′′_{2}} or C′′ = C \ {β_{1}, β_{2}} ∪ {β′′_{1}} and show that gcon(C) = gcon(C′′). Letting

β′′_{1} = (⟨a′′_{1}, x⟩ ≡_{f_{1}} b′′_{1}),  β′′_{2} = (⟨a′′_{2}, x⟩ ≡_{f_{2}} b′′_{2}),

we show that piv_{<}(a′′_{1}) = i and, if β′′_{2} is defined, then piv_{<}(a′′_{2}) < i. There are two cases.

1. Either β_{1} or β_{2} is an equality; without loss of generality, we assume that β_{1} is an equality so that f_{1} = 0. Then, we let β′′_{1} = β_{1} and, using Gaussian elimination,

a′′_{2} = a_{2} − (a_{2i}/a_{1i})a_{1}, b′′_{2} = b_{2} − (a_{2i}/a_{1i})b_{1}.

2. Both β_{1}, β_{2} are proper congruences; so that we can assume that f_{1} = f_{2} = 1. Let gcdext(a_{1i}, a_{2i}) = (r, (s, t)) Let

a′′_{1} = sa_{1} + ta_{2},  b′′_{1} = sb_{1} + tb_{2}
a′′_{2} = (−a_{2i}/r)a_{1} + (a_{1i}/r)a_{2},  b′′_{2} = (−a_{2i}/r)b_{1} + (a_{1i}/r)b_{2}.

In both cases, let C′′ = C \ {β_{1}, β_{2}} ∪ {β′′_{1}, β′′_{2}} if a′′_{2} ≠ 0 or b′′_{2} ≠ 0, and let C′′ = C \ {β_{1}, β_{2}} ∪ {β′′_{1}}, otherwise. Then gcon(C) = gcon(C′′). Note that these transformations require a computation for each coefficient of the considered congruences so that their complexity is O(n).

The proof of the result is by induction on i; where 0 ≤ i ≤ n is the maximum value for which there exist distinct congruences β_{1}, β_{2} ∈ C defined as in (9) such that piv_{<}(β_{1}) = piv_{<}(β_{2}) = i.

The base case is when i = 0 so that a_{1} = a_{2} = 0. In this case, if there exists (⟨0, x⟩ ≡_{f} b) ∈ C and b ≡_{f} 0 is inconsistent, let C′ = {⟨0, x⟩ ≡_{0} 1}; otherwise, let C′ = { β ∈ C | piv_{<}(β) ≠ 0 }.

For the step case, we keep applying the transformation (1) if either f_{1} = 0 or f_{2} = 0, and (2) otherwise until no more transformations are applicable for this index; that is when we obtain a congruence system C_{j} for which j < i is the maximum index such that there exist distinct congruences β_{1}, β_{2} ∈ C_{j} where piv_{<}(β_{1}) = piv_{<}(β_{2}) = j. We note that we will have to perform these transformations at most m times for each step, where # C = m, so that the complexity of each step is O(nm). By the inductive hypothesis, we can compute C′ in minimal form such that gcon(C′) = gcon(C_{j}). Therefore gcon(C′) = gcon(C). As we iterate at most n times over the step case, it can be seen that the algorithm has complexity O(n^{2}m). ⊓ ⊔

Proof (of Proposition 3). If P = ∅, then ggen(G) = ∅; in this case, let G′ = (∅, ∅, ∅). Suppose now that there exists a point p ∈ P. Let Q_{p} = {p′′ − p | p′′ ∈ P \ {p}} ∪ Q and let G_{p} = (L, Q_{p}, {p}). Then, by equation (1), ggen(G_{p}) = ggen(G).

To prove the result, we first define the key transformation step in the algorithm and show that the resulting generator system describes the same grid. Suppose there exist distinct generators v_{1}, v_{2} ∈ L ∪ Q_{p} such that piv_{>}(v_{1}) = piv_{>}(v_{2}) = i ≤ n. We will define a generator system G′′ = (L′′, Q′′_{p}, {p}) in ℚ^{n} where L′′ ∪ Q′′_{p} = L ∪ Q_{p} \ {v_{1}, v_{2}} ∪ ({v′′_{1}, v′′_{2}} \ {0}), ggen(G′′) = ggen(G_{p}), piv_{>}(v′′_{1}) = i and, if v′′_{2} ≠ 0, piv_{>}(v′′_{2}) > i. There are two cases.

1. Either v_{1} or v_{2} is in L; without loss of generality, we assume that v_{1} ∈ L. Then, using Gaussian elimination, let v′′_{1} = v_{1} and

v′′_{2} = v_{2} − (v_{2i}/v_{1i})v_{1}.

Let L′′ = L\{v_{2}}∪({v′′_{2}}\{0}), if v_{2} ∈ L, and Q′′_{p} = Q_{p}\{v_{2}}∪({v′′_{2}}\{0}), if v_{2} ∈ Q_{p}.

2. Both v_{1}, v_{2} are in Q_{p}. Let gcdext(v_{1i}, v_{2i}) = (r, (s, t)),

v′′_{1} = sv_{1} + tv_{2},  v′′_{2} = (−v_{2i}/r)v_{1} + (v_{1i}/r)v_{2},
L′′ = L,  Q′′_{p} = Q_{p} \ {v_{1}, v_{2}} ∪ ({v′′_{1}, v′′_{2}} \ {0}).

In both cases, ggen(G′′) = ggen(G_{p}). Note that these transformations require a computation for each coefficient of the considered generators so that their complexity is O(n).

The proof of the result is by induction on n + 1 − i where i is the least value such that, if there exists a pair of distinct generators v_{1}, v_{2} ∈ L ∪ Q_{p}, then i = piv_{>}(v_{1}) = piv_{>}(v_{2}) and i = n + 1, otherwise.

The base case is when i = n + 1, in which case G_{p} is already in minimal form so let G′ = G_{p}.

For the step case, we apply the transformations (1) and (2) until no more transformations are applicable with index i; that is when we obtain a generator system G_{j} = (L_{j}, Q_{j}, {p}) for which j > i is the least value such that, if there exists a pair of distinct generators v_{1}, v_{2} ∈ L_{j} ∪Q_{j}, then j = piv_{>}(v_{1}) = piv_{>}(v_{2}) and j = n + 1, otherwise. We note that we will have to perform these transformations at most m times for each step, where # L+# Q_{p} = m, so that the complexity of each step is O(nm). By the inductive hypothesis, we can compute G′ in minimal form such that ggen(G′) = ggen(G_{j}). Therefore ggen(G′) = ggen(G). As we iterate at most n times over the step case, the algorithm has complexity O(n^{2}m). ⊓ ⊔

Proof (of Proposition 5). If L_{1} = ∅ or L_{2} = ∅, then L_{1} ⊖ L_{2} = L_{1} \L_{2} = L_{1} and, by lines (1-2), the algorithm returns L = L_{1}. Similarly, if L_{1} ⊆ L_{2}, then L_{1} ⊖ L_{2} = L_{1} \ L_{2} = ∅ and, by lines (3-4), the algorithm returns L = ∅. Suppose now that L_{1} ≠ ∅, L_{2} ≠ ∅ and L_{1} ⊈ L_{2}. Let L′ be the empty grid in 𝔾_{n} defined on line (5). Then the algorithm executes lines (6-14). Notice that there are two lines in this range that return a value for L; line (13) when L = L_{1} and line (14) when L = L′.

Consider first the case when line (13) is executed so that L = L_{1}. By definition of grid difference, L ⊇ L_{1} ⊖ L_{2}; Hence it remains to show that L_{1} ⊆ L_{1} ⊖ L_{2}. If p ∈ L_{1} \ L_{2}, then, by the definition of grid difference, p ∈ L_{1} ⊖ L_{2}. Suppose now that p ∈ L_{1} ∩ L_{2}. As line (13) is only executed by following the else branches of the conditionals on lines (8) and (9), for some congruence β = (e ≡_{f} 0) ∈ C_{2}, there exists a point q ∈ L_{1} that does not satisfy (2e ≡_{f} 0) so that q does not satisfy β and hence q ∉ L_{2}. Consider the point r = p + 2(q − p). Then, as r is an integral affine combination of points in L_{1}, r ∈ L_{1}. Let e = (⟨a, x⟩ − b). Then, as p ∈ L_{2} satisfies β, ⟨a, p⟩ − b ≡_{f} 0. If r also satisfies β, then ⟨a, r⟩ − b ≡_{f} 0 and hence ⟨a, 2q⟩ − 2b ≡_{f} 0 so that q would satisfy (2e ≡_{f} 0); a contradiction. Thus r ∉ L_{2}. Therefore p = 2q − r is an integral affine combination of points in L_{1} \ L_{2} and hence p ∈ L_{1} ⊖ L_{2}. As p ∈ L_{1} = L was arbitrary, L ⊆ L_{1} ⊖ L_{2}.

Suppose now that line (13) is not executed. Then the loop iterates once for each congruence in C_{2} before executing line (14). Suppose # C_{2} = c and β_{i} = (e_{i} ≡_{f} 0) ∈ C_{2} is the congruence selected at line (6) in the i-th iteration of the loop, for 0 < i ≤ c. Let L′_{0} = ∅ and L′_{i} denote the grid L′ after the i-th iteration. Then we need to show that L′_{c} = L_{1} ⊖ L_{2}. We prove that L′_{c} ⊆ L_{1} ⊖ L_{2} and L′_{c} ⊇ L_{1} ⊖ L_{2} separately.

We first show that L′_{c} ⊇ L_{1} ⊖ L_{2}. Since L_{1} ⊖ L_{2} is the smallest grid containing L_{1} \ L_{2}, we just need to show that L′_{c} ⊇ L_{1} \ L_{2}. To do this, let p ∈ L_{1} \ L_{2}; then we prove that p ∈ L′_{c}. As p ∉ L_{2}, for some j = 1, …, c, p ∉ gcon({β_{j}}). Consider the j-th iteration of the loop. Then the test on line (8) will succeed and the execution continues with the test on line (9). Moreover, as we know that line (13) will not be executed, this test must succeed so that p ∈ gcon({2e_{j} ≡_{f} 0}) and lines (10-11) will be executed with β = β_{j}. As gcon({β_{j}}) and gcon({2e_{j} ≡_{f} f}) are disjoint and their set union is the grid gcon({2e_{j} ≡_{f} 0}), p must satisfy the congruence (2e_{j} − f ≡_{2f} 0). Let L_{β_{j}} = gcon(C_{1} ∪ {2e_{j} − f ≡_{2f} 0}) as on line (10). Then, as p ∈ L_{1}, we have p ∈ L_{β_{j}}; hence, after line (11), p ∈ L′_{j}. For each i = j + 1, …, c, either L′_{i} = L′_{i−1} or line (10) is executed, in which case L′_{i} ⊇ L′_{i−1}; hence p ∈ L′_{i}. In particular, p ∈ L′_{c}. As this holds for all p ∈ L_{1} \ L_{2}, L′_{c} ⊇ L_{1} \ L_{2}.

Finally we prove, by induction on i, that, for each i = 0, …, c, L′_{i} ⊆ L_{1} ⊖ L_{2}. Initially L′_{0} = ∅ and the result holds. Suppose now that i > 0 and that L′_{i−1} ⊆ L_{1} ⊖ L_{2}. If L_{1} ⊆ gcon({β_{i}}), then L′_{i} = L′_{i−1} is unchanged by the iteration. On the other hand, if L_{1} ⊈ gcon({β_{i}}), the test on line (8) will succeed and the execution continues with the test on line (9). Moreover, as we know that line (13) will not be executed, this test must succeed so that L_{1} ⊆ gcon({2e_{i} ≡_{f} 0}). Let L_{β_{i}} = gcon(C_{1} ∪ {2e_{i} − f ≡_{2f} 0}) as defined on line (10); then L_{β_{i}} ∩ L_{2} = ∅ so that L_{β_{i}} ⊆ L_{1} \ L_{2} ⊆ L_{1} ⊖ L_{2}. Since, on line (11), L′_{i} is assigned L′_{i−1} ⊕ L_{β_{i}}, by definition of grid join and grid difference, L′_{i} ⊆ L_{1} ⊖ L_{2}. Therefore, letting i = c, we have L′_{c} ⊆ L_{1} ⊖ L_{2}. ⊓ ⊔

Proof (of Proposition 6). Let C be a congruence system for L and let L_{R} be the smallest rectilinear grid that contains L. Suppose also that B defines the rectilinear grid L_{B}. Then we need to show that L_{B} = L_{R}. Let C_{B} be the congruence system for L_{B} defined by equation (2) and C_{R} be congruence system in minimal form for L_{R}. Observe that L_{R} = ggen(G_{R}) where G_{R} is a minimal form for the generator system

({ℓ_{i}e_{i}|ℓ ∈ L, ℓ_{i} ≠ 0}, {q_{i}e_{i}|q ∈ Q, q_{i} ≠ 0}, {p}).

To prove the result, we consider each dimension i where 1 ≤ i ≤ n separately and show that there exists γ_{i} = (⟨e_{i}, x⟩ ≡_{f} b) ∈ C_{R} for some f, b ∈ ℚ if and only if β_{i} exists and gcon(γ_{i}) = gcon(β_{i}). We consider each of the three cases (1), (2) and (3) separately.

Suppose, for some ℓ ∈ L, ℓ_{i} ≠ 0. Then case (1) applies and I_{i} = (0, 0). Moreover, there is a line e_{i}ℓ_{i} in G_{R}. Thus, there are no congruences ⟨a, x⟩ ≡_{f} b in C_{B} or C_{R} where piv_{<}(a) = i.

Suppose, for all ℓ ∈ L, ℓ_{i} = 0, and q_{1i} = ⋯ = q_{qi} = 0. Then case (2) applies and I_{i} = (p_{i}, ∞). Therefore there exists β_{i} = (⟨e_{i}, x⟩ = p_{i}) ∈ C_{B}. Moreover, there are no lines or parameters in G_{R} with a non-zero i-th coordinate. Thus there exists γ_{i} ∈ C_{R} where f = 0 and b = p_{i} so that gcon(γ_{i}) = gcon(β_{i}).

Suppose, for all ℓ ∈ L, ℓ_{i} = 0, and, for some j = 1, …, q, q_{ji} ≠ 0. Then case (3) applies and I_{i} = (p_{i}, p_{i} + |g|) where g = gcd({q_{1i}, …, q_{qi}}) ≠ 0. Therefore there exists β_{i} = (⟨e_{i}, x⟩ ≡_{g} p_{i}) ∈ C_{B}. Moreover, there are no lines in G_{R} with a non-zero i-th coordinate. Thus there exists γ_{i} ∈ C_{R} where f = g and b = p_{i} + mg for some m ∈ ℤ. Thus gcon(γ_{i}) = gcon(β_{i}). ⊓ ⊔

In the following lemma, letting H ∈ ℝ^{n×n}, we say that H is upper triangular form if, for all i = 1, …, n, H_{ii} ≠ 0 and, for all j where 1 ≤ j < i ≤ n, H_{ij} = 0. Similarly, a matrix H is in lower triangular form if, for all i = 1, …, n, H_{ii} ≠ 0 and, for all j where 1 ≤ i < j ≤ n, H_{ij} = 0. Thus the transpose of an upper triangular matrix is lower triangular and vice-versa.

Lemma 1. There exists a computable, invertible function that converts a generator system G = (L, Q, {p}) in ℚ^{n} in minimal form to a congruence system C = E ∪ F in ℚ^{n} in minimal form where E are equalities and F are proper congruences and such that 1. # Q = # F = n − # L − # E;

2. there exists q ∈ Q if and only if there exists β = (⟨a, x⟩ ≡_{1} 0) ∈ F, such that, for some k ∈ {1, …, n}, piv_{>}(q) = piv_{<}(a) = k and q_{k}a_{k} = 1;

3. for all ℓ ∈ L and β = (⟨a, x⟩ = 0) ∈ E, piv_{>}(ℓ) ≠ piv_{<}(a);

4. gcon(C) = ggen(G).

Proof. Let Ĝ be the homogeneous form for G and Ĉ = F̂ ∪Ê be the homogeneous form for C where G is as defined in equation (5) and Ê and F̂ are as defined in equations (3) and (4) in Section 5. Let Ê = Ê^{T}x̂ = 0 and F̂ = F̂^{T}x̂ ≡_{1} 0 where x̂ = (x_{0}, …, x_{n})^{T}. Then, as G and C are in minimal form, by the definition of homogeneous form, items (1), (2), (3), and (4) hold if and only if the following hold:

5. # Q̂ = # F̂ = n + 1 − # L̂ − # Ê;

6. there exists q ∈ Q̂ if and only if there exists β = (⟨a, x̂⟩ ≡_{1} 0) ∈ F̂, such that, for some k ∈ {0, …, n}, piv_{>}(q) = piv_{<}(a) = k and q_{k}a_{k} = 1;

7. for all ℓ ∈ L̂ and β = (⟨a, x̂⟩ = 0) ∈ Ê, piv_{>}(ℓ) ≠ piv_{<}(a);

8. gcon(Ĉ) = ggen(Ĝ)

It therefore remains to prove that (5), (6), (7), and (8) hold.

Let M, N be matrices in ℤ^{n+1}, whose vectors are of the form e_{i}, 1 ≤ i ≤ n, and such that (L̂, Q̂, M) and (N, F̂, Ê) are square and nonsingular. Then we let

(N, F̂, Ê) := ((L̂, Q̂, M)^{−1})^{T}

where # N = # L̂, # Q̂ = # F̂ and # Ê = # N; ensuring (5) holds. Note that matrix inversion and transposition are both computable and invertible.

Let ℓ = # L̂, q = # Q̂ − 1 and m = # M. Then

x̂ ∈ Ĝ ⟺ x̂ = Q̂π + L̂λ + M0,  for λ ∈ ℝ^{ℓ}, π ∈ ℤ^{q+1}, 0 ∈ ℤ^{m}
⟺ x̂ = (Q̂, L̂, M)(π^{T}, λ^{T}, 0^{T})^{T},  for λ ∈ ℝ^{ℓ}, π ∈ ℤ^{q+1}, 0 ∈ ℤ^{m}
⟺ (Q̂, L̂, M)^{−1}x̂ = (π^{T}, λ^{T}, 0^{T})^{T}, for λ ∈ ℝ^{ℓ}, π ∈ ℤ^{q+1}, 0 ∈ ℤ^{m}.

As Ĝ is in minimal form (so that the set of vectors in (L̂, Q̂) is a subset of the set of vectors in a triangular matrix) and M adds the missing columns with only the diagonal elements being non-zero, (Q̂, L̂, M) is a permutation of a lower triangular matrix. Thus the transposed inverse (N, F̂, Ê) is the same permutation of an upper triangular matrix so that Ĉ is also a congruence system in minimal form. The generating system Ĝ is the homogeneous form of a system that contains a point if and only if there is an element q ∈ Q̂ such that q_{0} = 1 and piv_{>}(q) = 1; that is, if and only if there is also an element a ∈ F̂ such that a_{0} = 1 and piv_{<}(a) = 1. Thus Ĉ is a congruence system in minimal homogeneous form if and only if Ĝ is a generator system in minimal homogeneous form.

More generally, for vectors v_{i} and w_{i} in (N, F̂, Ê) and (L̂, Q̂, M), respectively, piv_{<}(v_{i}) = piv_{>}(w_{i}) = k and v_{ik}w_{ik} = 1. Since # N = # L̂, # Q̂ = # F̂ and # Ê = # M we have proved items (6) and (7).

Since (N, F̂, Ê) = ((Q̂, L̂, M)^{−1})^{T} we have that x̂ ∈ ggen(Ĝ) if and only if N^{T}x̂ = λ, F̂^{T}x̂ = π, Ê^{T}x̂ = 0 for some λ ∈ ℝ^{ℓ}, π ∈ ℤ^{q+1}. Since λ ∈ ℝ^{ℓ}, the condition N^{T}x̂ = λ is vacuous and, as π ∈ ℤ^{q+1}, the condition F̂^{T}x̂ = π holds if and only if F̂^{T}x̂ ≡_{1} 0. Thus we have that x̂ ∈ ggen(Ĝ) if and only if F̂^{T}x̂ ≡_{1} 0, Ê^{T}x̂ = 0. Therefore gcon(Ĉ) = ggen(Ĝ) and (8) holds. ⊓ ⊔

Proof (of Propositions 7 and 8). These are direct consequences of Lemma 1. ⊓ ⊔

Proof (of Proposition 9). Suppose that C is not in strong minimal form. Then, by Definition 5, there exists a proper congruence β = (⟨a, x⟩ ≡_{1} b) ∈ C, such that the following holds:

1. there exists i > 0 and a proper congruence γ = (⟨c, x⟩ ≡_{1} d) ∈ C \{β} where piv_{<}(γ) = i and either 2a_{i} ≤ −c_{i} or 2a_{i} > c_{i}.

Suppose that 0 ≤ k ≤ n is the maximum value for the index i such that condition (1) holds.

We show, by induction on k, that there exists a sequence of at most n transformations, each of which having complexity O(n), from β to the congruence β′ = (⟨a′, x⟩ ≡_{1} b′), such that, if C′ := (C \{β})∪{β′}, then gcon(C′) = gcon(C) and condition (1) (when β is replaced by β′) does not hold.

If k = 0, then condition (1) does not hold for β. Therefore let β′ = β.

Suppose now that k > 0 so that condition (1) holds for i = k. As C is in minimal form, k < piv_{<}(a). Let

a′′ = \begin{cases}a − ⌈\frac{a_{k}}{c_{k}} ⌉ c and b′′ = b − ⌈\frac{a_{k}}{c_{k}} ⌉ d, & if a_{k} mod c_{k} > \frac{a_{k}}{2} ; \\ a − ⌊\frac{a_{k}}{c_{k}} ⌋ c and b′′ = b − ⌊\frac{a_{k}}{c_{k}} ⌋ d, & if a_{k} mod c_{k} ≤ \frac{a_{k}}{2} . \end{cases}

Then −c_{k} < 2a′′_{k} ≤ c_{k}. Also, for k + 1 ≤ j ≤ n, we have c_{j} = 0 so that a_{j} = a′′_{j} and piv_{<}(a′′) = piv_{<}(a). Letting β′′ := (⟨a′′, x⟩ ≡_{1} b′′) and C′′ := (C \ {β}) ∪ {β′′}, we have gcon(C′′) = gcon(C). Note that this transformation has a complexity O(n). As k′′, the maximum index such that condition (1) holds for β′′, is strictly less than k we can apply the inductive hypothesis to C′′ and β′′. Thus there is a sequence of at most n−1 transformations from β′′ to β′ such that, gcon((C′′ \ {β′′}) ∪ {β′}) = gcon(C′′) and condition (1) (when β is replaced by β′) does not hold. Thus there is a sequence of at most n transformations from β to β′ such that gcon((C \{β})∪{β′}) = gcon(C). As each of the individual steps

has complexity O(n), the sequence of transformations has complexity O(n^{2}).

We repeat this sequence of transformations for each proper congruence in C to obtain a congruence system C′ such that, for each proper congruence β′ ∈ C′, condition (1) does not hold. Thus, by Definition 5, C′ is in strong normal form. Thus, as there are at most n proper congruences in C since, by hypothesis, C is in minimal form, the complexity of computing the strong minimal form is O(n^{3}). ⊓ ⊔

The proofs of Propositions 10 and 11 depend on the following lemma. This shows that if one grid is a subset of another then the pivot elements of the proper congruences of the larger grid must be divisible by the corresponding pivot elements of the smaller grid.

Lemma 2. Let L_{1} = gcon(C_{1}), L_{2} = gcon(C_{2}) be non-empty grids in 𝔾_{n} such that L_{1} ⊆ L_{2} and the congruence systems C_{1} and C_{2} are in minimal form. Then, for each γ = (⟨c, x⟩ ≡_{g} d) ∈ C_{2}, there exists β = (⟨a, x⟩ ≡_{f} b) ∈ C_{1} such that piv_{<}(a) = piv_{<}(c) = k and either f = g = 0 or g ≠ 0 and ga_{k} | fc_{k}.

Proof. Suppose γ = (⟨c, x⟩ ≡_{g} d) ∈ C_{2} and k = piv_{<}(c). Then, as L_{1} ⊆ L_{2}, L_{1} ⊆ gcon({γ}). Let G_{1} = (L_{1}, Q_{1}, {p}) be a generator system for L_{1} in minimal form constructed as in Lemma 1 from C_{1}.

We first prove that there exists β = (⟨a, x⟩ ≡_{f} b) ∈ C_{1} such that piv_{<}(a) = k. To see this, suppose instead that, for all β = (⟨a, x⟩ ≡_{f} b) ∈ C_{1}, piv_{<}(a) ≠ k. Then, by Lemma 1, there must exist a line ℓ ∈ L_{1} such that piv_{>}(ℓ) = k; hence ⟨c, ℓ⟩ = c_{k}ℓ_{k} ≠ 0. Since L_{1} ⊆ gcon({γ}), this implies that

⟨c, (p + rℓ)⟩ = ⟨c, p⟩ + rc_{k}ℓ_{k} ≡_{g} d,

for all r ∈ ℝ; which is a contradiction.

We next show that if g = 0 then f = 0. To see this, suppose instead that g = 0 but f ≠ 0. Then, by Lemma 1, there exists q ∈ Q_{1} such that piv_{>}(q) = k; hence ⟨c, q⟩ = c_{k}q_{k} ≠ 0. Since L_{1} ⊆ gcon({γ}), this implies that

⟨c, (p + rq)⟩ = ⟨c, p⟩ + mc_{k}q_{k} ≡_{g} d,

for all m ∈ ℤ; which is a contradiction.

We now assume that g ≠ 0 and show that ga_{k} | fc_{k}. This is trivial if f = 0; therefore, suppose f ≠ 0. By Lemma 1, there exists a parameter q in Q_{1} such that piv_{>}(q) = k (so that q_{k}c_{k} ≠ 0) and q_{k} = fa_{k}^{−1}. Thus, as L_{1} ⊆ gcon({γ}), ⟨q, c⟩ = q_{k}c_{k} = mg, for some m ∈ ℤ\{0}. Therefore we must have ga_{k} | fc_{k}. ⊓ ⊔

Proof (of Proposition 4). Let C′ be the congruence system obtained, as in Lemma 1, from G. Let G = (L, Q, P) and let C = (F, E) and C′ = (F′, E′) where E, E′ are sets of equalities and F, F′ are sets of proper congruences. Then, by Lemma 1,

# Q = # F′ = n − # L − # E′.

By applying Lemma 2 twice where L_{1} = L_{2} = L, we obtain # E = # E′ and # F = # F′. Therefore

# Q = # F = n − # L − # E.
⊓ ⊔

The next lemma, which is needed in the proofs of Propositions 2, 10 and 11 shows that if, two grids, one a subset of the other are described by two congruence systems in strong minimal form that are pivot equivalent, then, relative to the affine hull of the grids, pivot equivalent congruences in these systems are the same.

Lemma 3. Let L_{1} = gcon(C_{1}), L_{2} = gcon(C_{2}) be non-empty grids in 𝔾_{n} where L_{1} ⊆ L_{2} and the congruence systems C_{1} and C_{2} are in strong minimal form. Suppose that C_{1} is pivot equivalent to C_{2}. Then, for each β ∈ C_{1} and γ ∈ C_{2} such that β ⇑ γ,

gcon({β}) ∩ affine.hull(L_{1}) = gcon({γ}) ∩ affine.hull(L_{1}).  (10)

Proof. Let β = (⟨a, x⟩ ≡_{f} b) ∈ C_{1}. By the definition of pivot equivalence for congruence systems in Section 3, as C_{1} ⇑ C_{2} there exists γ = (⟨c, x⟩ ≡_{g} d) ∈ C_{2} such that β ⇑ γ. We show that equation (10) holds. By the definition of pivot equivalence for congruences in Section 2, piv_{<}(a) = piv_{<}(c) = k and ga_{k} = fc_{k}. Thus as a_{k}, c_{k} ≠ 0, either f = g = 0 and β, γ are both equalities, or we have f, g ≠ 0 so that β, γ are both proper congruences and we can assume that f = g = 1.

Let E_{1} be the set of equalities in C_{1}. By Gaussian elimination, the set E_{1} can be transformed to the set of equalities E′_{1} such that gcon(E′_{1}) = gcon(E_{1}) = affine.hull(L_{1}) and has the following property: let E′_{1} = {β_{1}, …, β_{m}} such that, for each i ∈ {1, …, m}, piv_{<}(β_{i}) = k_{i} and β_{i} = (⟨a_{i}, x⟩ = b_{i}); then, for each i, j ∈ {1, …, m} where i ≠ j, we have a_{ik_{j}} = 0. Let

a′′ = a − c − ∑_{i=1}^{m} \frac{(a_{k_{i}} − c_{k_{i}})}{a_{k_{i}}}a_{i},  b′′ = b − d − ∑_{i=1}^{m} \frac{(a_{k_{i}} − c_{k_{i}})}{a_{k_{i}}}b_{i}. (11)

Let β′′ := (⟨a′′, x⟩ ≡_{1} b′′) and ℓ := piv_{<}(a′′). Then L_{1} ⊆ gcon({β′′}) ⊆ L_{2}. Moreover, for any equality β_{i} ∈ E′_{1}, piv_{<}(a_{i}) ≠ ℓ. Thus, if β, γ are equalities, a′′ = 0 and, as C_{1} is consistent, b′′ = 0. Therefore gcon(E′_{1}) ⊆ gcon({β}) and gcon(E′_{1}) ⊆ gcon({γ}). Hence equation (10) holds.

Consider now the case when β, γ are proper congruences. We first show that a′′ = 0 and b′′ ∈ ℤ. Without loss of generality we can assume that f = g = 1. Note that, as a_{k} = c_{k} we have ℓ < k. We show ℓ = 0; suppose, to the contrary that ℓ > 0. Since L_{1} ⊆ L_{2}, and γ ∈ C_{2}, we have L_{1} ⊆ gcon({γ}); so we can apply Lemma 2 to the grids L_{1} and gcon({γ}). Thus there exists a proper congruence β′ = (⟨a′, x⟩ ≡_{1} b′) ∈ C_{1} where piv_{<}(a′) = ℓ and a′_{ℓ} | a′′_{ℓ}. Note that the number of proper congruences p_{1} in C_{1} is equal to the number of proper congruences p_{2} in C_{2}; since by Lemma 2, p_{2} ≤ p_{1} and, by hypothesis, p_{1} ≤ p_{2}. Therefore, by Lemma 2, there must exist a proper congruence γ′ = (⟨c′, x⟩ ≡_{1} d′) ∈ C_{2} where piv_{<}(c′) = ℓ and c′_{ℓ} = a′_{ℓ}. Now as C_{1} and C_{2} are in strong minimal form, by Definition 5,

−\frac{a′_{ℓ}}{2} < a_{ℓ} ≤ \frac{a′_{ℓ}}{2} and − \frac{c′_{ℓ}}{2} < c_{ℓ} ≤ \frac{c′_{ℓ}}{2} .

Therefore −a′_{ℓ} < a′′_{ℓ} < a′_{ℓ}. It follows that, as a′_{ℓ}|a′′_{ℓ}, a′′_{ℓ} = 0, contradicting the assumption that piv_{<}(a′′) = ℓ > 0. Therefore a′′ = 0 and β′′ is the relation b′′ ≡_{1} 0 for some b′′ ∈ ℤ.

It follows that, by (11),

a − c = ∑_{i=1}^{m} \frac{(a_{k_{i}} − c_{k_{i}})}{a_{k_{i}}}a_{i}, b − d ≡_{1} ∑_{i=1}^{m} \frac{(a_{k_{i}} − c_{k_{i}})}{a_{k_{i}}}b_{i}.

Thus

gcon({γ, β_{1}, …, β_{m}}) ⊆ gcon({β}) gcon({β, β_{1}, …, β_{m}}) ⊆ gcon({γ})

so that

gcon({γ, β_{1}, …, β_{m}}) = gcon({β, β_{1}, …, β_{m}}).

Hence equation (10) holds. ⊓ ⊔

Proof (of Proposition 2). By Proposition 9, we can convert C_{1} and C_{2} to strong minimal form, C′_{1} and C′_{2} respectively, so that, for i = 1, 2, L_{i} = gcon(C′_{i}) and C′_{i} is pivot equivalent to C_{i}. By Lemma 3, for each β ∈ C′_{1} and γ ∈ C′_{2},

gcon({β}) ∩ affine.hull(L_{1}) = gcon({γ}) ∩ affine.hull(L_{1}).

Thus L_{1} = L_{2} as required. ⊓ ⊔

Proof (of Proposition 10). Suppose C_{1}, C_{2} ≠ ∅. By Lemma 2, for each β = (⟨a, x⟩ ≡_{f} b) ∈ C_{1} there exists γ = (⟨c, x⟩ ≡_{g} d) ∈ C_{2} such that piv_{<}(a) = piv_{<}(c) = k and either f = g = 0 or f ≠ 0 and fc_{k} | ga_{k}.

Also by Lemma 2, for γ = (⟨c, x⟩ ≡_{g} d) ∈ C_{2} there exists β′ = (⟨a′, x⟩ ≡_{f′} b′) ∈ C_{1}, such that piv_{<}(a′) = piv_{<}(c) = k and either f′ = g = 0 or g ≠ 0 and ga′_{k} | f′c_{k}. However as C_{1} is in strong minimal form β = β′. Therefore, by Definition 2, a_{k}, c_{k} > 0, hence fc_{k} = ga_{k}.

Therefore, for each β ∈ C_{1} there exists γ ∈ C_{2} such that β ⇑ γ so that we can apply Lemma 3 to prove the thesis. ⊓ ⊔

Proof (of Proposition 11). In order to show that ‘∇’ is a widening operator, we prove that conditions (1) and (2) in Definition 4 hold. Let L_{1} = gcon(C_{1}), L_{2} = gcon(C_{2}) ∈ 𝔾_{n}, where L_{1} ⊆ L_{2}, C_{1} is in minimal form and C_{2} is in strong minimal form.

By Definition 6, if L_{1} = ∅ or dim(L_{1}) < dim(L_{2}), then L_{1} ∇ L_{2} = L_{2}. Therefore, in this case, condition (1) holds. Clearly, the empty grid can occur only as the first element of a strictly increasing chain of grids; moreover, if L and L′ are any two successive and distinct grids in the increasing chain of condition (2) in Definition 4, then 0 ≤ dim(L) ≤ dim(L′) ≤ n. Hence, the case when L_{1} = ∅ or dim(L_{1}) < dim(L_{2}) hold can occur no more than a finite number of times in such a chain.

Suppose now that L_{1} ≠ ∅ and dim(L_{1}) = dim(L_{2}), so that the second case of the widening computation applies (note that, due to the inclusion hypothesis, dim(L_{1}) > dim(L_{2}) cannot hold), and let C_{S} be as given in equation (6) in Definition 6. Then, since C_{S} ⊆ C_{2}, condition (1) holds. By Proposition 2, if C_{S} = C_{2}, we have L_{1} = L_{2}; thus, if L_{1} ≠ L_{2}, we have # C_{S} < # C_{2}. By Lemma 2, as C_{1} and C_{2} are in minimal form, it follows that # C_{2} ≤ # C_{1} so that, if L_{1} ≠ L_{2}, # C_{S} < # C_{1}. Therefore condition (2) of Definition 4 holds. ⊓ ⊔