Counting integer points in polytopes is a fundamental problem. There are numerous applications in various areas of mathematical science, and fascinating mathematics behind; see e.g., [2, 4]. This problem is computationally intractable, i.e., it is P-hard . Approximate counting as well as exact counting under fixed parameter settings has been rich sources for developments in the theory of algorithms and computational complexity. A seminal work by Barvinok  showed that there is a polynomial time algorithm to count integer points in rational polytope when the dimension of is fixed. His algorithm computes a certain “compact” expression (Brion-Lawrence formula) of the generation function for multivariate indeterminate . This (extremely difficult) algorithm is now implemented in computer package LattE , and provides a useful tool to the study of geometric combinatorics.
for integer matrix and -dimensional vector . Let , and consider the Z-transform . Brion-Vergne  showed that admits a very simple closed formula , and that the wanted is recovered by the inverse Z-transformation, which is a multi-dimensional contour integration of . This reduces the counting problem to the residue computation of . By this approach, Lasserre-Zeron [13, 14] developed an -time algorithm to count integral points in , where is a function of matrix .
In this paper, we study the contour integration of the inverse Z-transformation from the numerical analysis point of view, and obtain a new algorithm to count integer points for an important class of polytopes. Our main result is as follows.
Suppose that is nonnegative. For , the number of integer points in the polytope can be computed in time and space.
Notice that there is a simple DP algorithm with the same time complexity: For , consider the matrix consisting of the first columns of , and the number of integer points of polyhedron . By the nonnegativity of , integers for are obtained from for in time. The resulting DP algorithm, however, requires an -space for the DP table. Thus our result is regarded as an improvement in terms of space complexity.
Our technique for proving Theorem 1 is as follows. Instead of the residue computation of , we apply the numerical integration to the inverse Z-transform of , where we use the trapezoidal rule, a basic and popular method of numerical integration. By the standard error analysis of the trapezoidal rule 
, we obtain an error estimate with respect to the numberof sampled points of the numerical integration and contour radius . Interestingly this estimate gives rise to a new inclusion-exclusion type formula for , and brings the algorithm in Theorem 1, which is quite simple and is easier to be implemented. A notable feature of our inclusion-exclusion is to use the cancellation structure of trigonometric function in the complex plane . This extends the usual inclusion-exclusion based on the cancellation of and in . Our algorithm computes the number of integer points in the expression for , and recovers the “true” value by algebraic computation on the group ring of cyclic group , which avoids numerical computation of .
Our result is applicable to packing-type polytopes, which ubiquitously arise from graph theory and combinatorial mathematics. Consider the particular case where is 0-1 valued and is the all-one vector . Then is the number of perfect matchings in the hypergraph corresponding to . Exact exponential time counting algorithms for matchings have been intensively studied in recent years [5, 6, 9, 11]. The formula of hypergraph matching derived from our result is viewed as a variant of that given by Björklund-Husfeldt . The time complexity matches that of their algorithm. By exploiting special properties, Björklund-Husfeldt  improve the time complexity to for -partite hypergraphs. This result is viewed as a multipartite extension of the classical result of Ryser  for bipartite matching. We also exploit a special structure of our formula in -partite hypergraphs, and prove the following -matching generalization of the results of Ryser and Björklund-Husfeldt.
Let be a -partite hypergraph, and let . The number of perfect -matchings in can be computed in time and space.
Counting perfect -matchings of a -partite graph (i.e., -dimensional -matchings) has many applications in a wide range of mathematical sciences that include combinatorics, representation theory, and statistics; see e.g., 
. For example, counting multiway contingency tables with prescribed margins, an important problem for statistical analysis on contingency tables, is nothing but-dimensional -matching counting.
The rest of this paper is organized as follows. In Section 2, we set up basic notation, and introduce Z-transformation, its inverse, and approximate inverse Z-transformation obtained by numerical integration. In Section 3, we present our algorithm to prove the main theorem. In Section 4, we discuss hypergraph matching and prove Theorem 2 in a further generalized form.
Let denote the set of positive real numbers. Let denote the set of nonnegative integers. For a matrix , let denote the -th column vector of the matrix. For an integer vector and a complex vector , define by
For a positive integer , define by
The following relation is well-known:
For a function and , let denote , such as
The symbol is meant as “undefined.” We use when the function value is defined via integration or infinite summation, possibly not converging.
For a function , define the Z-transform of by
The inverse of the Z-transformation is given as follows. For a function and , define by
where we change variables by in (9). Under an appropriate condition on and , map is actually the inverse of the Z-transformation:
We do not go into details under which conditions (10) holds. Instead, we consider an approximate inverse Z-transform by the numerical integration applied to (9). Here we use the trapezoidal rule, which is a basic and popular method of numerical integration; see e.g., . For a positive integer (the number of points in the numerical integration), define by
Recall notation . Our counting algorithm is based on .
3 Counting integral points in a polytope
Let be an integral matrix. We assume that there is no nonzero nonnegative vector with . This assumption ensures that the polytope is bounded for every . In the case where is a nonnegative matrix, this assumption is equivalent to the property that each column of has at least one nonzero entry.
As mentioned in the introduction, define function by
Our starting point is the following formula for .
For with , the Z-transform is given by
where the series absolutely converges.
For with , it holds
We establish an approximate version of the above theorem as follows.
Suppose that is nonnegative. For , , and , it holds
From the assumption that is nonnegative and each column of has nonzero entry, for with , it holds . By the previous theorem, the Z-transform is given by
Substituting this expression to (in (11)), we have
where the summations are interchangeable, thanks to the absolute convergence (by ). By the relation (3), is if , and zero otherwise. Notice that cannot occur since . Thus we have
Gathering with , we obtain (15). ∎
This proof is inspired by the standard argument to derive the exponential convergence of the trapezoidal rule applied to periodic functions, where the the cancellation technique using (3) in the proof is known as aliasing in numerical analysis; see [17, Theorem 2.1].
Notice that is a variant of the Ehrhart (quasi)polynomial, and is bounded by a polynomial in with degree . Hence the series in (15) actually absolutely converges for .
Suppose that is nonnegative. Let and . Then is equal to the coefficient of in
Here we regard as an indeterminate.
In the formula (14), if has , then for some , and does not contribute to the constant term . The claim follows from this fact and the observation . ∎
Our goal is to compute the coefficient of in (21). Instead of numerical computation of the trigonometric function , we develop an algebraic algorithm. Let denote the group ring of the cyclic group of order . Namely, consists of polynomials with variable , rational coefficients, and degree at most , in which the multiplication rule is given by . Consider the bivariate polynomial ring with variables . Then is a ring homomorphism from to . Letting and in (20), we obtain a polynomial in :
Let denote the coefficient of in (22). Then it holds
Our algorithm first computes , and then computes .
can be computed in time and space.
Let . From (21), we first consider the computation of the coefficient of in
for fixed . It suffices to compute the above polynomial (24) modulo , which can be written as
for some . The integers are obtained in time by computing and for . Expand (25) to the form , and obtain . This computation is done by the multiplication of polynomials with degree modulo , where their coefficients are polynomials with degree in multiplication rule . Now is the sum of over divided by . Thus is obtained in arithmetic operations (over ).
Finally we estimate the bit-size required for the computation. It suffices to estimate the size of , which is the sum of at most terms of form . Then the coefficients of have bit-length . Thus the required bit-size is at most . ∎
Next we consider how to compute from .
is written as
where the sum is taken over divisors of with some coefficient .
Regard as . Then the map defined by is a group homomorphism. Therefore the image of is the cyclic group for some divisor of . Also the number of inverse images of each is given by . Then, for with , it holds
Thus we have
where is the sum of over such that the image of is . Notice . ∎
According to this lemma, we obtain a simple algorithm to compute from as follows.
If for all , then output ; stop
Choose the minimum index with . Let for each index that is the multiple of the index , and go to step 1.
The correctness of the algorithm is clear from the above lemma: The chosen index in step 2 is a divisor of with . Hence the algorithm computes . After at most iterations, the algorithm terminates and outputs the correct answer .
Consider the following matrix and vector :
Then the polytope has three integer points:
Let us count the integer points according to Corollary 5. The coefficient of in
Therefore we obtain .
Consider the following matrix and vector :
Then the polytope has only one integer point:
Consider the coefficient of in
By calculation, is
Thus we obtain .
4 Hypergraph Matching
We next show Theorem 2 in a generalized form. Let be an nonnegative integer matrix. For each column index , consider subset consisting of row indices with . Let denote the hypergraph on vertex set and hyperedge set . By a stable set of we mean a vertex subset such that every hyperedge meets at most one vertex in .
Suppose that we are given a stable set of . For , we can compute in time and space.
where each is computable modulo in time and space.
By arranging indices of , we can assume that , and that is regarded as a block matrix , where consists of columns such that the corresponding hyperedge meets . Accordingly, vector is also partitioned as
We suppose that an -dimensional vector is embedded to by filling to the first components. Each is uniquely represented as for , where is the -th unit vector. Then we have