ON COUNTING INTEGRAL POINTS IN A CONVEX RATIONAL POLYTOPE

JEAN B. LASSERRE AND EDUARDO S. ZERON

Abstract. Given a convex rational polytope Ω(b) := {x ∈ ℝ_{+}^{n} | Ax = b}, we consider the function b ↦ f(b), which counts the nonnegative integral points of Ω(b). A closed form expression of its ℤ-transform z ↦ F(z) is easily obtained so that f(b) can be computed as the inverse ℤtransform of F. We then provide two variants of an inversion algorithm. As a by-product, one of the algorithm provides the Ehrhart polynomial of a convex integer polytope Ω. We also provide an alternative that avoids the complex integration of F(z) and whose main computational effort is to solve a linear system. This latter approach is particularly attractive for relatively small values of m, where m is the number of nontrivial constraints (or rows of A).

1. Introduction

In this paper, we are interested in the number f(b) of nonnegative integral points x ∈ ℤ^{n}∩Ω where Ω is the convex rational polytope {x ∈ ℝ_{+}^{n} | Ax = b} (that is, the entries of A and b are all in ℤ).

Counting integral points (or, more generally, lattice points) of a convex polytope Ω is an important problem in computational geometry (and operations research as well, in view of its connection with integer programming) which has received much attention in recent years (see e.g. the works of Barvinok [2, 3], Beck [5], Beck, Diaz and Robins [6], Brion [7], Brion and Vergne [8]), Kantor and Khovanskii [12], Khovanskii and Pukhlikov [13]). In particular, using generating functions, Brion and Vergne [8, p. 801] provide generalized residue formulae that yield closed form expressions for f(b) and further exploited in Baldoni-Silva and Vergne [1] for particular cases like flow polytopes. Beck [5] and Beck Diaz and Robins [6] also provide a complete analysis based on residue techniques for the case of a tetrahedron (m = 1) and also mention the possibility of evaluating f(b) for general polytopes by means of residues.

In principle, these theoretical results can be exploited to devise an algorithm to compute f(b) numerically. For instance, Barvinok [2] proposed a conceptual algorithm for rational polytopes with polynomial time computational complexity when the dimension n is fixed. This algorithm requires knowledge of the vertices v of Ω and uses a special representation due to Brion [7] of the generating function g(c) := ∑_{x∈Ω∩ℤ^{n}} c^{x} of the number of integral points in Ω. Essential in the procedure is the triangulation of certain closed convex cones into unimodular cones (see e.g. Barvinok and Pommersheim [3, Theor. 4.4]). Alternatively, Brion and Vergne ’s formula [8, p. 801] could also be exploited in this manner. As mentioned in Baldoni-Silva and Vergne [1] it requires many steps (in particular, one has to build up chambers (maximal cones) in a subdivision of a closed convex cone into polyhedral cones [8]). The approach. The proposed approach is conceptually very simple and does not require explicit knowledge of the vertices of Ω. Following ideas already developed in e.g. Beck [5] and Brion and Vergne [8] and instead of working in the space ℝ^{n} of primal variables x_{i}, we rather work in the space ℂ^{m} of dual variables z_{i} associated with the m nontrivial constraints that define the polytope (by trivial constraint we mean the nonnegativity constraint x ≥ 0). Thus, in a sense, this approach might be seen as “dual” to the ones working in ℝ^{n}. Indeed, one considers the generating function of f(b) instead of the generating function g(c) := ∑_{x∈Ω∩ℤ^{n}} c^{x} as in Barvinok [2].

Date: September 20, 2002.

Research partially supported by the ECOS-Nord (France)-ANUIES (México) Educational and Scientific Cooperation Program PM98M02; and Cinvestav (México).

- We consider the (multi-dimensional) ℤ-transform F(z) of f(b) (as in Brion and Vergne [8]). The dimension of the complex vectors z is m, the number of nontrivial constraints in the half-space description of the polytope.

- We do not exploit the general formulae for f(b) (e.g. as provided in Brion and Vergne [8]) but rather evaluate f(b) numerically by computing the inverse ℤ-transform of F at point b by m (one-dimensional) complex integrations with respect to one variable at a time.

The m-dimensional inversion with respect to (w.r.t.) (z_{1}, ⋯, z_{m}) is done in m steps, where step k consists of several one-dimensional complex integrations w.r.t. the variable z_{k}. These one-dimensional integrations are performed by repeated application of Cauchy ’s Residue Theorem. In fact, we provide two variants of an inversion algorithm and analyze their respective advantages. This approach is similar in spirit to the Laplace transform algorithm in Lasserre and Zeron [15] for computing the volume of a convex rational polytope. In particular, we provide concrete recepies on several possible ways on how to apply Cauchy ’s Residue Theorem as efficiently as possible.

In addition, and as a by product, one of the two algorithms that we propose, provides the Ehrhart polynomial of Ω for an integer polytope Ω. We also briefly describe an approximation procedure with no occurrence of multiple poles all along the (therefore simplified) integration process.

- Finally, we consider an alternative method still based on the ℤ-transform F(z) but which avoids its complex integration. It is purely algebraic and uses Hilbert Nullstellensatz Theorem to provide a special decomposition of the ℤ-transform of f(b) into certain partial fractions. As the integration of each term of the decomposition is easy, the main work is the computation of the coefficients of the polynomials involved in this decomposition, which reduces to solving a linear system as bounds on the degree of the polynomials in this decomposition are available (see e.g. Seidenberg [17] and Kollár [14]). This approach might be a viable alternative, particularly for relatively small values of m.

The paper is organized as follows. In §2 we provide an explicit expression of the ℤ-transform F of f(b). In §3 we describe and analyse an algorithm to invert the so-called associated ℤ-transform of f. For illustration purposes, a simple example is worked out in §4 and a general algorithm is outlined in §5. An approximate algorithm with a simplified integration process is also presented in §6. Finally, an alternative approach is presented in §7.

2. The ℤ transform of f

The notation ℝ_{+} stands for the usual positive closed cone of ℝ. As usual, ℤ denotes the set of relative integers and ℤ_{+} = ℕ = {0, 1, 2, …} the set of natural numbers. We denote by c′ and A′ the respective transpose of the vector c and the matrix A. Finally, given any two vectors z ∈ ℂ^{m} and u ∈ ℤ^{m}, the notation z^{u} and ln(z) stands (respectively) for

(2.1)  z^{u} := z_{1}^{u_{1}}z_{2}^{u_{2}} ⋯ z_{m}^{u_{m}},
(2.2)  ln(z) := [ ln(z_{1}), ln(z_{2}), …, ln(z_{m}) ].

As mentioned in the introduction, we consider the convex polytope

(2.3)  Ω(y) = {x ∈ ℝ^{n} | Ax = y; x ≥ 0},

where y ∈ ℤ^{m} and A ∈ ℤ^{n×m}, and we want to compute the number of points x ∈ ℕ^{n} of Ω(y), that is, the cardinality of the set

(2.4)  ℕ^{n} ∩ Ω(y).

We will actually calculate the following related function

(2.5)  y ↦ f(y) := ∑_{ℕ^{n}∩Ω(y)}e^{c′x}.

for a given vector c ∈ ℝ^{n}.

We trivially have that f(y) is equal to the cardinality of ℕ^{n} ∩ Ω(y) when c = 0 (and for more details on f(y) (with c = 0) the reader is referred to Beck [4]). Morever, as observed in Barvinok and Pommersheim [3], taking c very small and rounding f(y) to the nearest integer (or taking an appropriate residue) will give the number of points x ∈ ℕ^{n} of Ω(y).

Of course, computing the number of points x ∈ ℕ^{n} of the convex polytope

(2.6)  Ω_{1}(y) = {x ∈ ℝ_{+}^{n} | A_{1}x ≤ y}

for some A_{1} ∈ ℤ^{m×n}, reduces to computing the cardinality of ℕ^{n} ∩Ω(y) with Ω as in (2.3) and with A := [A_{1} | I] (I ∈ ℕ^{m×m} being the identity matrix).

Remark 2.1. The set Ω(y) in (2.3) is a convex polytope, that is, a compact convex polyhedron, if and only if its recession cone {x ∈ ℝ_{+}^{n} | Ax = 0} is the singleton {0} (see e.g. Hiriart-Urruty and Lemarechal [11, p. 109]). By a specialized version of Farkas Lemma due to Carver, the latter is in turn is equivalent to stating that the cone

(2.7)  𝕂 := {u ∈ ℝ^{m} | A′u > 0}

is not empty (see e.g. Schrijver [16, 7, (33), p. 95]).

We want to compute the exact value of f(y) and as we will see, an appropriate tool to analyse such a function is the ℤ-transform. Thus, the m-dimensional two-sided ℤ-transform, F : ℂ^{m}→ ℂ, of f is defined by the Laurent series

(2.8)  z ↦ F(z) := ∑_{y∈ℤ^{m}} f(y)z^{−y},

with z^{−y} as in (2.1).

Now let e_{m} = (1, 1, …) be the unit vector in ℝ^{m}. We have the following result:

Theorem 2.2. Let f and F be as in (2.5) and (2.8) respectively, with c ≤ A′u_{o} for some u_{o} ∈ 𝕂. Then :

(2.9)  F(z) = ∏_{k=1}^{n} \frac{1}{(1 − e^{c_{k}}z_{1}^{−A_{1k}}z_{2}^{−A_{2k}} ⋯ z_{m}^{−A_{mk}})},

on the domain

(2.10)  ( |z_{1}|, |z_{2}|, …, |z_{m}|) ∈ {β ∈ ℝ_{+}^{n} | A′ ln(β) > c}.

Moreover, for every y ∈ ℤ^{n},

(2.11)  f(y) = \frac{1}{(2πi)^{m}} ∫_{|z_{1}|=w_{1}} ⋯ ∫_{|z_{m}|=w_{m}} F(z)z^{y−e_{m}}dz,

where the real constant vector w ∈ ℝ_{+}^{m} satisfies A′ ln(w) > c.

Proof. Apply the definition (2.8) of F to obtain :

F(z) = ∑_{y∈ℤ^{m}} z^{−y} [ ∑_{Ax=y,x∈ℕ^{n}} e^{c′x}] = ∑_{x∈ℕ^{n}} e^{c′x} z^{−(Ax)}.

On the other hand, notice that

e^{c′x}z_{1}^{−(Ax)_{1}}z_{2}^{−(Ax)_{2}} ⋯ z_{m}^{−(Ax)_{m}} = ∏_{k=1}^{m} (e^{c_{k}}z_{1}^{−A_{1k}}z_{2}^{−A_{2k}} ⋯ z_{m}^{−A_{mk}})^{x_{k}}

Hence, the conditions |z_{1}^{A_{1k}}z_{2}^{A_{2k}} … z_{m}^{A_{mk}}| > e^{c_{k}} for k = 1, 2, ⋯ n, or equivalently, A′(ln |z_{1}|, ln |z_{2}|, …, ln |z_{m}|) > c, yields

F(z) = ∏_{k=1}^{n} ∑_{x_{k}=0}^{∞} (e^{c_{k}}z_{1}^{−A_{1k}}z_{2}^{−A_{2k}} ⋯ z_{m}^{−A_{mk}})^{x_{k}}

= ∏_{k=1}^{n} \frac{1}{(1 − e^{c_{k}}z_{1}^{−A_{1k}}z_{2}^{−A_{2k}} ⋯ z_{m}^{−A_{mk}})}

which is (2.9). Finally, equation (2.11) is obtained by analyzing the integral ∫_{|z|=r} z^{w}dz with r > 0. This integral is equal to 2πi only if w = −1, whereas if w is any integer different of −1, then the integral is equal to zero. It remains to show that indeed, the domain {β ∈ ℝ_{+}^{n} | A′ ln(β) > c} is not empty. But this follows directly from our choice of the vector c and from Remark 2.1 (take β := e^{2u_{0}}). □

3. Inversion of the ℤ-transform F

Theorem 2.2 allows us to compute f(y) for y ∈ ℤ^{m} via (2.11), that is, by computing the inverse ℤ-transform of F(z) at the point y. Moreover, we can directly calculate (2.11) by using Cauchy’s Residue Theorem because F(z) is a rational function with only a finite number of poles (with respect to one variable at a time). We will call this technique the direct ℤ inverse. We will get back to this in §5.2.

On the other hand, we can also simplify the inverse problem and invert what we call the associated ℤ-transform which yields some advantages when compared to the direct inversion (see the discussion in §5.2).

3.1. The associated ℤ-transform. Assume with no loss of generality that y ∈ ℤ^{m} is such that y_{1} ≠ 0. We may also suppose (without lost of generality) that each y_{i} is a multiple of y_{1} (taking 0 to be multiple of any other integer). Otherwise, we just need to multiply each constraint (Ax)_{i} = y_{i} by y_{1} ≠ 0 when i = 2, 3, …, m, so that the new matrix A and vector y still have entries in ℤ.

Hence, there exists a vector D ∈ ℤ^{m} with first entry D_{1} = 1 and such that y = Dy_{1}. Notice that D may have entries equal to zero or even negative, but not the first one. The inversion problem is thus reduced to evaluate, at the point t := y_{1}, the function g : ℤ→ ℕ defined by

(3.1)  g(t) := f(Dt) = \frac{1}{(2πi)^{m}} ∫_{|z_{m}|=w_{m}} ⋯ ∫_{|z_{1}|=w_{1}} F(z)z^{Dt−e_{m}}dz,

where F is given in (2.9), e_{m} := (1, 1, ⋯) is the unit vector in ℝ^{m} and the real (fixed) vector w ∈ ℝ_{+}^{m} satisfies A′ ln(w) > c. The following technique permits us to calculate (3.1).

Consider the following simple change of variables. Let p = z^{D} and d = w^{D} in (3.1), so that z_{1} = p ∏_{j=2}^{m} z_{j}^{−D_{j}} and

g(t) = \frac{1}{(2πi)^{m}} ∫_{|z_{m}|=w_{m}} … ∫_{|z_{2}|=w_{2}} [∫_{|p|=d} p^{t−1}F̂ dp] dz_{2} … dz_{m},

where

(3.2)  F̂ = F (p ∏_{j=2}^{m} z_{j}^{−D_{j}}, z_{2}, …, z_{m}) ∏_{j=2}^{m} z_{j}^{−D_{j}}.

We can rewrite g(t) as follows

(3.3)  g(t) = \frac{1}{2πi} ∫_{|p|=d} p^{t−1}G(p)dp,  with
(3.4)  G(p) := \frac{1}{(2πi)^{m−1}} ∫_{|z_{2}|=w_{2}} … ∫_{|z_{m}|=w_{m}} F̂ dz_{2} … dz_{m},

and G is called the associated ℤ-transform of f with respect to D.

Recall that F(z) is well defined on the domain given by conditions (2.10), so the domain of definition of F̂ is given by the conditions

(3.5)  ( |p| ∏_{j=2}^{m} |z_{j}|^{−D_{j}}, |z_{2}|, …, |z_{m}|) ∈ {β ∈ ℝ_{+}^{n} | A′ ln(β) > c}.

3.1.1. The Ehrhart polynomial. Consider the convex polytope Ω_{1}(y) defined in (2.6) with y = D and the dilated polytope tΩ_{1}(D) := {tx | x ∈ Ω_{1}(D)}. Then with A := [A_{1} | I] and c := 0, consider the function t ↦ g(t) defined in (3.1).

When Ω_{1}(D) is an integer polytope (an integer polytope has all its vertices in ℤ^{m}) then g is a polynomial in t, that is,

(3.6)  g(t) = a_{n}t^{n} + a_{n−1}t^{n−1} + ⋯ + a_{0},

with a_{0} = 1 and a_{n} = volume(Ω_{1}(D)). It is called the Ehrhart polynomial of Ω_{1}(D) (see e.g. Barvinok and Pommersheim [3], Ehrhart [10]).

On the other hand, if Ω_{1}(D) is not an integer polytope, then the function t ↦ g(t) is a quasipolynomial, that is, a polynomial in the variable t like in (3.6) but now with coefficients {a_{i}(t)} that are periodic functions of t, for instance as in the illustrative example (4.5) below.

3.1.2. Inversion of the associated ℤ-transform. To compute g(t) we proceed as follows.

- We first compute G(p) in (3.4), in m − 1 steps, where each step k is a one-dimensional integration w.r.t. z_{k}, k = 2, …, m.

- We then compute g(t) in (3.3), a one-dimensional integration w.r.t. p. At each of these steps, the one-dimensional complex integrals are eval-

uated by Cauchy ’s Residue Theorem that we now describe in the present context.

3.2. Calculating integrals by residues. One of the easiest ways of calculating (3.3) and (3.4) is to use Cauchy’s Residue Theorem. In our context, we have to use this theorem at each of the m one-dimensional integration steps. We are going to see (cf. example in §4 below) that we have to integrate several times (at each step and along a circle |z| = w) a rational function of the following kind:

(3.7)  R(z) = \frac{α_{0} + α_{1}z + … + α_{o_{1}}z^{o_{1}}}{∏_{k=1}^{o_{2}}(z^{d_{k}} − β_{k})^{δ_{k}}} ,

where each d_{k} and δ_{k} are positive integers. This rational function can obviously be re-written as follows:

(3.8)  R(z) = \frac{α_{0} + α_{1}z + … + α_{o_{1}}z^{o_{1}}}{∏_{k=1}^{o_{3}}(z − β_{k}^{∗})^{η_{k}}} ,

where each η_{k} is a positive integer and the coefficients β_{k}^{∗} are pairwise distinct. We can integrate (3.8) by using Cauchy ’s Residue Theorem, which can be done in several ways, for instance, by the three different techniques (a), (b) and (c) proposed below.

(a) One way to proceed is as follows. Suppose that w > 0 and |β_{k}^{∗}| ≠ w for k = 1, 2, …, o_{3}. Then

(3.9)  \frac{1}{2πi} ∫_{|z|=w} R(z) dz = ∑_{|β_{k}^{∗}|<w} Res (R, β_{k}^{∗}).

However, calculating residues in (3.9) is not always a simple task, mainly when some η_{k} is large. Therefore, we next propose an alternative technique.

(b) Consider a real number w^{∗} > 0 big enough to ensure |β_{k}^{∗}| < w^{∗} for every k = 1, 2, …, o_{3}. Then

∑_{k=1}^{o_{3}} Res (R, β_{k}^{∗}) = \frac{1}{2πi} ∫_{|z|=w^{∗}} R(z) dz = \frac{1}{2πi} ∫_{|v|=1/w^{∗}} \frac{R(1/v) dv}{v^{2}} .

Notice that the change of variable z = 1/v gives us a negative sign in last integral, but this sign gets canceled because we also change the orientation of the integration path. Moreover, the function

(3.10)  \frac{R(1/v)}{v^{2}} = \frac{α_{0}v^{o_{1}} + α_{1}v^{o_{1}−1} + … + α_{o_{1}}}{v^{o_{4}} ∏_{k=1}^{o_{3}}(1 − vβ_{k})^{η_{k}}} ,
(3.11)  o_{4} = 2 + o_{1} − (η_{1} + η_{2} + … + η_{o_{3}}),

is analytic inside the circle |v| = 1/w^{∗} when o_{4} ≤ 0 because |β_{k}^{∗}| < w^{∗} for every k = 1, 2, …, o_{3}. Hence:

∑_{k=1}^{o_{3}} Res (R, β_{k}^{∗}) = \{ \begin{matrix}0 & if & o_{4} ≤ 0, \\ α_{o_{1}} & if & o_{4} = 1. \end{matrix}

We may not wish calculate the above sum when o_{4} ≥ 2, because it can become too complicated; namely, we have to calculate the (o_{4} −1) derivative of R(1/v)v^{o_{4}−2}, and it is obviously prohibitive for large o_{4}. On the other hand, combining the above sum with (3.9) allows us to calculate the integral of R along the circle |z| = w by considering residues of poles outside the integration path. Recall that o_{4} ≤ 1 is given by (3.11). Hence,

(3.12) \frac{1}{2πi} ∫_{|z|=w} R(z) dz = − ∑_{|β_{k}^{∗}|>w} Res (R, β_{k}^{∗}) + \{ \begin{matrix}0 & if & o_{4} ≤ 0, \\ α_{o_{1}} & if & o_{4} = 1. \end{matrix}

(c) Alternatively, one may also calculate directly the integral of (3.7) by using (3.9) and (3.12). Suppose that the real numbers |β_{k}|^{1/d_{k}}, for k = 1, ...o_{2}, are all different by pairs, then the function R has the following expansion into simple fractions,

(3.13)  R(z) = Q_{0}(z) + ∑_{k=1}^{o_{2}} \frac{Q_{k}(z)}{(z^{d_{k}} − β_{k})^{δ_{k}}};

where Q_{0} and each Q_{j} are polynomials; degree of Q_{k} is strictly less than d_{k}δ_{k} for every k = 1, 2, …, o_{2}; and exponents d_{k} and δ_{k} are given by (3.7). Previous expansion of R can be done in time polynomial with respect to its numerator degree, and in time exponential with respect to its denominator degree. On the other hand, we now integrate each term of R along a circle |z| = w. Suppose that |β_{k}| ≠ w^{δ_{k}} for every k = 1, 2, …, o_{2}. Notice that the integral of Q_{0}(z) is equal to zero because it is a polynomial. Moreover, from (3.9) and (3.12), we have that:

∫_{|z|=w} \frac{Q_{k}(z)}{(z^{δ_{k}} − β_{k})^{d_{k}}} dz = 0  if |β_{k}| > w^{δ_{k}} or Degree(Q_{k}) ≤ δ_{k}d_{k} − 2.

Therefore,

(3.14)  \frac{1}{2πi} ∫_{|z|=w} R(z) dz = ∑_{|β_{k}|<w^{δ_{k}}}(δ_{k}d_{k} − 1)-coefficient of Q_{k},

where each Q_{k} is given by (3.13) and we suppose that the (δ_{k}d_{k} − 1)- coefficient of Q_{k} is equal to zero if deg(Q_{k}) ≤ δ_{k}d_{k} − 2.

Before we outline the algorithm in the general case we briefly illustrate the method on a simple example.

4. Illustrative example

Consider the following convex rational polytope with three (m = 3) nontrivial constraints.

Ω_{1}(te_{3}) := {x ∈ ℝ_{+}^{2} | x_{1} + x_{2} ≤ t; −2x_{1} + 2x_{2} ≤ t and 2x_{1} − x_{2} ≤ t}.

With

c = 0, A := [ \begin{matrix}1 & 1 & 1 & 0 & 0 \\ −2 & 2 & 0 & 1 & 0 \\ 2 & −1 & 0 & 0 & 1 \end{matrix} ] ,

and by Theorem 2.2, we have to calculate the inverse ℤ-transform of :

(4.1)  F = \frac{z_{1}z_{2}z_{3}}{(z_{1} − 1)(z_{2} − 1)(z_{3} − 1)(1 − z_{1}^{−1}z_{2}^{2}z_{3}^{−2})(1 − z_{1}^{−1}z_{2}^{−2}z_{3})},

where

(4.2)  \{ \begin{matrix}|z_{j}| > 1 for j = 1, 2, 3, \\ |z_{1}z_{2}^{−2}z_{3}^{2}| > 1, \\ |z_{1}z_{2}^{2}z_{3}^{−1}| > 1. \end{matrix}

We wish to work with rational functions whose denominator’s degree is the smallest possible, so we are going to fix z_{1} = p/(z_{2}z_{3}) and divide by z_{2}z_{3} (see 3.2) because z_{1} has the exponents with smallest absolute value.

F̂ = \frac{z_{3}p^{2}}{(z_{2}^{−1}p − z_{3})(z_{2} − 1)(z_{3} − 1)(z_{3} − z_{2}^{3}p^{−1})(z_{2}p − z_{3}^{2})},

where

\{ \begin{matrix}|p| > |z_{2}z_{3}|, \\ |z_{3}p| > |z_{2}^{3}| > 1, \\ |z_{2}p| > |z_{3}^{2}| > 1. \end{matrix}

Notice that z_{2}^{∗} = z_{3}^{∗} = 2 and p^{∗} = 5 is a solution of the previous system of inequalities. We are going to calculate (3.3) and (3.4) by fixing w_{2} = z_{2}^{∗}, w_{3} = z_{3}^{∗} and d = p^{∗}. Let us integrate F̂ along the circle |z_{3}| = z_{3}^{∗} with a positive orientation. It is easy to see that (taking p := p^{∗} and z_{2} := z_{2}^{∗} constant) F̂ has two poles located on the circle of radius |z_{2}p|^{1/2} > z_{3}^{∗}, and three poles located on the circles of radii 1 < z_{3}^{∗}, |z_{2}^{−1}p| > z_{3}^{∗} and |z_{2}^{3}p^{−1}| < z_{3}^{∗}. Whence, we can consider poles inside the circle |z_{3}| = z_{3}^{∗}, in order to avoid analyzing the pole z_{3} := (z_{2}p)^{1/2} with fractional exponent. This yields

(4.3)  I_{1}(p, z_{2}) = \frac{z_{2}p^{2}}{(p − z_{2})(z_{2} − 1)(p − z_{2}^{3})(z_{2} − p^{−1})}

+ \frac{z_{2}^{3}p^{5}}{(p^{2} − z_{2}^{4})(z_{2} − 1)(z_{2}^{3} − p)(p^{3} − z_{2}^{5})}. (4.4)

Next, we integrate I_{1} along the circle |z_{2}| = z_{2}^{∗}. Taking p := p^{∗} as a constant, the first term of I_{1} has poles on circles of radii |p| > z_{2}^{∗}, 1 < z_{2}^{∗}, |p|^{1/3} < z_{2}^{∗} and |p|^{−1} < z_{2}^{∗}. We consider poles outside the circle |z_{2}| = z_{2}^{∗}, in order to avoid analysing the pole z_{2} := p^{1/3} (recall equation (3.12)). We obtain

\frac{− p^{3}}{(p − 1)(p^{2} − 1)^{2}} .

The second term of I_{1} has poles on circles of radii |p|^{1/2} > z_{2}^{∗}, 1 < z_{2}^{∗}, |p|^{1/3} < z_{2}^{∗} and |p|^{3/5} > z_{2}^{∗}. Notice that we have poles with fractional exponents inside and outside the integration path |z_{2}| = z_{2}^{∗}, so we use equation (3.14). Expanding the second term of I_{1} into simple fractions, yields

\frac{p^{5}}{(z_{2} − 1)(p^{2} − 1)(1 − p)(p^{3} − 1)}

+ \frac{(2p^{3} + 2p^{4} + 3p^{5} + p^{6} + p^{7})z_{2}^{2} + α_{1}(p)z_{2} + α_{0}(p)}{(z_{2}^{3} − p)(p − 1)^{3}(p + 1)^{2}(p^{2} + 1)}

+ \frac{Q_{1}(z, p)}{(z_{2}^{4} − p^{2})} + \frac{Q_{2}(z, p)}{(z_{2}^{5} − p^{3})},

obviously, we only calculate those terms with denominator (z_{2}−1) or (z_{2}^{3}−p). As a result, the associated ℤ-transform G(p) is given by

G(p) = \frac{− p^{3}}{(p − 1)^{3}(p + 1)^{2}} − \frac{p^{5}}{(p − 1)^{2}(p + 1)(p^{3} − 1)}

+ \frac{2p^{3} + 2p^{4} + 3p^{5} + p^{6} + p^{7}}{(p − 1)^{3}(p + 1)^{2}(p^{2} + 1)} .

Finally, it remains to integrate G(p)p^{t−1} along the circle |p| = p^{∗}. Observe that G(p)p^{t−1} has poles when p is equal to 1, −1, i = \sqrt{−1}, −i, σ = e^{2πi/3} and σ̄. Hence, for t ∈ ℤ_{+}, we finally obtain

(4.5)  g(t) = \frac{17t^{2}}{48} + \frac{41t}{48} + \frac{139}{288} + \frac{t(−1)^{t}}{16} + \frac{9(−1)^{t}}{32}

+\frac{\sqrt{2}}{8} cos (\frac{πt}{2} + \frac{π}{4} ) + \frac{2}{9} cos (\frac{2πt}{3} + \frac{π}{3} ) .

Notice that we are using the following two identities which can be easily calculated.

\sqrt{2} cos (\frac{πt}{2} + \frac{π}{4} ) = \frac{(i − 1)i^{t}}{2i} + \frac{(i + 1)(−i)^{t}}{2i} ,
2 cos (\frac{2πt}{3} + \frac{π}{3} ) = \frac{(σ − 1)σ^{t}}{σ − σ̄} + \frac{(σ̄ − 1)σ̄^{t}}{σ̄ − σ} .

It is now easy to count the number f(te_{3}) = g(t) of points x ∈ ℕ^{2} of Ω(te_{3}).

f(0) = \frac{139}{288} + \frac{9}{32} + \frac{1}{8} + \frac{1}{9} = 1
f(e_{3}) = \frac{17}{48} + \frac{41}{48} + \frac{139}{288} − \frac{1}{16} − \frac{9}{32} − \frac{1}{8} − \frac{2}{9} = 1
f(2e_{3}) = \frac{68}{48} + \frac{82}{48} + \frac{139}{288} + \frac{2}{16} + \frac{9}{32} − \frac{1}{8} + \frac{1}{9} = 4

The Ehrhart polynomial.

Observe that in (4.5) the function t ↦ g(t) is a quasipolynomial because the linear and constant term have periodic coefficients. However, with t = 12, one may check that the rational polytope Ω_{1}(te_{3}) is in fact an integer polytope, that is, all its vertices are in ℤ^{n}.

Therefore, Ω_{1}(12te_{3}) is the t-dilated polytope of Ω_{1}(12e_{3}) and g̃(t) := g(12t) with g(t) as in (4.5), is thus the Ehrhart polynomial of Ω_{1}(12e_{3}) (see Ehrhart [10], Barvinok and Pommersheim [3]). We obtain

(4.6)  g̃(t) = 51t^{2} + 11t + 1,

and indeed, 51 is the volume of Ω_{1}(12e_{3}) and the constant term is 1 as it should be (see the discussion in §4.6).

5. A general algorithm

We here describe a general algorithm, and in fact, two variants. The first one is via the inversion of the associated ℤ-transform whereas the second one is the direct inversion of the ℤ-transform.

5.1. The inverse associated ℤ-transform algorithm.

5.1.1. Sketch of the algorithm. We want to calculate the inverse ℤ-transform f of an analytic function F which is well defined on the open domain

E_{1} = {(z_{k} := β_{k}e^{iθ_{k}}) ∈ ℂ^{m} | β, θ ∈ ℝ^{m}, A′ ln(β) > c};

see (2.9), (2.10) and (2.11). Moreover, we simplify our problem of calculating f(y) with y ∈ ℤ^{m} by supposing, let us say, that y_{1} ≠ 0 divides every other entry y_{j}. That is, y = Dy_{1} where the first entry of the vector D ∈ ℤ^{m} is equal to one. Then, we calculate the associated ℤ-transform F̂ by doing the change of variable p = z^{D} and dividing by z_{2}z_{3} ⋯ z_{m}, as in (3.2). From this change of variable we can deduce that

F̂(z_{2}, z_{3}, …, z_{m}, p) is well defined on a domain E_{2} ⊂ ℂ^{m}.

The algorithm consists of

- m − 1 intermediate steps to compute G(p) in (3.4) by successive onedimensional integrations w.r.t. z_{2}, z_{3}, …, z_{m}, respectively, using either equation (3.9), or (3.12), or (3.14) (cf. (a), (b) and (c) in §3.2).

- a final step to integrate G(p) in (3.3) (integration w.r.t. p).

We obviously fix a real vector (w_{2}, w_{3}, …, w_{m}, d) ∈ E_{2} ∩ ℝ^{m} in order to compute previous integrals; it is easy to see that E_{2} ∩ ℝ^{m} ≠ ∅ from the definition (3.5) of E_{2}.

Of course, at each step, we want to use among the equations (3.9), (3.12) and (3.14), the less time consuming one. That is, we first try (3.9) or (3.12) and rather choose to use (3.14) if we have to deal with fractional exponents (recall that Cauchy ’s Residue Theorem cannot be applied to non-analytic functions).

5.1.2. Removable singularities. Let us analyze integral (3.4) because it needs some care to handle removable singularities. Suppose we have calculated the following partial integration, 2 ≤ s ≤ m,

I_{s}^{∗} = \frac{1}{(2πi)^{s−1}} ∫_{|z_{2}|=w_{2}} … ∫_{|z_{s}|=w_{s}} F̂ dz_{2} … dz_{s}.

I_{s}^{∗} is analytic on E_{2} ∩ ({(w_{2}, …, w_{s})} × ℂ^{m−s+1}), because F̂ is analytic on the open domain E_{2}, and each integration path |z_{j}| = w_{j} is a compact circle. However, some care is needed. For instance, had we set w_{2} = w_{3} = 2 and d = 8 (instead of w_{2} = w_{3} = 2 and d = 5) in the above example, both terms (4.3) and (4.4) of I_{1} have three poles given by z_{2}^{3} := p, exactly on the path of integration (recall that |z_{2}| = |p^{1/3}| = 2). Those poles are called removable singularities because I_{1} is actually analytic there (we have just shown that), and we only need to add together both terms of I_{1}.

Fortunately, removable singularities is a pathological case which happens with probability zero in a set of problems with randomly generated data A ∈ ℤ^{m×n} and this issue could be ignored in practice.

However, for the sake of mathematical rigor, we must provide a solution. This has actually be done already. Indeed, if any term of I_{s}^{∗} has poles on the integration path |z_{s+1}| = w_{s+1} or |p| = d, we only need to add together all terms which have poles there. For example, adding both terms of I_{1} together, we get

I_{1}(p, z_{2}) = \frac{p^{2}z_{2}(p^{4} + p^{2}z_{2}^{2} − pz_{2}^{3} − z_{2}^{6})}{(p − z_{2})(z_{2} − 1)(z_{2} − p^{−1})(p^{2} − z_{2}^{4})(p^{3} − z_{2}^{5})}.

Whence, taking p := d as a constant, we have that I_{1} has actually poles only on circles of radii |p| > w_{2}, 1 < w_{2}, |p|^{−1} < w_{2}, |p|^{1/2} > w_{2} and |p|^{3/5} > w_{2}. We consider poles inside the circle |z_{2}| = z_{2}^{∗}, getting the associated ℤ-transform

G(p) = \frac{p^{3}(1 + 2p + p^{2} + p^{3})}{(p − 1)^{2}(p + 1)(p^{3} − 1)} − \frac{p^{6}}{(p − 1)^{3}(p + 1)^{2}(p^{2} + 1)},

which can be easily checked to be equal to the one obtained some paragraphs above, we only need to move terms around.

5.2. Direct ℤ inverse algorithm. We calculated the inverse ℤ-transform of (4.1) via the associated ℤ-transform, because it seems to be the easier way. However, we could have instead computed the integral (2.11) directly. That is, we could have computed the inverse ℤ-transform (4.1) at the point y.

5.2.1. Sketch of the algorithm. In the general case, we want to compute the inverse ℤ-transform f of an analytic function F which is well defined on an open domain E_{1} ⊂ ℂ^{m}; see (2.9), (2.10) and (2.11).

As for the associated ℤ inverse algorithm of §5.1, the algorithm consists of m successive steps to compute (2.11), where step k is a one-dimensional integration w.r.t. z_{k}, for all k = 1, …, m.

In order to integrate F in (2.11), we first define the integration path by fixing a real vector (w_{1}, w_{2}, …, w_{m}) ∈ E_{1} ∩ ℝ^{m}. From the definition (2.10) of E_{1} it follows that E_{1} ∩ ℝ^{m} ≠ ∅.

Of course, at each step, we want to use among the equations (3.9), (3.12) and (3.14), the less time consuming one. However, in contrast to the associated transform algorithm, in the present context, equations (3.12) and (3.14) may become really time prohibitive, so we are almost restricted to use only equation (3.9).

Moreover, the function I_{s}^{∗} obtained in the s-th steep (1 ≤ s < m) is analytic on E_{1} ∩ ({(w_{1}, …, w_{s})} × ℂ^{m−s}), because F is analytic on E_{1}, and each integration path |z_{j}| = w_{j} is a compact circle (cf. §5.1.2). Hence, we add together all terms of I_{s}^{∗} which have removable singularities, where in the present case, removable singularities means either poles on the path of integration |z_{s+1}| = w_{s+1} or variables with fractional exponent.

We illustrate below this algorithm by analyzing again the convex polytope Ω(te_{3}) presented in §4. Example 5.1. Consider the polytope Ω_{1}(te_{3}) already treated in §4 (with c = 0). We have to calculate (2.11) when F is given by (4.1). Notice that z_{1}^{∗} = 3 and z_{2}^{∗} = z_{3}^{∗} = 2 is a solution of (4.2), so we are going to integrate (2.11) by fixing w_{1} = 3 and w_{2} = w_{3} = 2. Let us integrate

F(z)z^{(t−1)e_{3}} = \frac{z_{1}^{t+1}z_{2}^{t+2}z_{3}^{t+2}}{(z_{1} − 1)(z_{2} − 1)(z_{3} − 1)(z_{1}z_{3}^{2} − z_{2}^{2})(z_{2}^{2} − z_{1}^{−1}z_{3})}

first along |z_{2}| = w_{2}. Supposing z_{1} and z_{3} constant, we have poles located on circles of radii 1 < w_{2}, |z_{1}^{1/2}z_{3}| > w_{3} and |z_{1}^{−1}z_{3}|^{1/2} < w_{3}. Now, we do not want to consider poles outside the integration path, as in (3.12), because the work required increases with t; see the paragraph just before equation (3.12). Moreover, we do not want to decompose F(z)z^{(t−1)e_{3}} into a sum of simple fractions either, as in (3.13), because this expansion is done in time polynomial in the parameter t as well. Therefore, we may wish to consider only inside poles (the alternative (c) in §3.2), which yields

I_{1}(z_{1}, z_{3}) = \frac{z_{1}^{t+1}z_{3}^{t+2}}{(z_{1} − 1)(z_{3} − 1)(z_{1}z_{3}^{2} − 1)(1 − z_{1}^{−1}z_{3})}

+ \frac{z_{1}^{(t+1)/2}z_{2}^{(3t+5)/2}}{2 (z_{1} − 1)(z_{3} − 1)(z_{1}^{−1/2}z_{3}^{1/2} − 1)(z_{1}z_{3}^{2} − z_{1}^{−1}z_{3})}

− \frac{(−1)^{t}z_{1}^{(t+1)/2}z_{2}^{(3t+5)/2}}{2 (z_{1} − 1)(z_{2} − 1)(−z_{1}^{−1/2}z_{3}^{1/2} − 1)(z_{1}z_{3}^{2} − z_{1}^{−1}z_{3})}.

We cannot work with fractional exponents, for neither z_{1}^{1/2} nor z_{3}^{1/2} are analytic at zero, so we adopt the same remedy suggested in §5.1.2, i.e., we add together the last two terms, and we obtain

I_{1}(z_{1}, z_{3}) = \frac{z_{1}^{t+2}z_{3}^{t+2}}{(z_{1} − 1)(z_{3} − 1)(z_{1}z_{3}^{2} − 1)(z_{1} − z_{3})}

+ \frac{z_{1}^{(t+5)/2} z_{3}^{(3t+3)/2} Q(z_{1}, z_{3})}{2 (z_{1} − 1)(z_{3} − 1)(z_{3} − z_{1})(z_{1}^{2}z_{3} − 1)},

where

Q(z_{1}, z_{3}) = z_{1}^{−1/2}z_{3}^{1/2} + 1 + (−1)^{t}(z_{1}^{−1/2}z_{3}^{1/2} − 1).

Notice that I_{1} has no more fractional exponents, for its second term’s numerator is equal to

\{ \begin{matrix}2 z_{1}^{t/2+2}z_{3}^{3t/2+2} & if & t is even, \\ 2 z_{1}^{(t+5)/2} z_{3}^{(3t+3)/2} & if & t is odd. \end{matrix}

To simplify the exposition, we fix the value of t, say t = 120 for instance. Then

I_{1}(z_{1}, z_{3}) = \frac{z_{1}^{122}z_{3}^{122}}{(z_{1} − 1)(z_{3} − 1)(z_{1}z_{3}^{2} − 1)(z_{1} − z_{3})}

+ \frac{z_{1}^{62}z_{3}^{182}}{(z_{1} − 1)(z_{3} − 1)(z_{3} − z_{1})(z_{1}^{2}z_{3} − 1)}.

We integrate the first term of I_{1} along |z_{1}| = w_{1}, so we suppose z_{3} := 2 constant, and we have poles on circles of radii 1 < w_{1}, |z_{3}^{−2}| < w_{1} and |z_{3}| < w_{1}. We get

I(z_{3}) := \frac{− z_{3}^{122}}{(z_{3} − 1)^{3}(z_{3} + 1)} + \frac{z_{3}^{−120}}{(z_{3} − 1)^{2}(z_{3} + 1)(z_{3}^{3} − 1)} + \frac{z_{3}^{244}}{(z_{3} − 1)^{2}(z_{3}^{3} − 1)}.

Finally, we integrate I(z_{3}) along |z_{3}| = w_{3}. The integral of its second term is equal to zero because of (3.12), whereas the first and third integral yield

− \frac{29281}{8} + \frac{1}{8} + \frac{88208}{9} + \frac{1}{9} = 6141.

Next, we integrate the second term of I_{1}(z_{1}, z_{3}) along the circle |z_{3}| = w_{3}, so we have poles on circles or radii 1 < w_{3}, |z_{1}| > w_{3} and |z_{1}^{−2}| < w_{3}. This yields

I(z_{1}) := \frac{− z_{1}^{62}}{(z_{1} − 1)^{3}(z_{1} + 1)} + \frac{z_{1}^{−300}}{(z_{1} − 1)^{2}(z_{1} + 1)(z_{3}^{3} − 1)}.

Again, the integral of the second term of I(z_{1}) along |z_{1}| = w_{1} is equal to zero because of (3.12), whereas the first term gives us

− \frac{7441}{8} + \frac{1}{8} = −930

Thus, we have have obtained g(120), that we also compare with the one obtained in (4.5) with t = 120, by the inverse associated transform algorithm.

f(120e_{3}) = 6141 − 930  [by direct inversion]
= \frac{17 ∗ 120^{2}}{48} + \frac{41 ∗ 120}{48} + \frac{139}{288} + \frac{120}{16} + \frac{9}{32} + \frac{1}{8} + \frac{1}{9} [by (4.5)].

5.2.2. Which ℤ-inverse should be used? As a final comparison between the direct inversion of the ℤ-transform and the inversion of the associated ℤtransform, we can point out the following facts:

i) To compute f(y) for y ∈ ℤ^{m} via the direct inversion of the ℤtransform F, we need not suppose that there exists an entry y_{1} ≠ 0 which divides every other entry. However, for practical efficiency of the algorithm, we are condamned to use the only integration technique (c) in §3.2.

ii) On the other hand, the inversion of the associated ℤ-transform gives us an explicit formula for g(t) = f(Dt) (under the above restriction on the entries of y). Moreover, it is likely to be more time efficient because at each integration step, and depending on the data on hand, we may choose between the three alternative integration techniques (a), (b) or (c) in §3.2. As a by-product, we also obtain the Ehrhart polynomial of an integer polytope.

5.2.3. Computational complexity. We have not succeeded in getting an exact evaluation of the computational complexity of both algorithms (direct inversion and inversion of the associated ℤ-transform). One important parameter to evaluate this complexity is the total number of poles that are considered at each of the m one-dimensional integration steps.

Basically, it turns out that in the worst case, the number of poles at step k is equal to ∑_{i_{1},...i_{k}} Δ^{1,2,...k}(i_{1}i_{2} … i_{k}), where

(5.1)  Δ^{1,2,...k}(i_{1}i_{2} … i_{k}) := det [ \begin{matrix}A_{1i_{1}} & A_{1i_{2}} & . & A_{1i_{k}} \\ A_{2i_{1}} & A_{2i_{2}} & . & A_{2i_{k}} \\ . & . & . & . \\ A_{ki_{1}} & A_{ki_{2}} & . & A_{ki_{k}} \end{matrix} ] ,

where i_{j} ∈ {1, … n} and the i_{j}’s are pairwise distinct. Therefore, as we have at most \binom{n}{k} such determinants, we roughly have at most O(Kn^{k}) poles with K := max_{i_{1},...i_{k}} |Δ^{1,2,...k}(i_{1}i_{2} … i_{k})| (which is a polynomial in the data). Therefore, the total number of poles is O(Mn^{m}) with M a polynomial in the data (exponential in the input size in computational complexity terminology). We also must evaluate the total work required in the integration techniques (a), (b) or (c) in §3.2, which is difficult because it involves expansions in fractions. Moreover, the extra work required for the eventual derivations needed in case of multiple poles may become complex very quickly. However, it is worth noticing that a careful choice of the vector c permits to avoid the occurence of multiple poles during the integration procedure, as shown in the next section.

6. An approximate inversion of the ℤ-transform

We have seen in section §3.2 that the integration procedure can become complex, particularly when multiple poles are encountered at some stage of the algorithm. Fortunately, one may avoid the occurrence of multiple poles by invoking an approximation of the cost function c in the expression of F(z).

The main idea of this algorithm is to calculate an approximate value of (2.5). That is, instead of considering a fixed value of the weighting vector c, we will calculate the inverse ℤ-transform of

(6.1)  F̃(z) = ∏_{k=1}^{n} \frac{1}{1 − e^{d_{k}}z_{1}^{−A_{1k}} ⋯ z_{m}^{−A_{mk}}} ,

where the vector d is allowed to lie in an open ball centered at c and with radius small enough.

In order to simplify the exposition let us assume that every entry of A is positive. Recall that we are calculating integral in equation (2.11) along the integration paths |z_{j}| = w_{j} where the fixed vector w ∈ ℝ_{+}^{m} satisfies A′ ln(w) > c. So fix a radius r > 0 small enough so as to ensure that A′ ln(w) > d for every d in the open ball B(c, r) ⊂ ℝ^{n}. Then, we calculate the integral

(6.2)  f̃(y) = \frac{1}{(2πi)^{m}} ∫_{|z_{1}|=w_{1}} ⋯ ∫_{|z_{m}|=w_{m}} F̃(z)z^{y−e_{m}}dz,

by choosing a vector d ∈ B(c, r) that makes calculations easier.

For example, we can always choose d ∈ B(c, r) in such a way that the real numbers:

|e^{d_{k}}z_{2}^{A_{2k}}z_{3}^{A_{3k}} ⋯ z_{m}^{A_{mk}}|^{1/A_{1k}}, k = 1, …, n

are pairwise distinct. This condition automatically implies that the poles of z_{1} ↦ F̃(z) are simple. Moreover, we can calculate the following decomposition :

(6.3)  F̃(z) = ∑_{k=1}^{n} \frac{Q_{k}(z_{1}^{−1})}{1 − e^{d_{k}}z_{1}^{−A_{1k}} ⋯ z_{m}^{−A_{mk}}} ,

for some real-valued polynomials

Q_{k}(z_{1}^{−1}) = ∑_{j=1}^{A_{1k}−1} Q_{kj}z_{1}^{−j}, k = 1, … n,

where each coefficient Q_{jk} is constant with respect to variable z_{1}. Hence :

F̃(z) = ∑_{k=1}^{n} ∑_{j=1}^{A_{1k}−1} ∑_{x∈ℕ}Q_{kj}e^{d_{k}x}z_{1}^{−j−A_{1k}x}z_{2}^{−A_{2k}x} ⋯ z_{m}^{−A_{mk}x}.

Integrating with respect to z_{1}, we obtain

\frac{1}{2πi} ∫_{|z_{1}|=w_{1}} F̃(z)z^{y−e_{m}}dz_{1} = ∑_{k=1}^{n} ∑_{j=1}^{A_{1k}−1} Q_{kj} ×

z_{2}^{y_{2}−A_{2k}x_{0}−1} ⋯ z_{m}^{y_{m}−A_{mk}x_{0}−1} \{ \begin{matrix}e^{d_{k}x_{0}} & if x_{0} = \frac{y−j}{A_{jk}} ∈ ℕ \\ 0 & otherwise. \end{matrix}

Notice that each Q_{jk} is a real-valued rational fraction of the variables z_{2}, z_{3}, … z_{n}. Therefore, we can write them as

Q_{jk} = P (z_{2}^{−1}) ∏_{i=1}^{o_{jk}} (1 − z_{2}^{−α_{i}}β_{i})^{−1} ,

where P is a real-valued polynomial of z_{2}^{−1}, each α_{i} ∈ ℕ and each β_{i} is constant with respect to z_{2} and a real-valued rational function of variables z_{3}, z_{4}, … z_{n}.

Thus, as we already did for z_{1}, we can choose again a vector d ∈ B(c, r) such that Q_{jk} has a decomposition as in (6.3), and repeat the same procedure until we obtain the value of f̃(y) for some vector d ∈ ℝ^{n} close enough to c.

7. An alternative method

In this section we propose an alternative method which avoids complex integration and is purely algebraic. It is particularly attractive for small values of m.

7.1. The method. With e_{m} = (1, 1, …) and doing the change of variables z = p^{−e_{m}} in the functions (2.8) and (2.9), we obtain

(7.1)  F(p^{−e_{m}}) = ∑_{y∈ℤ^{m}} f(y)p^{y} = ∏_{k=1}^{n} [\frac{1}{1 − e^{c_{k}}p_{1}^{A_{1k}}p_{2}^{A_{2k}} ⋯ p_{m}^{A_{mk}}} ] .

Partition the matrix A ∈ ℤ^{m×n} into its positive and negative parts A^{+}, A^{−} ∈ ℕ^{m×n}, defined by

A_{ij}^{+} := \{ \begin{matrix}A_{ij} & if A_{ij} ≥ 0, \\ 0 & otherwise, \end{matrix} and A^{−} := A^{+} − A.

The notation A_{k} stands for the k-th column of A for all k = 1, …, n. The same notation A_{k}^{+} and A_{k}^{−} applies to A^{+} and A^{−} respectively. Besides,

recall that p^{A_{k}^{+}} = p_{1}^{A_{1k}^{+}} ⋯ p_{m}^{A_{mk}^{+}} as in (2.1). Hence, we can rewrite F as

(7.2)  F(p^{−e_{m}}) = ∏_{k=1}^{n} [\frac{p^{A_{k}^{−}}}{p^{A_{k}^{−}} − e^{c_{k}}p^{A_{k}^{+}}} ] .

In order to simplify the exposition, we make the following weak assumption.

Assumption 7.1. Let n > m and o := n−m. The solutions s ∈ ℝ^{m} (whenever they exist) of the \binom{n}{m} systems of linear equations

(s_{1}, …, s_{m}) [A_{(j_{1})}|A_{(j_{2})}| ⋯ |A_{(j_{o})}] = (c_{(j_{1})}, … c_{(j_{o})}) ,

1 ≤ j_{1} < j_{2} < … < j_{o} ≤ n, are pairwise distinct.

Obviously, Assumption 7.1 implies in particular that the \binom{n}{m} polynomials

∏_{k≠j_{1}...j_{m}} (p^{A_{k}^{−}} − e^{c_{k}}p^{A_{k}^{+}}) 1 ≤ j_{1} < j_{2} < … < j_{m} ≤ n

have no common zero in ℂ^{m}.

Therefore, by Hilbert’s Nullstellensatz Theorem, there exist \binom{n}{m} realvalued polynomials Q_{j_{1},...j_{m}} ∈ ℝ[p_{1}, …, p_{m}] such that

(7.3)  1 = ∑_{1≤j_{1}<...<j_{m}≤n} Q_{j_{1}...j_{m}}(p) ∏_{k≠j_{1},...,j_{m}} (p^{A_{k}^{−}} − e^{b_{k}}p^{A_{k}^{+}}) .

Moreover, bounds on the degree of the Q_{j_{1},...,j_{m}} are available (see e.g. Seidenberg [17, Theorems 55-57, p.296-297] or Kollár [14]). For instance, in our present context, with r := \binom{n}{m}, let d_{1} ≥ d_{2} ≥ … ≥ d_{r} be the degrees of the polynomials

q_{j_{1},...,j_{m}}(p) := ∏_{k≠j_{1},...,j_{m}} (p^{A_{k}^{−}} − e^{c_{k}}p^{A_{k}^{+}}) , 1 ≤ j_{1} < … < j_{m} ≤ n.

Then, for all 1 ≤ j_{1} < … < j_{m} ≤ n, the degree of

Q_{j_{1}...j_{m}}(p)q_{j_{1},...j_{m}}(p),

in (7.3) is bounded by d_{1}d_{2} ⋯ d_{m−1}d_{r} (see Kollár [14]).

So one may compute the coefficients of the polynomials Q_{j_{1},...,j_{m}} by solving a linear system that matches the terms of same power in both sides of (7.3). In view of the above bound on the degree, this approach is realistic for relatively small values of m because the size of the resulting linear system grows rapidly with m. Indeed, a real-valued polynomial of degree r in m variables may have as many as \binom{m+r}{m} coefficients.

Multiplying together equations (7.2) and (7.3), we can deduce that

F_{1}(p^{−e_{m}}) = ∑_{1≤j_{1}<...<j_{m}≤n} \frac{Q_{j_{1}...j_{m}}(p) ∏_{k=1}^{n} p^{A_{k}^{−}}}{∏_{k=j_{1},...,j_{m}} (p^{A_{k}^{−}} − e^{c_{k}}p^{A_{k}^{+}})}

= ∑_{1≤j_{1}<...<j_{m}≤n} \frac{Q_{j_{1}...j_{m}}(p)p^{B(j_{1}...j_{m})}}{∏_{k=j_{1},...,j_{m}} (1 − e^{c_{k}}p^{A_{k}})},

where

B(j_{1} … j_{m}) := ∑_{k≠j_{1},...,j_{m}} A_{k}^{−}.

Notice that |e^{c_{k}}p^{A_{k}}| < 1 for k = 1, …, n, because of (2.10) and p = z^{−e_{m}}, so that we have the expansion

∏_{k=j_{1},...,j_{m}} \frac{1}{1 − e^{c_{k}}p^{A_{k}}} = ∑_{x∈ℤ^{m}} ∏_{k=1}^{m} e^{c_{j_{k}} x_{k}}p^{A_{(j_{k})}x_{k}}.

Moreover, writing

Q_{j_{1}...j_{m}}(p) = ∑_{β∈ℕ^{m},β≤M} Q_{j_{1}...j_{m}}^{(β)}p^{β},

for some constant bound M ∈ ℕ^{m}, we get that

(7.4)  F(p^{−e_{m}}) = ∑_{1≤j_{1}<...<j_{m}≤n} ∑_{{β∈ℕ^{m}, \atop β≤M}}  ∑_{x∈ℕ^{m}} Q_{j_{1}...j_{m}}^{(β)} ×

e^{∑_{k} c_{j_{k}} x_{k}} p^{β+B(j_{1}...j_{m})+A(j_{1}...j_{m})x},

with the square submatrices

A(j_{1} … j_{m}) := [A_{(j_{1})}|A_{(j_{2})}| ⋯ |A_{(j_{m})}]

for all 1 ≤ j_{1} < j_{2} … < j_{m} ≤ n.

Finally, notice that the sums in equations (7.1) and (7.4) are equal. Hence, if we want to deduce the exact value of f(y) from equation (7.4), we only have to sum up all the terms for which the exponent β+B(j_{1} … j_{m})+A(j_{1} … j_{m})x is equal to y. That is, given the condition

(7.5)  β_{x} := y − B(j_{1} … j_{m}) − A(j_{1} … j_{m})x ∈ [0, M]^{m},

we have

(7.6) f(y) = ∑_{1≤j_{1}<...<j_{m}≤n} ∑_{x∈ℕ^{m}} Q_{j_{1}...j_{m}}^{(β_{x})} \{ \begin{matrix}e^{∑_{k} c_{j_{k}} x_{k}} & if (7.5) holds, \\ 0 & otherwise. \end{matrix}

Moreover, with the additional assumption that all the square submatrices A(j_{1} … j_{m}) are non-singular, (7.6) simplifies, that is, considering the condition :

(7.7)  x_{β} := A(j_{1} … j_{m})^{−1} [y − β − B(j_{1} … j_{m})] ∈ ℕ^{m},

we now have

(7.8) f(y) = ∑_{1≤j_{1}<...<j_{m}≤n} ∑_{{β∈ℕ^{m}, \atop β≤M}} Q_{j_{1}...j_{m}}^{(β)} \{ \begin{matrix}e^{∑_{k} c_{j_{k}} [x_{β}]_{k}} & if (7.7) holds, \\ 0 & otherwise. \end{matrix}

Example 7.2. Consider the matrices A = \binom{1,2,1}{2,1,1} in ℝ^{2×3}, c = ln(2)(1, 1, 1) in ℝ^{3}, and the convex polytope

Ω(b) := {x ∈ ℝ^{3} | x_{1} + 2x_{2} + x_{3} = b_{1}, 2x_{1} + x_{2} + x_{3} = b_{2}, x ≥ 0}.

We obtain

F(p^{−e_{3}}) = \frac{1}{(1 − 2p_{1}p_{2}^{2})(1 − 2p_{1}^{2}p_{2})(1 − 2p_{1}p_{2})},

and

1 = Q_{23}(p)(1 − 2p_{1}p_{2}^{2}) + Q_{13}(p)(1 − 2p_{1}^{2}p_{2}) + Q_{12}(z)(1 − 2p_{1}p_{2}),

where

Q_{23}(p) = 1 + 2p_{1}^{2}p_{2}, Q_{13}(p) = 1 + 2p_{1}p_{2}^{2}
and Q_{12}(p) = −(1 + 2p_{1}p_{2} + 4p_{1}^{2}p_{2}^{2}).

Hence,

f(y) = ∑_{Ax=y,x∈ℕ^{3}} 2^{x_{1}+x_{2}+x_{3}}

= \{ \begin{matrix}2^{y_{2}} & if y_{1} − y_{2} ≥ 0 and 2y_{2} − y_{1} ≥ 0 \\ 0 & otherwise \end{matrix}

+ \{ \begin{matrix}2^{y_{2}} & if y_{1} − y_{2} − 1 ≥ 0 and 2y_{2} − y_{1} ≥ 0 \\ 0 & otherwise \end{matrix}

+ \{ \begin{matrix}2^{y_{1}} & if y_{2} − y_{1} ≥ 0 and 2y_{1} − y_{2} ≥ 0 \\ 0 & otherwise \end{matrix}

+ \{ \begin{matrix}2^{y_{1}} & if y_{2} − y_{1} − 1 ≥ 0 and 2y_{1} − y_{2} ≥ 0 \\ 0 & otherwise \end{matrix}

− \{ \begin{matrix}2^{(y_{1}+y_{2})/3} & if \frac{2y_{2}−y_{1}}{3} ∈ ℕ and \frac{2y_{1}−y_{2}}{3} ∈ ℕ \\ 0 & otherwise \end{matrix}

− \{ \begin{matrix}2^{(y_{1}+y_{2}+5)/3} & if \frac{2y_{2}−y_{1}+1}{3} ∈ ℕ and \frac{2y_{1}−y_{2}+1}{3} ∈ ℕ \\ 0 & otherwise \end{matrix}

− \{ \begin{matrix}2^{(y_{2}+y_{2}+10)/3} & if \frac{2y_{2}−y_{1}+2}{3} ∈ ℕ and \frac{2y_{1}−y_{2}+2}{3} ∈ ℕ \\ 0 & otherwise \end{matrix}

7.2. Discussion. It worth noticing that in this approach, most of the computational effort is concentrated on the determination of the coefficients of the polynomials ℚ_{j_{1},...,j_{m}} in the Hilbert NullStellensatz decomposition (7.3), which does not depend on the right-hand-side b of Ax = b; indeed b is only involved in the final (and easy) evaluation (7.7)-(7.8). Moreover, once the Q_{j_{1},...,j_{m}} are calculated by solving a linear system, one may evaluate f(b) for any value of b at a cheap cost (like Brion and Vergne’ s formula [8] which is valid for all b in the same chamber). As already mentioned, this approach is particularly attractive for relatively small values of m, otherwise the linear system to solve is too large.

At this stage, it is not clear whether this general method whose main work does not depend on b is preferrable to the more specific associated ℤ-inverse method that relies on Cauchy ’s Residue Theorem but where all calculations for b are useless for another b′ ≠ b.

8. Conclusion

We have presented an algorithm (in fact, three) for computing the number of nonnegative integral points in a convex rational polytope {x ∈ ℝ_{+}^{n} | Ax ≤ b}. They are all based on the inversion of some ℤ-transform by means of residues. In contrast to the algorithm proposed by Barvinok [3] which works in the space ℝ^{n} of primal variables x, we rather work in the space ℝ^{m} of dual variables z associated with the m nontrivial constraints Ax = b. In addition, we need not know the vertices of the polytope explicitly. As such, it provides an alternative method. Despite we have not completely characterized the computational complexity of the algorithm, it might work for potentially large values of n and relatively small values of m, a context “dual” to that of Barvinok’s algorithm which is polynomial in the problem size for fixed dimension n. Finally the (algebraic) alternative approach described in §7 does not use residues “directly” and is particularly attractive for relatively small values of m.

References

[1] W. Baldoni-Silva, M. Vergne. Residues formulae for volumes and Ehrhart polynomials of convex polytopes, arXiv:math.CO/0103097 v1, (2001).

[2] A.I. Barvinok. A polynomial time algorithm for counting integral points in polyhedra when the dimension is fixed, Math. Oper. Res. 19 (1994), 769–779.

[3] A.I. Barvinok, J.E. Pommersheim. An algorithmic theory of lattice points in polyhedra, in New Perspectives in Algebraic Combinatorics, MSRI Publications 38 (1999), 91–147.

[4] M. Beck. Multidimensional Ehrhart reciprocity, J. Comb. Theory Ser. A 97 (2002), 187–194.

[5] M. Beck. Counting Lattice Points by means of the Residue Theorem. Ramanujan Journal 4 (2000), 399–310.

[6] M. Beck, R. Diaz and S. Robins. The Frobenius problem, rational polytopes, and Fourier-Dedekind sums, J. Numb. Theor., to appear.

[7] M. Brion. Points entiers dans les polyèdres convexes, Ann. Ecol. Norm. Sup. (Sér. 4) 21 (1988), 653–663.

[8] M. Brion, M. Vergne. Residue formulae, vector partition functions and lattice points in rational polytopes, J. Amer. Math. Soc. 10 (1997), 797–833.

[9] J.B. Conway. Functions of a complex variable I, 2nd ed., Springer, New York, (1978).

[10] E. Ehrhart. Sur un problème de géométrie diophantienne linéaire II, J. Reine Angew. Math. 227 (1967), 25–49.

[11] J.-B. Hiriart-Urruty, C. Lemarechal. Convex Analysis and Minimization Algorithms I, Springer-Verlag, Berlin, 1993.

[12] J-M. Kantor, A, Khovanskii. Une application du théorème de Riemann-Roch combinatoire au polynôme d’Ehrhart des polytopes entiers, C.R. Acad. Sci. Paris (Série I) 317 (1993), 501–507.

[13] A. Khovanskii, A. Pukhlikov. A Riemann-Roch theorem for integrals and sums of quasipolynomials over virtual polytopes, St-Petersburg Math. J. 4 (1993), 789–812.

[14] J. Kollár. Sharp effective Nullstellensatz, J. Am. Math. Soc. 1 (1988), 963–975.

[15] J.B. Lasserre, E.S. Zeron. A Laplace transform algorithm for the volume of a convex polytope, J. ACM 48 (2001), 1126–1140.

[16] A. Schrijver. Theory of Linear and Integer Programming, John Wiley & Sons, Chichester, 1986.

[17] A. Seidenberg. Constructions in algebra, Trans. Am. Math. Soc. 197 (1974), 273–313.

LAAS-CNRS, 7 Avenue du Colonel Roche, 31077 Toulouse Cédex 4, France.

E-mail address: lasserre@laas.fr

Depto. Matemáticas, CIVESTAV-IPN, Apdo. Postal 14740, Mexico D.F. 07000, México.

E-mail address: eszeron@math.cinvestav.mx