Let K be a polyhedron in Rd, given by a system of m linear inequalities, with rational number coe cients bounded over in absolute value by L. We propose an algorithm for computing an irredundant representation of the integer points of K, in terms of \simpler" polyhedra, each of them having at least one integer point. Using the terminology of W. Pugh: for any such polyhedron P , no integer point of its grey shadow extends to an integer point of P . We show that, under mild assumptions, our algorithm runs in exponential time w.r.t. d and in polynomial w.r.t m and L. We report on a software experimentation.
The integer points of polyhedral sets are of interest in many areas of
mathematical sciences, see for instance the landmark textbooks of A. Schrijver [
In [
Projections of polyhedra and parametric programming are tightly related
problems, see [
In [
Our algorithm takes as input a system of linear inequalities Ax ≤ b where A is a matrix over Z with m rows and d columns, x is the unknown vector and b is a vector of m coe cients in Z. The points x ∈ Rd satisfying Ax ≤ b form a polyhedron K and our algorithm decomposes its integer points (that is, K ∩ Zd) into a disjoint union K1 ∩ Zd ⊍ ⋯ ⊍ Ke ∩ Zd, where K1; : : : ; Ke are \simpler" polyhedra such that Ki ∩ Zd ≠ ∅ holds, for 1 ≤ i ≤ e. To use the terminology introduced by W. Pugh for the Omega test, for 1 ≤ i ≤ e, no integer point of the grey shadow of the polyhedron Ki extends to an integer point of Ki. As a consequence, applying our algorithm to Ki would return Ki itself, for 1 ≤ i ≤ e. Let us present the key principles and features of our algorithm through an example. Consider the polyhedron K of R4 given below: ⎧⎪ 2x + 3y − 4z + 3w ≤ 1 ⎪ ⎪⎪⎪⎪⎪ −2x − 3y + 4z − 3w ≤ −1 ⎪ ⎪⎪⎪⎪⎪−13x − 18y + 24z − 20w ≤ −1 : ⎪⎪⎨−26x − 40y + 54z − 39w ≤ 0 ⎪ ⎪ ⎪⎪⎪⎪−24x − 38y + 49z − 31w ≤ 5 ⎪ ⎪ ⎪⎪⎪ 54x + 81y − 109z + 81w ≤ 2 ⎩ A rst procedure, called IntegerNormalize, detect implicit equations and solve them using techniques based on Hermite normal form, see Sect. 3. In our example 2x+3y−4z+3w = 1 is an implicit equation and IntegerNormalize(Ax ≤ b) returns a triple (t; x = Pt+ q; Mt ≤ v) where t is a new unknown vector, the linear system x = Pt+q gives the general form of an integer solution of the implicit equation(s) and Mt ≤ v is obtained by substituting x = Pt + q into Ax ≤ b. In our example, the systems x = Pt + q and Mt ≤ v are given by: ⎪⎪⎧⎪ x = −3t1 + 2t2 − 3t3 + 2 ⎪⎪⎧⎪ 3t1 − 2t2 + t3 ≤ 7 ⎪⎪⎪⎪ y = 2t1 + t3 − 1 and ⎪⎪⎪⎪−2t1 + 2t2 − t3 ≤ 12 : ⎪⎪⎪⎩⎪⎨⎪⎪⎪wz == tt32 ⎩⎪⎪⎪⎪⎪⎪⎪⎨−4t1 + t2 + −3tt23 ≤≤ 1−525 A second procedure, called DarkShadow, takes Mt ≤ v as input and returns a couple (t′; ) where t′ stands for all t-variables except t1, and is a linear system in the t′-variables such that any integer point solving of extends to an integer point solving Mt ≤ v. In our example, t′ = {t2; t3} and is given by: The polyhedron D of R2 de ned by , and the inequalities of Mt ≤ v not involving t1, is called the dark shadow of the polyhedron de ned by Mt ≤ v. On the left-hand side of Fig. 1, one can see the polyhedron de ned in R3 by Mt ≤ v together with its dark shadow D (shown in dark blue) as well as its projection on the (t2; t3)-plane, denoted by R and called real shadow by W. Pugh. The righthand side of Fig. 1 gives a planar view of D and R. It turns out that, if M′t′ ≤ v′ is the linear system generated by applying Fourier-Motzkin elimination (without removing redundant inequalities) to Mt ≤ v (in order to eliminate t1), then is given by a linear system of the form M′t′ ≤ w′. This explains why, on the right-hand side of Fig. 1, each facet of the dark shadow D is parallel to a facet of the real shadow R. While this property is observed on almost all practical problems, in particular in the area of analysis and transformation of computer programs, it is possible to build examples where this property does not hold.
On the right-hand side of Fig. 1, one observes that the region R ∖ D, called grey shadow, contains integer points. Some of them, like (t2; t3) = (29; 9), do not extend to an integer solution of Mt ≤ v. Indeed, plugging (t2; t3) = (29; 9) into Mt ≤ v yields 327 ≤ t1 ≤ 3
56 , which has no integer solutions. However, other integer points of R ∖ D may extend to integer solutions of Mt ≤ v. In order to determine them, a third procedure, called GreyShadow, is needed. We give a disjoint decomposition of grey shadow by negating each inequality of in turn. Then, compute the integer points in each disjoint part. Details are given in Sect. 4.
Returning to our example, the negation of the inequality 2t2 − t3 ≤ 48 from , combined with the system Mt ≤ v, yields the following
⎪⎪⎪⎪⎩⎨⎪⎪⎪⎧⎪ttt321 === tt−542+t41+ 2t5 + 1 ; and ⎪⎪⎨⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎪⎪⎧−1−0t24t4+−−7tttt5554 ≤≤≤≤ 18−−12448 ;
where t4; t5 are new variables. Continuing in this manner with the GreyShadow
procedure, a decomposition of the integer points of Mt ≤ v is given by:
⎧⎪ 3t1 − 2t2 + t3 ≤ 7
⎪
⎪
⎪⎪⎪⎪−2t1 + 2t2 − t3 ≤ 12
⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪−4t−15+t2t22t+2+1−−33tttt3332 ≤≤≤≤ 146−58725 ; ⎪⎪⎪⎪⎩⎪⎨⎪⎧⎪⎪ttt321 === 121675 ; ⎪⎪⎪⎪⎩⎪⎨⎪⎧⎪⎪ttt321 === 131838 ; ⎪⎪⎪⎪⎩⎪⎨⎪⎧⎪⎪ttt213 === 121554 ; ⎪⎪⎪⎪⎪⎪⎪⎩⎪⎨⎪⎪⎧⎪⎪⎪⎪−25ttt≤321t===6 551≤009−++16t26:t6 :
⎪
⎩⎪⎪⎪ 2 ≤ t3 ≤ 17
Denoting these ve systems respectively by S1; : : : ; S5 the integer points of K are
nally given by the union of the integer points of the systems x = Pt + q ∪ Si,
for 1 ≤ i ≤ 5. The systems S2; : : : ; S5 look simple enough to be considered as
solution sets. What about S1? The system S1, as well as S2; : : : ; S5, satis es a
\back-substitution" property which is similar to that of a regular chain in the
theory of polynomial system solving [
With respect to the original Omega Test [
Omega Test may search for integer points in intersecting polyhedral regions. 3. The original Omega Test uses an ad-hoc routine for computing the integer solutions of linear equation systems, while we rely on Hermite normal form for this task. Consequently, we deduce complexity estimates for that task. 4. We also provide complexity estimates for the procedures GreySahdow and DarkShadow under realistic assumptions. From there, we derive complexity estimates for the entire algorithm, whereas no complexity estimates were known for the original Omega Test. Our algorithm is discussed in Sect. 3 and 4 while complexity estimates are stated in Sect. 5. Sect. 2 gathers background materials on polyhedra. 2
This section is a review of the theory of polyhedral sets. It is based on the books
of B. Grunbaum [
Given a positive integer d, we consider the d-dimensional Euclidean space Rd equipped with the Euclidean topology. Let K be a subset of Rd. The dimension dim(K) of K is a − 1 where a is the maximum number of a nely independent pointsxin∈ KR.dLSetaTax∈ =Rbd, .leWteb ∈saRy athnadt dtehneotheypbeyrpHlantheeHhyspueprppolratnse Kdeifneeidthbeyr H = { } sup{aT x S x ∈ K} = b or inf{aT x S x ∈ K} = b holds, but not both.
From now on, let us assume that K is convex. A set F ⊆ K is a face if either F = ∅ or F = K, or if there exists a hyperplane H supporting K such that we have F = K ∩ H. The set of all faces of K is denoted by F (K). We say that F ∈ F (K) is proper if we have F ≠ ∅ or F ≠ K. We note that the intersection of any family of faces of K is itself a face of K.
We say that K is a polyhedral set or polyhedron if it is the intersection of nitely many closed half-spaces of Rd. We say that K is full-dimensional, if we have dim(K) = d, that is, if the interior of K is not empty. The proper faces of K that are ⊆-maximal are called facets and those of dimension zero are called vertices. We observe that every face of K is also a polyhedral set.
Let H1; : : : ; Hm be closed half-spaces such that the intersection ∩ii==1m Hi is irredundant, that is, ∩ii==1m Hi ≠ ∩i==1m;j≠i Hi for all 1 ≤ j ≤ m. We observe that this i intersection is closed and convex. For each i, where 1 ≤ i ≤ m, let ai ∈ Rd and bi ∈ R such that Hi is de ned by aiT x ≤ bi. We denote by A the m × d matrix (aiT ; 1 ≤ i ≤ m) and by b the vector (b1; : : : ; bd)T .
From now on, we assume that K = ∩ii==1m Hi holds. Such irredundant
decomposition of a polyhedral set can be computed from an arbitrary intersection of
nitely many closed half-spaces, in time polynomial in both d and m, using linear
programming; see L. Khachian in [
For a ∈ Rd and b ∈ R, we say that aT x ≤ b is an implicit equation in Ax ≤ b if for all x ∈ Rd we have
Ax ≤ b
Ô⇒
aT x = b :
(1)
Following [
A=x = b=
and A+x < b+ : The facets of K are in 1-to-1 correspondence with the inequalities of A+x ≤ b+. In addition, if K is full-dimensional, then A+ = A and b+ = b both hold; moreover the system of inequalities Ax ≤ b is a unique representation of K, up to multiplication of inequalities by positive scalars.
From now on and in the sequel of this paper, we assume that variables are
ordered as x1 > ⋯ > xd. We call initial coe cient, or simply initial, of an
inequality aiT x ≤ bi, for 1 ≤ i ≤ m, the coe cient of aiT x in its largest variable.
Following the terminology of W. Pugh in [
Canonical representation. Recall that we assume that none of the inequalities of Ax ≤ b is redundant. If K is full-dimensional and if the initial of each inequality in Ax ≤ b is 1 or −1, then we call Ax ≤ b the canonical representation of K w.r.t. the variable ordering x1 > ⋯ > xd and we denote it by can(K; x1; : : : ; xd).
We observe that the notion of canonical representation can also be expressed in a more geometrical and less algebraic way, that is, independently of any coordinate system. Assume again that K is full-dimensional and that the intersection ∩ii==1n Hi = K of closed half-spaces H1; : : : ; Hn is irredundant. Since K is full-dimensional, the supporting hyperplane of each facet of K must be the frontier of one half-space among H1; : : : ; Hn. Clearly, two (or more) half-spaces among H1; : : : ; Hn may not have the same frontier without contradicting one of our hypotheses (K is full-dimensional, ∩ii==1n Hi is irredundant). Therefore, the half-spaces H1; : : : ; Hn are in one-to-one correspondence with the facets of K. This implies that there is a unique irredundant intersection of closed half-spaces equaling K and we denote it by can(K).
Projected representation. Let again Ax ≤ b be a canonical representation of
the polyhedral set K w.r.t. the variable ordering x1 > ⋯ > xd. We denote by
Ax1 (resp. A<x1 ) and bx1 (resp. b<x1 ) the rows of A and b corresponding to
the inequalities whose largest variable is x1 (resp. less than x1). For each upper
bound cx1 ≤ of x1 and each lower bound −ax1 ≤ − of x1 (where c > 0,
a > 0, ∈ R[x2; : : : ; xd] and ∈ R[x2; : : : ; xd] hold), we have a new inequality
c −a ≤ 0. Augmenting A<x1 with all inequalities obtained in this way, we obtain
a new linear system which represents a polyhedral set which is the standard
projection of K on the d − 1 least coordinates of Rd, namely (x2; : : : ; xd); hence
we denote this latter polyhedral set by x2;:::;xd K and we call it the real shadow
of K, following the terminology of [
We review how Hermite normal forms [
We say that H is the column Hermite normal form of A if H is a column Hermite normal form and there exists a uni-modular d × d-matrix U over Z such that we have H = AU . When those properties hold, we call {f (d − r + 1); : : : ; f (d)} the pivot row set of A.
Remark 1. The matrix A admits a unique column Hermite normal form. Let H be this column Hermite normal form and let U be the uni-modular d × d-matrix given in De nition 1. Let us decompose U as U = [UL; UR] where UL(resp. UR) consist of the rst d − r (resp. last r) columns of U . Then we de ne HL ∶= AUL and HR ∶= AUR. We have HL = 0m;d−r, where 0m;d−r is the zero-matrix with m rows and d − r columns. We observe that UR is a full column-rank matrix. Moreover, if A is full row-rank, that is, if r = m holds, then UR is non-singular.
Lemma 2 shows how to compute the integer solutions of the system of linear
equations Ax = b when A is full rank. In the general case, one can use Lemma 1
to reduce to the hypothesis of Lemma 2. While the construction of this latter
lemma relies on the HNF, alternative approaches are available. For instance, one
can use the equation elimination procedure of the Omega Test [
Notation 1 For I ⊆ {1; : : : ; m}, we denote by AI (resp. bI ) the sub-matrix (resp. vector) of A (resp. b) consisting of the rows of A (coe cients of b) with indices in I.
Lemma 1. Let I be the pivot row set of A, as given in De nition 1. Assume that Ax = b admits at least one solution in Rd. Then, for any x ∈ Rd, we have Ax = b ⇐⇒
AI x = bI : Lemma 2. We use the same notations as in De nition 1 and Remark 1. We assume that HR is non-singular. Then, the system Ax = b has an integer solution if and only if HR−1b is integral. In this case, all integral solutions to Ax = b are given by x = P t + q where 1. the columns of P consist of a Z-basis of the space {x ∶ Ax = 0}, 2. q is a particular solution of Ax = b, and 3. t = (t1; : : : ; td−r) is a vector of d − r unknowns.
The maximum absolute value of any coe cient in P (resp. q) can be bounded over by rr+1L2r (resp. rr+1L2r), where L is the maximum absolute value of a coe cient in A (resp. in either A or b). Moreover, P and q can be computed within O(mdr2(log r + log L)2 + r4(log r + log L)3) bit operations.
The matrix H is the column HNF of A, with unimodular matrix U and pivot row set [2; 4; 5]. We denote by HR the sub-matrix of H whose coe cients are in bold fonts. Applying Lemma 1, we deduce that for any vector b such that Ax = b admits one rational solution, we have: ⎧⎪ 3x1 + 4x2 − 4x3 − x4 = b1 ⎪ ⎪⎪⎪⎪⎪⎪2x1 − 2x2 + 8x3 + 4x4 = b2 ⎪⎧⎪2x1 − 2x2 + 8x3 + 4x4 = b2 ⎪⎨⎪5x1 + 2x2 + 4x3 + 3x4 = b3 ⇔ ⎨⎪⎪3x1 + 5x2 − 5x3 − 2x4 = b4 : ⎪⎪⎪⎪⎪⎪⎪3x1 + 5x2 − 5x3 − 2x4 = b4 ⎪⎪⎪⎪⎩2x1 − 3x2 + 9x3 + 5x4 = b5 ⎪⎪⎪2x1 − 3x2 + 9x3 + 5x4 = b5 ⎩ We apply Lemma 2: if Ax = b is consistent over Q and if HR−1[b2; b4; b5]T is integral, then all the integer solutions of the second equation system in Relation (2) are given by x = Pt + q, where P = [−1; 1; 0; 1]T , q = [ 53 b2 − 139 b4 − 1535 b5; − 1387 b2 + 793 b4 + 5795 b5; − 1189 b2 + 397 b4 + 2996 b5]T , t = (t1) and t1 is a new variable. 4
In this section, we present an algorithm for computing the integer points of a
polyhedron K ⊆ Rd, that is, the set K ∩ Zd. As mentioned in the introduction,
we adapt the Omega Test invented by W. Pugh [
With the notations of Section 2, consider again the polyhedron K: (2) ⎧ a11x1 + ⋯ + a1dxd ≤ b1 ⎪ ⎪ ⎪⎨ ⋮ ⎪⎪⎪ am1x1 + ⋯ + amdxd ≤ bm ⎩ We aim at eliminating x1, so we consider two inequalities in x1: { one with ai1 > 0, called upper bound:
li ∶ ai1x1 + ⋯ + aidxd ≤ bi, { another one with aj1 < 0, called lower bound:
lj ∶ aj1x1 + ⋯ + ajdxd ≤ bj.
The real projection makes a linear combination rij of li and lj:
− aj1(ai2x2 + ⋯ + aidxd) + ai1(aj2x2 + ⋯ + ajdxd) ≤ −aj1bi + ai1bj. The dark projection translates rij towards the center of K into dsij: −aj1(ai2x2+⋯+aidxd)+ai1(aj2x2+⋯+ajdxd) ≤ −aj1bi+ai1bj−(ai1 − 1)(−aj1 − 1). Lemma 3. For any integer point (x2; : : : ; xd) ∈ Zd−1 satisfying dsij, there exists at least one integer x1 ∈ Z satisfying both li and lj.
Let S<x1 be the subset of the inequalities de ning K in which x1 does not occur. Combining all the real projections with S<x1 , we obtain the so-called real shadow of the polyhedron K w:r:t: x1, that we denote by R. Meanwhile, combining all the dark projections with S<x1 , we obtain the so-called dark shadow of the polyhedron K w:r:t: x1, that we denote by D. Without loss of generality, we assume D = S<x1 ∩ ds1 ∩ ⋯ ∩ dst, where t is a positive integer and each dsi is a dark projection introduced above for 1 ≤ i ≤ t. Obviously, we have D = R ∩ ds1 ∩ ⋯ ∩ dst since D ⊂ R ⊂ S<x1 .
We de ne the grey shadow of K w.r.t. x1 as G ∶= R ∖ D. Then, it is easy to deduce that G = Gi, where " means disjoint union and Gi = R ∩ ds1 ∩ 1 "i t ≤ ≤ ⋯ ∩ dsi−1 ∩ ¬dsi, here ¬dsi is the negation of dsi for 1 ≤ i ≤ t. We call each Gi obtained by this way a grey shadow part of K.
We use the following steps to introduce the ideas of procedure GreyShadow: 1. let result ∶= {}; 2. for each lower bound aj1x1 + ⋯ + ajdxd ≤ bj of x1 in K1 3. for each upper bound ai1x1 + ⋯ + aidxd ≤ bi of x1 in K1 4. for each integer ` from 1 to ⌊(−aj1ai1 + aj1 − ai1)~ai1⌋ do 5. let gji` ∶ aj1x1 + ⋯ + ajdxd = bj − `; 6. solving for the integer solutions of gji`, we obtain a solution of the form x = Pt + q; 7. substitute x = Pt + q into K ∩ ¬dsji, we obtain a new polyhedron K1; 8. let result = result ∪ {x = Pt + q ∩ K1}; 9. end for; 10. let K = K ∩ dsji; 11. end for; 12. end for; 13. return result.
Example 2. Given a polyhedron K:
, ⎪⎪⎨−4x1 + x2 + 3x3 ≤ 15 ⎪ ⎪ ⎩⎪⎪⎪ −x2 ≤ −25 we have one upper bound and two lower bounds for x1. For the lower bound −2x1 + 2x2 − x3 ≤ 12 and the upper bound 3x1 − 2x2 + x3 ≤ 7, we have ⌊(−aj1ai1 + aj1 − ai1)~ai1⌋ = 0. Thus, one recursive call with: ⎧⎪−2x1 + 2x2 − x3 = 12 ⎪ ⎪ ⎪⎪⎪ 3x1 − 2x2 + x3 ≤ 7 ⎪ . ⎪⎪⎨−4x1 + x2 + 3x3 ≤ 15 ⎪ ⎪ ⎪⎪⎩⎪ −x2 ≤ −25 Solving −2x1 + 2x2 − x3 = 12 via Hermite normal form leads to: ⎧⎪⎪x1 = x4 ⎪ ⎪⎨x2 = x5 + 1 ⎪ ⎪⎪⎪⎩x3 = −2x4 + 2x5 + 1 , where x4; x5 are new variables. Thus, the grey shadow part G20 is: ⎪⎧⎪ x4 ≤ 8 ⎪ ⎪⎨−10x4 + 7x5 ≤ 11 . ⎪ ⎪⎪⎩⎪ −x5 ≤ −24 Befor compute the other grey shadow parts generated by any other lower bounds, we add the dark projection 2x2 − x3 ≤ 48 to the polyhedron K. For the lower bound −4x1 + x2 + 3x3 ≤ 15, we have ⌊(−aj1ai1 + aj1 − ai1)~ai1⌋ = 1, yielding two other grey shadow parts, which are disjoint with G20. 5
Notation 2 Fig. 2 illustrates the tree of recursive calls for the main algorithm, the IntegerSolve procedure. The root of the tree is labelled with S, which stands for the input system. The left (resp. right) child of a node, other than a leaf, is labelled by D (resp. G) which stands for the output of the DarkShadow procedure (resp. the GreyShadow procedure). Since the DarkShadow procedure generates one input, we use a simple → arrow as an edge to a D-node. However, the GreyShadow procedure may generate several inputs, leading to several recursive calls. Thus, we use a ⇒ arrow as an edge to a G-node. The numbers on the right-hand side of Fig. 2 stand for the levels in the tree.
Hypothesis 1 We assume that during the execution of the function call IntegerSolve(K), for any polyhedral set K′, input of a recursive call, each facet of the dark shadow of K′ is parallel to a facet of the real shadow of K′.
D G
d − 1
To state our main result, we need a notation for the running time of solving
a linear program. Indeed, linear programming is an essential tool for removing
redundant inequalities generated by Fourier-Motzkin elimination, see [
Lemma 4. Let 0 ≤ k < d be an integer. Let fd;m;k be the number of k-dimensional faces of K. Then, we have m fd;m;k ≤ d − k: Therefore, we have
fd;m;0 + fd;m;1 + ⋯ + fd;m;d−1 ≤ md: Lemma 5. Any path in Fig. 2 can be implemented within O(m2r+2d3+"r10(log d+ log L)3) + O(rm2r+2LP(d; dmrr3(log d + log L))) bit operations.
Lemma 6. We have: Tr+1 ≤ mr+1drr2r2 L3r2 Tr for 0 ≤ r ≤ d − 2. Thus, we have Td−1 ≤ md2 d3d3 L3d3 .
Theorem 1. Under Hypothesis 1, the call function IntegerSolve(K) runs within O(m2d2 d4d3 L4d3 LP(d; mdd4(log d + log L))) bit operations. 6
We have implemented the algorithm presented in the rst paper within the Polyhedra library in Maple. This library is publicly available in source on the download page of the RegularChains library at www.regularchains.org
In this section, we gather experimental results obtained with our implementations. Our objectives were two-fold: Example