Computing the Integer Points of a Polyhedron (Extended Abstract) Rui-Juan Jing1,2 and Marc Moreno Maza2 1 KLMM, UCAS, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, rjing8@uwo.ca, 2 University of Western Ontario, moreno@csd.uwo.ca. Abstract. Let K be a polyhedron in Rd , given by a system of m linear inequalities, with rational number coefficients bounded over in absolute value by L. We propose an algorithm for computing an irredundant rep- resentation 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. 1 Introduction The integer points of polyhedral sets are of interest in many areas of mathemat- ical sciences, see for instance the landmark textbooks of A. Schrijver [19] and A. Barvinok [3], as well as the compilation of articles [4]. One of these areas is the analysis and transformation of computer programs. For instance, integer programming [7] is used by P. Feautrier in the scheduling of for-loop nests [8] Barvinok’s algorithm [2] for counting integer points in polyhedra is adapted by M. Köppe and S. Verdoolaege in [16] to answer questions like how many memory locations are touched by a for-loop nest. In [17], W. Pugh proposes an algorithm, called the Omega Test, for testing whether a polyhedron has integer points. In the same paper, W. Pugh shows how to use the Omega Test for performing de- pendence analysis [17] in for-loop nests. Then, in [18], he uses the Omega Test for deciding Presburger arithmetic formulas. In [18], W. Pugh also suggests, without stating a formal algorithm, that the Omega Test could be used for quantifier elimination on Presburger formulas. This observation is a first motivation for the work presented in this paper: we adapt the Omega Test so as to describe the integer points of a polyhedron via a projection scheme, thus performing elimination of existential quantifiers on Presburger formulas. Projections of polyhedra and parametric programming are tightly related problems, see [12]. Since the latter is essential to the parallelization of for-loop nests [7], which is of interest to the authors [5], we had here a second motivation for developing the proposed algorithm. 2 In [9], M. J. Fischer and M. O. Rabin show that any algorithm for decid- ing Presburger arithmetic formulas has a worst case running time which is a doubly exponential in the length of the input formula. However, this worst case scenario is based on a formula alternating existential and universal quantifiers. Meanwhile, in practice, the original Omega Test (for testing whether a polyhe- dron has integer points) can solve “difficult problems” as shown by W. Pugh in [18] and others, e.g. D. Wonnacott in [20]. This observation brings our third motivation: determining realistic assumptions under which our algorithm, based on the Omega Test, could run in a single exponential time. 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 coefficients 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 first 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 ⎪ ⎪ ⎪−2t1 + 2t2 − t3 ≤ 12 ⎨ and ⎨ . ⎪ ⎪ ⎪ z = t2 ⎪ ⎪ ⎪−4t1 + t2 + 3t3 ≤ 15 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩w = t3 ⎪ ⎩ −t2 ≤ −25 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: 3 ⎧ ⎪ 2t2 − t3 ≤ 48 ⎪ ⎪ ⎪ ⎨−5t2 + 13t3 ≤ 67 ⎪ ⎪ ⎪ ⎪ ⎩ −t2 ≤ −25 The polyhedron D of R2 defined by Θ, and the inequalities of Mt ≤ v not involving t1 , is called the dark shadow of the polyhedron defined by Mt ≤ v. On Fig. 1: The real, the dark and the grey shadows of a polyhedron. the left-hand side of Fig. 1, one can see the polyhedron defined 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 right- hand 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 37 2 ≤ t1 ≤ 56 3 , 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 4 ⎧ ⎪−2t1 + 2t2 − t3 = 12 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 3t1 − 2t2 + t3 ≤ 7 ⎪ ⎪ ⎪ ⎨−4t1 + t2 + 3t3 ≤ 15 , ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ −t2 ≤ −25 ⎪ ⎪ ⎪ ⎪ ⎩ −2t 2 t3 ≤ −49 + which, by means of IntegerNormalize, rewrites to: ⎧ ⎪ t4 ≤ 8 ⎧ ⎪t = t ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 1 4 ⎪ ⎪ ⎪−10t4 + 7t5 ≤ 11 ⎨t2 = t5 + 1 , and ⎨ , ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ −t5 ≤ −24 ⎪ ⎩t3 = −2t4 + 2t5 + 1 ⎪ ⎪ ⎪ ⎪ ⎩ −2t4 − t5 ≤ −48 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 ⎪ ⎪ ⎪ ⎧ ⎪ t1 = 19 ⎪ ⎪ ⎪−4t + t + 3t ≤ 15 ⎧ ⎪t = 15 ⎧ ⎪t = 18 ⎧ ⎪t = 14 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 1 2 3 ⎪ ⎪ ⎪ 1 ⎪ ⎪ ⎪ 1 ⎪ ⎪ ⎪ 1 ⎪ ⎪ ⎪ t2 = 50 + t6 ⎨ 2t2 − t3 ≤ 48 , ⎨t2 = 27 , ⎨t2 = 33 , ⎨t2 = 25 , ⎨ . ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ t3 = 50 + 2t6 ⎪ ⎪ ⎪ −5t + 13t ≤ 67 ⎪ ⎩3 t = 16 ⎪ ⎩3 t = 18 ⎪ ⎩3 t = 15 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 2 3 ⎪ ⎩−25 ≤t6 ≤ −16. ⎪ ⎪ ⎪ −t ≤ −25 ⎪ ⎪ ⎪ 2 ⎪ ⎩ 2 ≤ t3 ≤ 17 Denoting these five systems respectively by S1 , . . . , S5 the integer points of K are finally 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 , satisfies a “back-substitution” property which is similar to that of a regular chain in the theory of polynomial system solving [1]. This property, when applied to S1 , says that for all 1 ≤ i ≤ 2, every integer point of Ri solving all the inequalities of S1 involving t1 , . . . , ti only, extends to an integer point of Ri+1 solving all the inequalities of S1 involving t1 , . . . , ti+1 . With respect to the original Omega Test [17], our contributions are as follows. 1. We turn the decision procedure of the Omega Test into an algorithm decom- posing all the integer points of a polyhedron. 2. Our decomposition is disjoint whereas the recursive calls in the original 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. 5 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 Polyhedral Sets This section is a review of the theory of polyhedral sets. It is based on the books of B. Grünbaum [10] and A. Schrijver [19], where proofs of the statements below can be found. 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 affinely independent points in K. Let a ∈ Rd , let b ∈ R and denote by H the hyperplane defined by H = {x ∈ Rd ∣ aT x = b}. We say that the hyperplane H supports K if either sup{aT x ∣ x ∈ K} = b or inf{aT x ∣ 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 finitely 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 ∩i=mi=1 Hi is irredundant, that is, ∩i=m i=1 Hi ≠ ∩i=1,j≠i Hi for all 1 ≤ j ≤ m. We observe that this i=m intersection is closed and convex. For each i, where 1 ≤ i ≤ m, let ai ∈ Rd and bi ∈ R such that Hi is defined by aTi x ≤ bi . We denote by A the m × d matrix (aTi , 1 ≤ i ≤ m) and by b the vector (b1 , . . . , bd )T . From now on, we assume that K = ∩i=m i=1 Hi holds. Such irredundant decom- position of a polyhedral set can be computed from an arbitrary intersection of finitely many closed half-spaces, in time polynomial in both d and m, using linear programming; see L. Khachian in [15]. The following property is essential. For every face F of K, there exists a subset I of {1, . . . , m} such that F corresponds to the set of solutions to the system of equations and inequalities aTi x = bi for i ∈ I, and aTi x ≤ bi for i ∈/ I . This latter property has several important consequences. For each i 1 ≤ i ≤ m, the set Fi = K ∩ Hi is a facet of K and the border of K equals ∪i=m i=1 Fi . In particular, each proper face of K is contained in a facet of K. Each facet of a facet of K is the intersection of two facets of K. Moreover, if the m × d-matrix A has rank d (full column rank) then the ⊆-minimal faces are the vertices. The set F(K) is finite and has at most 2m elements 6 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 [19], we denote by A= (resp. A+ ) and b= (resp. b+ ) the rows of A and b corresponding to the implicit equations (resp. inequalities). The following properties are easy to prove. If K is not empty, then there exists x ∈ K satisfying both 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 coefficient, or simply initial, of an in- equality aTi x ≤ bi , for 1 ≤ i ≤ m, the coefficient of aTi x in its largest variable. Following the terminology of W. Pugh in [17], if v is the largest variable of the inequality aTi x ≤ bi , we say that this inequality is an upper (resp. lower) bound of v whenever the initial c of aTi x ≤ bi is positive (resp. negative); indeed, we have v ≤ γc (resp. v ≥ γc ) where γ = bi − aTi x + c v. 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 co- ordinate system. Assume again that K is full-dimensional and that the inter- section ∩i=n i=1 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, ∩i=n i=1 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 0, a > 0, γ ∈ R[x2 , . . . , xd ] and α ∈ R[x2 , . . . , xd ] hold), we have a new inequality cα−aγ ≤ 0. Augmenting A ⋯ > xd and denote by proj(K; x1 , . . . , xd ) the linear system given by Ax1 x ≤ bx1 if d = 1 and, by the conjunction of Ax1 x ≤ bx1 and proj(Π x2 ,...,xd K; x2 , . . . , xd ), otherwise. 3 Integer Solutions of Linear Equation Systems We review how Hermite normal forms [6, 19] can be used to represent the integer solutions of systems of linear equations. Let A = (ai,j ) and H = (hi,j ) be two matrices over Z with m rows and d columns, and let b be a vector over Z with d coefficients. We denote by r the rank of A and by h the maximum bit size of coefficients in the matrix [A b]. Definition 1 is taken from [13], see also [11]. Definition 1. The matrix H is called a column Hermite normal form (abbr. column HNF) if there exists a strictly increasing map f from [d − r + 1, d] ∩ Z to [1, m] ∩ Z satisfying the following properties for all j ∈ [d − r + 1, d] ∩ Z: 1. for all integer i such that i > f (j) holds, we have hi,j = 0, 2. for all integer k such that j < k ≤ d holds, we have hf (j),j > hf (j),k ≥ 0, 3. the first d − r columns of H are equal to zero. 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 Definition 1. Let us decompose U as U = [UL , UR ] where UL (resp. UR ) consist of the first d − r (resp. last r) columns of U . Then we define 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 [17], However, no running-time estimates is known for that procedure. 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 (coefficients of b) with indices in I. Lemma 1. Let I be the pivot row set of A, as given in Definition 1. Assume that Ax = b admits at least one solution in Rd . Then, for any x ∈ Rd , we have Ax = b ⇐⇒ AI x = bI . 8 Lemma 2. We use the same notations as in Definition 1 and Remark 1. We assume that HR is non-singular. Then, the system Ax = b has an integer solution −1 if and only if HR b 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 coefficient in P (resp. q) can be bounded over by rr+1 L2r (resp. rr+1 L2r ), where L is the maximum absolute value of a coefficient 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. Example 1. Let A, H and U be as follows: ⎛ 3 4 −4 −1 ⎞ ⎛ 0 −18 −1 −15 ⎞ ⎜ 2 −2 8 4 ⎟ ⎜ 0 18 2 16 ⎟ ⎛ −1 30 −3 −25 ⎞ ⎜ ⎟ ⎜ ⎟ ⎜ 1 −37 4 31 ⎟ A = ⎜5 2 4 3 ⎟ , H = ⎜0 0 1 1 ⎟ , U =⎜ ⎟ . ⎜ ⎟ ⎜ ⎟ ⎜ 0 −19 2 16 ⎟ ⎜ 3 5 −5 −2 ⎟ ⎜0 0 1 0 ⎟ ⎝ 1 0 0 0 ⎠ ⎝ 2 −3 9 5 ⎠ ⎝0 0 0 1 ⎠ 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 coefficients 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 . (2) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪3x1 + 5x2 − 5x3 − 2x4 = b4 ⎪ ⎩2x1 − 3x2 + 9x3 + 5x4 = b5 ⎪ ⎪ ⎪ ⎪ ⎩2x1 − 3x2 + 9x3 + 5x4 = b5 −1 We apply Lemma 2: if Ax = b is consistent over Q and if HR [b2 , b4 , b5 ]T is inte- gral, 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 = [ 35 b2 − 19 b − 155 3 4 3 5 b , − 18 37 b2 + 73 b + 9 b5 , − 18 b2 + 9 b4 + 9 b5 ] , t = (t1 ) and t1 is a new variable. 9 4 575 19 37 296 T 4 Integer Solutions of Inequality Systems 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 [17] for deciding whether or not a polyhedral set has an integer point. With the notations of Section 2, consider again the polyhedron K: ⎧ ⎪ a x + ⋯ + a1d xd ≤ b1 ⎪ ⎪ 11 1 ⎨ ⋮ ⎪ ⎪ ⎪ ⎩ am1 x1 + ⋯ + amd xd ≤ bm 9 We aim at eliminating x1 , so we consider two inequalities in x1 : – one with ai1 > 0, called upper bound: li ∶ ai1 x1 + ⋯ + aid xd ≤ bi , – another one with aj1 < 0, called lower bound: lj ∶ aj1 x1 + ⋯ + ajd xd ≤ bj . The real projection makes a linear combination rij of li and lj : − aj1 (ai2 x2 + ⋯ + aid xd ) + ai1 (aj2 x2 + ⋯ + ajd xd ) ≤ −aj1 bi + ai1 bj . The dark projection translates rij towards the center of K into dsij : −aj1 (ai2 x2 +⋯+aid xd )+ai1 (aj2 x2 +⋯+ajd xd ) ≤ −aj1 bi +ai1 bj −(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