1 Introduction
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 Phard [18]. 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 [1] 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 (BrionLawrence formula) of the generation function for multivariate indeterminate . This (extremely difficult) algorithm is now implemented in computer package LattE [15], and provides a useful tool to the study of geometric combinatorics.
A “dual” generatingfunction approach was initiated by BrionVergne [8], Beck [3], and LasserreZeron [13, 14]; see [12]. Suppose now that the input polytope is given by
(1) 
for integer matrix and dimensional vector . Let , and consider the Ztransform . BrionVergne [8] showed that admits a very simple closed formula , and that the wanted is recovered by the inverse Ztransformation, which is a multidimensional contour integration of . This reduces the counting problem to the residue computation of . By this approach, LasserreZeron [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 Ztransformation 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.
Theorem 1.
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 Ztransform of , where we use the trapezoidal rule, a basic and popular method of numerical integration. By the standard error analysis of the trapezoidal rule [17]
, we obtain an error estimate with respect to the number
of sampled points of the numerical integration and contour radius . Interestingly this estimate gives rise to a new inclusionexclusion 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 inclusionexclusion is to use the cancellation structure of trigonometric function in the complex plane . This extends the usual inclusionexclusion 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 packingtype polytopes, which ubiquitously arise from graph theory and combinatorial mathematics. Consider the particular case where is 01 valued and is the allone 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örklundHusfeldt [6]. The time complexity matches that of their algorithm. By exploiting special properties, BjörklundHusfeldt [6] improve the time complexity to for partite hypergraphs. This result is viewed as a multipartite extension of the classical result of Ryser [16] 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örklundHusfeldt.
Theorem 2.
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., [10]
. 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 Ztransformation, its inverse, and approximate inverse Ztransformation 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.
2 Preliminaries
2.1 Notation
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
(2) 
The following relation is wellknown:
(3) 
For a function and , let denote , such as
(4)  
(5)  
(6) 
The symbol is meant as “undefined.” We use when the function value is defined via integration or infinite summation, possibly not converging.
2.2 Ztransformation
For a function , define the Ztransform of by
(7) 
The inverse of the Ztransformation is given as follows. For a function and , define by
(8)  
(9) 
where we change variables by in (9). Under an appropriate condition on and , map is actually the inverse of the Ztransformation:
(10) 
We do not go into details under which conditions (10) holds. Instead, we consider an approximate inverse Ztransform 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., [7]. For a positive integer (the number of points in the numerical integration), define by
(11) 
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 .
Theorem 3 ([3, 8, 13]).

For with , the Ztransform is given by
(12) where the series absolutely converges.

For with , it holds
(13)
We establish an approximate version of the above theorem as follows.
Theorem 4.
Suppose that is nonnegative. For , , and , it holds
(14)  
(15) 
Proof.
From the assumption that is nonnegative and each column of has nonzero entry, for with , it holds . By the previous theorem, the Ztransform is given by
(16) 
Substituting this expression to (in (11)), we have
(17)  
(18)  
(19) 
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 .
Corollary 5.
Suppose that is nonnegative. Let and . Then is equal to the coefficient of in
(20)  
(21) 
Here we regard as an indeterminate.
Proof.
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 :
(22) 
Let denote the coefficient of in (22). Then it holds
(23) 
Our algorithm first computes , and then computes .
Lemma 6.
can be computed in time and space.
Proof.
Let . From (21), we first consider the computation of the coefficient of in
(24) 
for fixed . It suffices to compute the above polynomial (24) modulo , which can be written as
(25) 
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 bitsize 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 bitlength . Thus the required bitsize is at most . ∎
Next we consider how to compute from .
Lemma 7.
is written as
where the sum is taken over divisors of with some coefficient .
Proof.
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
(26)  
(27) 
Thus we have
(28)  
(29) 
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.

Let .

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 .
Example 8.
Consider the following matrix and vector :
(30) 
Then the polytope has three integer points:
(31) 
Let us count the integer points according to Corollary 5. The coefficient of in
(32) 
is
(33)  
(34)  
(35) 
Therefore we obtain .
Example 9.
Consider the following matrix and vector :
(36) 
Then the polytope has only one integer point:
(37) 
Consider the coefficient of in
(38) 
By calculation, is
(39)  
(40)  
(41)  
(42) 
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 .
Theorem 10.
Suppose that we are given a stable set of . For , we can compute in time and space.
Proof.
Let , , and . As before, it suffices to compute (22) modulo . We show that (22) admits the following factorization:
(43) 
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
Comments
There are no comments yet.