On Barvinok’s Algorithm for Counting Lattice Points in Fixed Dimension

Martin Dyer^{∗} School of Computer Studies University of Leeds, Leeds, UK

Ravi Kannan Department of Computer Science Yale University, New Haven CT, USA

March 25, 2002

Abstract

We describe a simplification of a recent polynomial-time algorithm of A. I. Barvinok for counting the number of lattice points in a polyhedron in fixed dimension. In particular, we show that only very elementary properties of exponential sums are needed to develop a polynomial-time algorithm.

1 Introduction

In [1], Barvinok gives the first polynomial time algorithm for counting the number of lattice points in a convex polyhedron in any fixed dimension d. This is a significant achievement, improving dramatically on the previously known cases for d ≤ 4 [3]. In fact, Barvinok’s algorithm counts the number of lattice points in a simplex with integer vertices, since it is known [3] that the general problem can be reduced to this case. To do this, the method employs exponential sums. For a given polyhedron P ⊆ R^{d} and vector c ∈ R^{d} this is an expression of the form

∗Supported in part by Esprit Working Group RAND.

σ(P, c) = ∑_{x∈P ∩Z^{d}} e^{−c.x}.

(Note that we have made a sign change from the notation of [1].) In particular, if K is a pointed polyhedral cone generated by the vectors u_{i} (i = 1, …, k), the exponential sum converges provided c.u_{i} > 0 for i = 1, …, k. Barvinok’s solution uses two deeper properties of these sums. The first property is that the sum (regarded as a function of c ∈ C^{d}) can be continued to a define a meromorphic function on C^{d}. The second is an identity of Brion [2] which relates the exponential sum over a polytope to the sums over the cones generated by the edges at each vertex. The analytic continuation is crucial here, since it is impossible to find a single c for which all the required sums converge. This introduces an element of symbolic computation, in order to avoid the complication caused by poles of the sums.

The remainder of Barvinok’s procedure uses an inclusion-exclusion method to replace a sum over an arbitrary cone with a sum over a polynomial (in the size of the data) number of primitive cones. The sum over a primitive cone can be evaluated explicitly. We discuss this in more detail below.

The purpose of this note is to indicate that the non-elementary properties of exponential sums invoked in Barvinok’s algorithm are unnecessary to obtain a polynomial time algorithm for this problem. We give an algorithm which uses only the reduction to primitive cones, replacing the symbolic computation with arithmetic computation. This results from the fact that, in our algorithm, we can compute a c for which all the sums involved converge. We should emphasise, however, that our method still relies heavily on Barvinok’s ideas.

2 Definitions and Notation

The notation we use mainly follows [1], and the reader may seek further information there on some of what follows.

Throughout, K ⊆ R^{d} will be a cone pointed at the origin with linearly independent integral generators u_{i} (i = 1, …, k). Then u_{ij} will denote the jth component of u_{i}. Different cones will be indicated by superscripts. Note that if C = K + v is a cone pointed at an integer point v, then σ(C, c) = e^{−c.v}σ(K, c).

A primitive cone K is a cone having a particular unique minimal set of generators such that x ∈ Z^{d} ∩ K if and only if x can be expressed as a linear combination of these u_{i} with nonnegative integer weights. It follows that, if K is primitive,

σ(K, c) = ∏_{i=1}^{k} \frac{1}{1 − e^{−c.u_{i}}} .

For given cone K, the index ind K of K can be defined in several equivalent ways (see [1]). We will use a computationally useful characterizarion of ind K, as follows. Let U be the d × k matrix [u_{1}u_{2} ⋯ u_{k}] formed by the generators of K. From the Hermite normal form consruction [5], there exists a d × d unimodular matrix T such that

T U = ( \begin{matrix}R \\ 0 \end{matrix} ) ,

where R is a nonsingular upper triangular matrix. Then ind K = | det R|. These computations can be carried out in polynomial time [5]. It follows easily that K is primitive if and only if ind K = 1 and that, if K′ is any face of K, ind K′ ≤ ind K.

Following [1], for integers ϵ_{m} and polyhedra P^{m} (m ∈ M), we will write

P = ∑_{m∈M} ϵ_{m}P^{m}

to mean the corresponding identity on the characteristic functions of the polyhedra, where the addition and scalar multiplication of functions is pointwise. Let us call the right hand side expression a composition of the polyhedra. Note that these do not correspond to the usual operations of addition and scalar multiplication of convex sets. In general, a composition of convex sets would not necessarily be a characteristic function, let alone that of a convex set. However, this will always be the case here.

3 The algorithm

We will first show that an integral simplex S ⊆ R^{d−1} with vertices v_{1}, …, v_{d} can be expressed as a composition of cones pointed at the vertices, such that all cones lie in the interior of some half-space. We will assume without loss that S is contained in the non-negative orthant of R^{d−1}. This can always be arranged with a suitable translation. We wish to determine N, the number of integer points in S. Note that, if γ_{1} = 1 + max_{i,j} v_{ij}, then N < γ_{1}^{d}.

While we could work directly in R^{d−1}, it is easier for exposition to embed the proof in R^{d}. (This is why we choose to start from a simplex in R^{d−1}.) Thus we consider the cone K^{∗} ⊆ R^{d} with generators u_{i}^{∗} = (1, v_{i}) (i = 1, …, d). Then, in an obvious notation, (1, S) = K^{∗} ∩ H_{1}, where H_{1} = {x : x_{1} = 1}.

Let us call K ⊆ R^{d} a standard cone if its generators satisfy u_{11} = 1 and u_{i1} = 0 (i = 2, …, k), and each of the vectors u_{i} is lexicographically positive. Note that, if K is a standard cone, and (1, C) = K ∩ H_{1} then C is a cone with vertex v where (1, v) = u_{1}.

Lemma 1 The cone K^{∗} is a composition of standard cones K^{m} (m ∈ M) with |M| ≤ 3^{d−1}, ϵ_{m} = ±1 and ind K^{m} ≤ ind K (m ∈ M). Moreover, the required composition can be determined in polynomial time.

Proof We will prove the Lemma by induction on r, for a cone K with k generators satisfying u_{i1} = 1 (i = 1, …, r), u_{i1} = 0 (i = r + 1, …, k) where 1 ≤ r ≤ k ≤ d. (We require the case r = k = d.) The basis for the induction is r = 1, when there is nothing to prove.

Assume the truth of the Lemma for (r − 1) and suppose r ≥ 2. We assume without loss that the generators of K are in lexicographically decreasing order. Let w = u_{1} − u_{2}, so w is lexicographically positive and w_{1} = 0. We now use the “inclusion-exclusion” method, as in [3, 1], to “insert” the vector w. Now K′ = K ∪ K′′ ,where K′ has generators w, u_{2}, …, u_{k} and K′′ has generators u_{1}, w, …, u_{k}. Let K′′′ = K ∩ K′′ be the (k − 1)-dimensional cone with generators u_{1}, u_{3}, …, u_{k}. Then clearly

K′ = K + K′′ − K′′′, i.e. K = K′ − K′′ + K′′′.

Thus K is a composition of three cones, each with only (r − 1) generators having nonzero first component. Note that ind K′ = ind K′′ = ind K, since they all generate the same lattice, and ind K′′′ ≤ ind K, since the generators of K′′′ are a subset of those of K. Applying the inductive hypothesis gives the Lemma. The proof clearly indicates a polynomial time algorithm. □

We remark that this type of decomposition is well known. For example, it is the starting point in [6] and it is mentioned there as a “folklore result”. It also leads to an an elementary proof of the Brion’s identity in the spirit of [4].

Now, we follow Barvinok’s method to express each of the cones K^{m} from Lemma 1 as a composition of primitive cones. This is done in the following way. For a given cone K with ind K > 1, determine the matrix R as in section 2 using a Hermite normal form algorithm. Compute the shortest nonzero vector λ, in ℓ_{∞} norm, in the lattice generated by the columns of R^{−1}. This can be done in polynomial time, in fixed dimension, using the basis reduction algorithm [5] followed by enumeration. By Minkowski’s theorem (see e.g. [5, p 71]),

||λ||_{∞} ≤ | det R^{−1}|^{1/k} = | det R|^{−1/k} = (ind K)^{−1/k}.

Now λ = R^{−1}z, for some integral z so z = Rλ and hence

w = T^{−1} ( \begin{matrix}z \\ 0 \end{matrix} ) = T^{−1} ( \begin{matrix}R \\ 0 \end{matrix} ) λ = Uλ,

is a nonzero integer vector which is a linear combination of the generators of K with weights at most (ind K)^{−1/k} in absolute value. We may clearly insist that w is lexicographically positive. We now, using inclusion-exclusion, we can express K as a composition of faces of the cones K^{i} with generators u_{1}, …, u_{i−1}, w, u_{i+1}, …, u_{k}. See [1] for details. There are at most k2^{k} cones of this form. Suppose U_{i} is the matrix of generators of K^{i}. Note that w = Uλ. Let Λ_{i} be the matrix which is a k × k identity, except that the ith column is λ. Then U_{i} = UΛ_{i}. Let T_{i}, T be the unimodular matrices which reduce U_{i}, U respectively to Hermite normal form, then

( \begin{matrix}R_{i} \\ 0 \end{matrix} ) = T_{i}U_{i} = T_{i}UΛ_{i} = T_{i}T^{−1} ( \begin{matrix}RΛ_{i} \\ 0 \end{matrix} ) ,

from which it follows that | det R_{i}| = | det(RΛ_{i})|. Hence

ind K^{i} = λ_{i}ind K ≤ (ind K)^{−1/k}ind K = (ind K)^{(k−1)/k}.

From this it follows that iterating this O(log log ind K) times we obtain cones with ind K = 1, i.e. primitive cones. The total number of cones generated is then only polynomial in log ind K, for fixed k.

At the end of this process, we have expressed the simplex S as a composition of a polynomial number of primitive cones, such that all cones have vertex in the non-negative orthant and all cone generators are lexicographically positive. We now show that there is a polynomial-sized c > 0 such that c.u_{i} > 0 for all generators of all cones.

Lemma 2 Let u_{r} ∈ R^{d} (r = 1, 2, …) be any collection of lexicographically positive integer vectors, such that γ_{2} = 1 + max_{r,j} |u_{r,j}| is bounded. Then the vector c_{0} = (γ_{2}^{d−1}, γ_{2}^{d−2}, …, γ_{2}, 1) satisfies c_{0}.u_{r} ≥ 1 (∀r).

Proof Suppose u_{r,s} ≥ 1 is the first nonzero element of u_{r}. Then

c_{0}.u_{r} ≥ γ_{2}^{d−s} − ∑_{j=0}^{d−s−1}(γ_{2} − 1).γ_{2}^{j} = 1

. □

Thus there exists a c_{0} such that for any c = δc_{0} (0 < δ ≤ 1), σ(S, c) can be expressed as a sum of τ terms of the form

±\frac{e^{−c.v}}{∏_{i=1}^{d}(1 − e^{−c.u_{i}})},  (1)

where c.v > 0 and c.u_{i} > 0 for i = 1, …, d and τ is bounded by a polynomial in the size of the description of S. We now show that we find a suitable δ such that σ(S, c) approximates N closely enough. Let γ = max{γ_{1}, γ_{2}}.

Lemma 3 Let δ = min{1/(4dγ^{2d}), 1/(8τ)} and c = δc_{0}, then

σ(S, c) ≤ N < σ(S, c) + \frac{1}{4}.

Proof For any x ∈ S we have

c.x ≤ δ(dγ_{2}^{d−1})γ_{1} < \frac{1}{4}γ^{−d}

Hence 1 − \frac{1}{4}γ^{−d} < 1 − c.x ≤ e^{−c.x} ≤ 1, and thus

N − \frac{1}{4} < ∑_{x∈S} e^{−c.x} ≤ N,

since N < γ^{d}. □

Clearly the c determined in Lemma 3 is a rational vector of polynomial size in the data.

We now determine an approximation σ̂(S, c) to σ(S, c) such that |σ̂(S, c) − σ(S, c)| ≤ \frac{1}{4}. Then

σ̂(S, c) − \frac{1}{4} ≤ N < σ̂(S, c) + \frac{1}{2},

and thus N is σ̂(S, c) rounded to the nearest integer. Thus it suffices to determine each term (1) to within 1/(4τ).

Let c.v = α, c.u_{i} = β_{i}. Then α ≤ \frac{1}{4}, and since

1 ≤ c_{0}.u_{i} ≤ ∑_{j=0}^{d−1}(γ_{2} − 1)γ_{2}^{j} < γ_{2}^{d},

and δ < \frac{1}{4}, we have

δ ≤ β_{i} ≤ δγ_{2}^{d} < δ(4dδ)^{−1} ≤ \frac{1}{2} \sqrt{δ} < \frac{1}{4}.

Now let a(y) be the approximation to e^{−y} using the first (2d + 3) terms of its series expansion. This series is alternating and we have 0 ≤ α ≤ δ, δ ≤ β_{i} ≤ \frac{1}{2} \sqrt{δ}. It follows that a(α) approximates e^{−α}, and (1 − a(β_{i})) approximates (1 − e^{−β_{i}}) with relative error less than \frac{1}{2}δ^{d+1}/(2d + 3)!. Thus using the approximation a(y) in (1) leads to an approximation with relative error less than (d + 1)δ^{d+1}/(2d + 3)!. Now the term (1) is at most (2/δ)^{d}. Hence the absolute error due to approximating (1) in this manner will be will be at most 2^{d}δ/(2d + 3)! < δ ≤ 1/(8τ). Thus it suffices to replace each exponential by the fixed polynomial a(x) of degree 2(d + 1). It now follows that we can compute the required σ̂(S, c) in polynomial time.

Acknowledgement

We are grateful to Alexander Barvinok for helpful comments.

References

[1] A. I. Barvinok, A polynomial time algorithm for counting integral points in polyhedra when the dimension is fixed, in Proceedings of the 34th Annual Symposium on Foundations of Computer Science, IEEE Computer Society Press, 1993, pp 566–572.

[2] M. Brion, Polyèdres et réseaux, L’enseignement Mathématique 38 (1992), 71–78.

[3] M. E. Dyer, On counting lattice points in polyhedra, SIAM Journal on Computing 20 (1991), 695–707.

[4] A.G. Khovanskii and A.V. Puhlikov, A Riemann-Roch theorem for integrals and sums of quasipolynomials on virtual polytopes (Russian), Algebra i analiz, 4 (1992), N 4, 188–216, translated in St. Petersburg Mathematical Journal, 4 (1992), 789–812.

[5] A. Schrijver, Theory of linear and integer programming, John Wiley, Chichester, 1986.

[6] A.N. Varchenko, Combinatorics and topology of the disposition of affine hyperplanes in real space (Russian), Funktsional. Anal. i Prilozhen., 21 (1987), N 1, 11-22, translated in Funct. Anal. Appl., 21 (1987), 9-19.