V Baldoni
email: baldoni@mat.uniroma2.it

Nicole Berline
email: nicole.berline@math.polytechnique.fr

Jesús A De Loera
email: deloera@math.ucdavis.edu

Brandon E Dutra
email: bedutra@ucdavis.edu

Matthias Koeppe
email: mkoeppe@math.ucdavis.edu

Michele Vergne
email: vergne@math.jussieu.fr

O F Coeffi- Cients

Denumerant 2014 Sylvester's

Hal-01087813

M K Öppe

M Vergne

Matthias Köppe

Michèle Vergne

à la diffusion de documents scientifiques de niveau

Keywords: Denumerants, Ehrhart quasi-polynomials, restricted partitions, asymptotic behavior, polynomial-time algorithms

recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Let α = [α 1 , α 2 , . . . , α N , α N +1 ] be a sequence of positive integers. If t is a non-negative integer, we denote by E(α)(t) the number of solutions in non-negative integers of the equation N +1 i=1 α i x i = t. In other words, E(α)(t) is the same as the number of partitions of the number t using the parts α 1 , α 2 , . . . , α N , α N +1 (with repetitions allowed). Let us begin with some background and history before stating the precise results:

The combinatorial function E(α)(t) was called by J. Sylvester the denumerant. The denumerant E(α)(t) has a beautiful structure: it has been known since the times of Cayley and Sylvester that E(α)(t) is in fact a quasi-polynomial, i.e., it can be written in the form E(α)(t) = N i=0 E i (t)t i , where E i (t) is a periodic function of t (a more precise description of the periods of the coefficients E i (t) will be given later). In other words, there exists a positive integer Q such that for t in the coset q + QZ, the function E(α)(t) coincides with a polynomial function of t. This paper presents a new algorithm to compute individual coefficients of this function and uncovers new structure in generating functions that allows one to compute their periodicity. The study of the coefficients E i (t), in particular determining their periodicity, is a problem that has occupied various authors and it is the key focus of our investigations here. Sylvester and Cayley first showed that the coefficients E i (t) are periodic functions having period equal to the least common multiple of α 1 , . . . , α N +1 (see [START_REF] Beck | The polynomial part of a restricted partition function related to the Frobenius problem[END_REF][START_REF] Bell | Interpolated denumerants and Lambert series[END_REF] and references therein). In 1943, E. T. Bell gave a simpler proof and remarked that the period Q is in the worst case given by the least common multiple of the α i , but in general it can be smaller. A classical observation that goes back to I. Schur is that when the list α consist of relatively prime numbers, then asymptotically

E(α)(t) ≈ t N N ! α 1 α 2 • • • α N +1
as the number t → ∞.

Thus, in particular, there is a large enough integer F such that for any t ≥ F , E(α)(t) > 0 and there is a largest t for which E(α)(t) = 0. Let us give a simple example:

Example 1.1. Let α = [START_REF] Baldoni | Intermediate sums on polyhedra: Computation and real Ehrhart theory[END_REF][START_REF] Agnarsson | On the Sylvester denumerants for general restricted partitions[END_REF][START_REF] Andrews | The theory of partitions[END_REF]. Then on each of the cosets q + 6Z, the function E(α)(t) coincides with a polynomial E [q] (t). Here are the corresponding polynomials.

E [0] (t) = 1 72 t 2 + 1 4 t + 1, E [1] (t) = 1 72 t 2 + 1 18 t -5 72 , E [2] (t) = 1 72 t 2 + 7 36 t + 5 9 , E [3] (t) = 1 72 t 2 + 1 6 t + 3 8 , E [4] (t) = 1 72 t 2 + 5 36 t + 2 9 , E [5] (t) = 1 72 t 2 + 1 9 t + 7 72 . Naturally, the function E(α)(t) is equal to 0 if t does not belong to the lattice N +1 i=1 Zα i ⊂ Z generated by the integers α i . So if g is the greatest common divisor of the α i (which can be computed in polynomial time), and α/g = [ α 1 g , α 2 g , . . . , α N +1 g ] the formula E(α)(gt) = E(α/g)(t) holds, and we may assume that the numbers α i span Z without changing the complexity of the problem. In other words, we may assume that the greatest common divisor of the α i is equal to 1.

Our primary concern is how to compute E(α)(t), a problem has received a lot of attention. Computing the denumerant E(α)(t) as a close formula or evaluating it for specific t is relevant in several other areas of mathematics. In the combinatorics literature the denumerant has been studied extensively (see e.g., [START_REF] Agnarsson | On the Sylvester denumerants for general restricted partitions[END_REF][START_REF] Beck | The polynomial part of a restricted partition function related to the Frobenius problem[END_REF][START_REF] Comtet | The art of finite and infinite expansions[END_REF][START_REF] Lisoněk | Denumerants and their approximations[END_REF][START_REF] Riordan | An introduction to combinatorial analysis[END_REF] and the references therein). The denumerant plays an important role in integer optimization too [START_REF] Kellerer | Knapsack problems[END_REF][START_REF] Martello | Knapsack problems[END_REF], where the problem is called an equality-constrained knapsack. In combinatorial number theory and the theory of partitions, the problem appears in relation to the Frobenius problem or the coin-change problem of finding the largest value of t with E(α)(t) = 0 (see [START_REF] Einstein | Frobenius numbers by lattice point enumeration[END_REF][START_REF] Kannan | Lattice translates of a polytope and the Frobenius problem[END_REF][START_REF] Ramírez Alfonsín | The Diophantine Frobenius problem[END_REF] for details and algorithms). Authors in the theory of numerical semigroups have also investigated the so called gaps or holes of the function (see [START_REF] Hemmecke | Computing holes in semi-groups and its application to transportation problems[END_REF] and references therein), which are values of t for which E(α)(t) = 0, i.e., those positive integers t which cannot be represented by the α i . For N = 1 the number of gaps is (α 1 -1)(α 2 -1)/2 but for larger N the problem is quite difficult.

Unfortunately, computing E(α)(t) or evaluating it are very challenging computational problems. Even deciding whether E(α)(t) > 0 for a given t, is a well-known (weakly) NPhard problem. Computing E(α)(t), i.e., determining the number of solutions for a given t, is #P -hard. Computing the Frobenius number is also known to be NP-hard [START_REF] Ramírez Alfonsín | The Diophantine Frobenius problem[END_REF]. Likewise, for a given coset q + QZ, computing the polynomial E [q] (t) is NP-hard. Despite the difficulty to compute the function, in some special cases one can compute information efficiently. For example, the Frobenius number can be computed in polynomial time when N + 1 is fixed [START_REF] Kannan | Lattice translates of a polytope and the Frobenius problem[END_REF][START_REF] Barvinok | Short rational generating functions for lattice point problems[END_REF]. At the same time for fixed N + 1 one can compute the entire quasi-polynomial E(α)(t) in polynomial time as a special case of a well-known result of Barvinok [START_REF] Barvinok | Polynomial time algorithm for counting integral points in polyhedra when the dimension is fixed[END_REF]. There are several papers exploring the practical computation of the Frobenius numbers (see e.g., [START_REF] Einstein | Frobenius numbers by lattice point enumeration[END_REF] and the many references therein).

We are certainly not the first to use generating functions to compute E(α)(t). Already Ehrhart obtained formulas for E(α)(t) in terms of binomial coefficients using partial fraction decomposition. Similary, in [START_REF] Sills | Formulae for the number of partitions of n into at most m parts (using the quasi-polynomial ansatz)[END_REF] the authors propose another way to recover the coefficients of the quasi-polynomial by a method they named rigorous guessing. In [START_REF] Sills | Formulae for the number of partitions of n into at most m parts (using the quasi-polynomial ansatz)[END_REF] quasi-polynomials are represented as a function f (t) given by q polynomials f [1] (t), f [2] (t), . . . , f [q] (t) such that f (t) = f [i] (t) when t ≡ i (mod q). To find the coefficients of the f [i] their method finds the first few terms of the Maclaurin expansion of the partial fraction decomposition to find enough evaluations of those polynomials and then recovers the coefficients of the f [i] as a result of solving a linear system. Here we are able to prove good complexity results and produced faster practical algorithms using the number-theoretic nature of the question.

It should be noted that the polynomial-time complexity results for fixed N were achieved using a powerful geometric interpretation of E(α)(t) (which was the original way we encountered the problem too): The function E(α)(t) can also be thought of as the number of integral points in the

N -dimensional simplex in R N +1 defined by ∆ α = { [x 1 , x 2 , . . . , x N , x N +1 ] : x i ≥ 0, N +1 i=1 α i x i = t } with rational vertices s i = [0, . . . , 0, t α i , 0, . . . , 0]
. In this context, E(α)(t) is a very special case of the Ehrhart function (in honor of French mathematician Eugène Ehrhart who started its study [START_REF] Ehrhart | Polynômes arithmétiques et méthode des polyèdres en combinatoire[END_REF]). Ehrhart functions count the lattice points inside a convex polytope P as it is dilated t times. All of the results we mentioned about E(α)(t) are in fact special cases of theorems from Ehrhart theory [START_REF]Integer points in polyhedra[END_REF]. For example, the asymptotic result of I. Schur can be recovered from seeing that the highest-degree coefficient of E α (t) is just the normalized N -dimensional volume of the simplex ∆ α . Our coefficients are very special cases of Ehrhart coefficients.

This paper is about the computation of the coefficients of E(α)(t). Here are our main results:

(1) It is clear that the leading coefficient is given by Schur's result. Our main result is a new algorithm for computing explicit formulas for more coefficients.

Theorem 1.2. Given any fixed integer k, there is a polynomial time algorithm to compute the highest k + 1 degree terms of the quasi-polynomial E(α)(t), that is

Top k E(α)(t) = k i=0 E N -i (t)t N -i .
The coefficients are recovered as step polynomial functions of t.

Note that the number Q of cosets for E(α)(t) can be exponential in the binary encoding size of the problem, and thus it is impossible to list, in polynomial time, the polynomials E [q] (t) for all the cosets q + QZ. That is why to obtain a polynomial time algorithm, the output is presented in the format of step polynomials, which we now introduce:

(i) We first define the function {s} = s -s ∈ [0, 1) for s ∈ R, where s denotes the largest integer smaller or equal to s. The function {s+1} = {s} is a periodic function of s modulo 1. (ii) If r is rational with denominator q, the function T → {rT } is a function of T ∈ R periodic modulo q. A function of the form T → i c i {r i T } will be called a (rational) step linear function. If all the r i have a common denominator q, this function is periodic modulo q.

(iii) Then consider the algebra generated over Q by such functions on R. An element φ of this algebra can be written (not in a unique way) as

φ(T ) = L l=1 c l J l j=1 {r l,j T } n l,j .
Such a function φ(T ) will be called a (rational) step polynomial. (iv) We will say that the step polynomial φ is of degree (at most) u if j n l,j ≤ u for each index l occurring in the formula for φ. 1 We will say that φ is of period q if all the rational numbers r j have common denominator q. In Example 1.1, instead of the Q = 6 polynomials E [0] (t), . . . , E [5] (t) that we wrote down, we would write a single closed formula, where the coefficients of powers of t are step polynomials in t:

1 72 t 2 + 1 4 - {-t 3 } 6 - { t 2 } 6 t+ 1 - 3 2 {-t 3 } - 3 2 { t 2 } + 1 2 {-t 3 } 2 + {-t 3 }{ t 2 } + 1 2 { t 2 } 2 .
For larger Q, one can see that this step polynomial representation is much more economical than writing the individual polynomials for each of the cosets of the period Q.

Our results come after an earlier result of Barvinok [9] who first proved a similar theorem valid for all simplices. Also in [START_REF] Baldoni | Computation of the highest coefficients of weighted Ehrhart quasi-polynomials of rational polyhedra[END_REF], the authors presented a polynomialtime algorithm to compute the coefficient functions of Top k E(P )(t) for any simple polytope P (given by its rational vertices) in the form of step polynomials defined as above. We note that both of these earlier papers use the geometry of the problem very strongly; instead our new algorithm is different as it uses more of the number-theoretic structure of the special case at hand. There is a marked advantage of our algorithms over the work in [START_REF]Computing the Ehrhart quasi-polynomial of a rational simplex[END_REF]: We compute in a closed formula using the step polynomials all the possibilities of E [q] (t) while [START_REF]Computing the Ehrhart quasi-polynomial of a rational simplex[END_REF] recovers a single polynomial E [q] (t) for a given q. More importantly, our new algorithm is much easier to implement. Another relevant prior work (also useful for comparison) is our algorithm LattE Top-Ehrhart presented in [START_REF] Baldoni | Computation of the highest coefficients of weighted Ehrhart quasi-polynomials of rational polyhedra[END_REF]. In that paper we extend Barvinok's results of [START_REF]Computing the Ehrhart quasi-polynomial of a rational simplex[END_REF] to weighted Ehrhart quasipolynomials via variation of his original approach. The other important ingredient used in the efficient computation of the top coefficients is the reinterpretation of some generating functions in terms of lattice points in cones. This allows us to apply the polynomial-time signed cone decomposition of Barvinok for simplicial cones of fixed dimension k [START_REF] Barvinok | Polynomial time algorithm for counting integral points in polyhedra when the dimension is fixed[END_REF].

(2) Although the main result is computational, interesting mathematics comes into play:

the new algorithm uses directly the residue theorem in one complex variable, which can be applied more efficiently as a consequence of a rich poset structure on the set of poles of the associated rational generating function for E(α)(t) (see Subsection 2.3). By Schur's result, it is clear that the coefficient E N (t) of the highest degree term is just an explicit constant. Our analysis of the high-order poles of the generating function associated to E(α)(t) allows us to decide what is the highest-degree coefficient of E(α)(t) that is not a constant function of t (we will also say that the coefficient is strictly periodic).

Theorem 1.3. Given a list of non-negative integer numbers α = [α 1 , . . . , α N +1 ], let be the greatest integer for which there exists a sublist α J with |J| = , such that its greatest common divisor is not 1. Then for k ≥ the coefficient of degree k is a constant while the coefficient of degree -1 of the quasi-polynomial E(α)(t) is strictly periodic. Moreover, if the numbers α i are given with their prime factorization, then detecting can be done in polynomial time.

Example 1.4. We apply the theorem above to investigate the question of periodicity of the denumerant coefficients in the case of the classical partition problem E([1, 2, 3, . . . , m])(t). It is well known that this coincides with the classical problem of finding the number of partitions of the integer t into at most m parts, usually denoted p m (t) (see [START_REF] Andrews | The theory of partitions[END_REF]). In this case, Theorem 1.3 predicts indeed that the highest-degree coefficient of the partition function p m (t) which is non-constant is the coefficient of the term of degree m/2 . This follows from the theorem because the even numbers in the set {1, 2, 3, . . . , m} form the largest sublist with gcd two.

(3) The paper closes with an extensive collection of computational experiments (Section 5). We constructed a dataset of over 760 knapsacks and show our new algorithm is the fastest available method for computing the top k terms in the Ehrhart quasipolynomial. Our implementation of the new algorithm is made available as a part of the free software LattE integrale [START_REF] Baldoni | A user's guide for LattE integrale v1.7.2[END_REF], version 1.7.2.2

The residue formula for E(α)(t)

Let us begin fixing some notation. If φ(z) dz is a meromorphic one form on C, with a pole at z = ζ, we write

Res z=ζ φ(z) dz = 1 2πi C ζ φ(z) dz,
where C ζ is a small circle around the pole ζ. If φ(z) = k≥k 0 φ k z k is a Laurent series in z, we denote by res z=0 the coefficient of z -1 of φ(z). Cauchy's formula implies that res z=0 φ(z) = Res z=0 φ(z) dz.

A residue formula for

E(α)(t). Let α = [α 1 , α 2 , . . . , α N +1 ] be a list of integers. Define F (α)(z) := 1 N +1 i=1 (1 -z α i )
.

Denote by Note that because the α i have greatest common divisor 1, we have ζ = 1 as a pole of order N + 1, and the other poles have order strictly smaller.

P = N +1 i=1 { ζ ∈ C : ζ α i = 1 }
Theorem 2.1. Let α = [α 1 , α 2 , . . . , α N +1 ] be a list of integers with greatest common divisor equal to 1, and let

F (α)(z) := 1 N +1 i=1 (1 -z α i ) .
If t is a non-negative integer, then

E(α)(t) = - ζ∈P Res z=ζ z -t-1 F (α)(z) dz (2.1)
and the ζ-term of this sum is a quasi-polynomial function of t with degree less than or equal to p(ζ) -1.

Proof. For |z| < 1, we write 1 1-z α i = ∞ u=0 z uα i so that F (α)(z) = t≥0 E(α)(t)z t . For a small circle |z| = of radius around 0, the integral of z k dz is equal to 0 except if k = -1, when it is 2πi. Thus E(α)(t) = 1 2πi |z|= z -t F (α)(z) dz z = 1 2πi |z|= z -t N +1 i=1 1 (1 -z α i ) dz z .
Because the α i are positive integers, and t a non-negative integer, there are no residues at z = ∞ and we obtain Equation (2.1) by applying the residue theorem (for a reference about computational complex analysis see [START_REF] Henrici | Applied and computational complex analysis[END_REF][START_REF]Power series-integration-conformal mapping-location of zeros[END_REF][START_REF]Special functions-integral transforms-asymptotics-continued fractions[END_REF].

) Write E ζ (t) := -Res z=ζ z -t F (α)(z) dz z ; then the dependence in t of E ζ (t) comes from the expansion of z -t near z = ζ. We write z = ζ + y, so that E ζ (t) = -Res y=0 (ζ + y) -t F (α)(ζ + y) dy ζ + y .
As the pole of F (α)(ζ + y) at y = 0 is of order p(ζ), to compute the residue at y = 0, we only need to expand in y the function (ζ + y) -t-1 and take the coefficient of y p(ζ)-1 . Now from the generalized Newton binomial theorem, for

k = t + 1 the function (ζ + y) -k = ∞ n=0 n+k-1 n ζ -k-n (-y) n
. From this expression one can recover the desired coefficient. One can easily check that the dependence in t of our residue is a quasi-polynomial with degree less than or equal to p(ζ) -1. We thus obtain the result.

2.2.

Poles of high and low order. Given an integer 0 ≤ k ≤ N , we partition the set of poles P in two disjoint sets according to the order of the pole:

P >N -k = { ζ : p(ζ) ≥ N + 1 -k }, P ≤N -k = { ζ : p(ζ) ≤ N -k }.
ζ 2 6 , ζ 3 6 = -1, ζ 4
6 are poles of order 3. Thus

P >N -k = P >2 is the union of { ζ : ζ 2 = 1 } = {-1, 1} and { ζ : ζ 3 = 1 } = {ζ 2 6 , ζ 4
6 , ζ 6 6 = 1}. According to the disjoint decomposition P = P ≤N -k ∪ P >N -k , we write

E P >N -k (t) = - ζ∈P >N -k Res z=ζ z -t-1 F (α)(z) dz and E P ≤N -k (t) = - ζ∈P ≤N -k Res z=ζ z -t-1 F (α)(z) dz.
The following proposition is a direct consequence of Theorem 2.1.

Proposition 2.3. We have

E(α)(t) = E P >N -k (t) + E P ≤N -k (t),
where the function E P ≤N -k (t) is a quasi-polynomial function in the variable t of degree strictly less than N -k.

Thus for the purpose of computing Top k E(α)(t) it is sufficient to compute the function E P >N -k (t). This function is computable in polynomial time, as stated in the main result of our paper: Theorem 2.4. Let k be a fixed number. Then the coefficient functions of the quasi-polynomial function E P >N -k (t) are computable in polynomial time as step polynomials of t.

We prove the theorem in the rest of this section and the next.

2.3.

The poset of the high-order poles. We first rewrite our set P >N -k . Note that if ζ is a pole of order ≥ p, this means that there exist at least p elements α i in the list α so that

ζ α i = 1. But if ζ α i = 1 for a set I ⊆ {1, . . . , N + 1} of indices i, this is equivalent to the fact that ζ f = 1, for f the greatest common divisor of the elements α i , i ∈ I.
Now let I >N -k be the set of subsets of {1, . . . , N + 1} of cardinality greater than N -k. Note that when k is fixed, the cardinality of I >N -k is a polynomial function of N . For each subset I ∈ I >N -k , define f I to be the greatest common divisor of the corresponding sublist α i , i ∈ I. Let G >N -k (α) = { f I : I ∈ I >N -k } be the set of integers so obtained and let G(f ) ⊂ C × be the group of f -th roots of unity,

G(f ) = { ζ ∈ C : ζ f = 1 }. The set { G(f ) : f ∈ G >N -k (α) } forms a poset P>N-k (partially ordered set) with respect to reverse inclusion. That is, G(f i ) P>N-k G(f j ) if G(f j ) ⊆ G(f i ) (the i and j become swapped). Notice G(f j ) ⊆ G(f i ) ⇔ f j divides f i . Even if P>N-k
has a unique minimal element, we add an element 0 such that 0 G(f ) and call this new poset P >N -k .

In terms of the group G(f ) we have thus

P >N -k = f ∈G >N -k (α) G(f )
. This is, of course, not a disjoint union, but using the inclusion-exclusion principle, we can write the indicator function of the set P >N -k as a linear combination of indicator functions of the sets G(f ):

[P >N -k ] = f ∈G >N -k (α) µ >N -k (f )[G(f )],
where µ >N -k (f ) := -µ >N -k ( 0, G(f )) and µ >N -k (x, y) is the standard Möbius function for the poset P >N -k :

µ >N -k (s, s) = 1 ∀s ∈ P >N -k , µ >N -k (s, u) = - s t≺u µ >N -k (s, t) ∀s ≺ u in P >N -k .
For simplicity, µ >N -k will be called the Möbius function for the poset P >N -k and will be denoted simply by µ(f ). We also have the relationship

µ(f ) = -µ >N -k ( 0, G(f )) = 1 + 0≺G(t)≺G(f ) µ >N -k ( 0, G(t)) = 1 - 0≺G(t)≺G(f ) -µ >N -k ( 0, G(t)) = 1 - 0≺G(t)≺G(f ) µ(t).
G >N -k (α) = {1, 1, 2, 1, 1} = {1, 2}. Accordingly, P >N -k = G(1) ∪ G(2). The poset P >2 is G(1)
G( 2)

0
The arrows denote subsets, that is G(1) ⊂ G(2) and 0 can be identified with the unit circle. The Möbius function µ is simply given by µ(1) = 0, µ(2) = 1, and so

[P >N -k ] = [G(2)]. (b) Now I >N -k = I >2 = {, , }, {, , }, . . . , {, , }, {, , , }, {, , , }, {, , , }, {, , , }, {, , , }, {, , , , } and thus G >N -k (α) = {2, 3, 1, 1} = {1, 2, 3}. Hence P >N -k = G(1) ∪ G(2) ∪ G(3) = {1} ∪ {-1, 1} ∪ {ζ 3 , ζ 2 3 , 1}, where ζ 3 = e 2πi/3 is a primitive 3rd root of unity. G(1) G(2) G(3) 0 The Möbius function µ is then µ(3) = 1, µ(2) = 1, µ(1) = -1, and thus [P >N -k ] = -[G(1)] + [G(2)] + [G(3)].
Theorem 2.6. Given a list α = [α 1 , . . . , α N +1 ] and a fixed integer k, then the values for the Möbius function for the poset P >N -k can be computed in polynomial time.

Proof. First find the greatest common divisor of all sublists of the list α with size greater than N -k. Let V be the set of integers obtained from all such greatest common divisors. We note that each node of the poset P >N -k is a group of roots of unity G(v). But it is labeled by a non-negative integer v.

Construct The number of sublists of α considered is

return 3: end if 4: L ← { v ∈ V : n | v } \ {n} 5: if L = ∅ then 6: M [n] ← 1 7: return 8: end if 9: M [n] ← 0 10: for all v ∈ L do 11: findMöbius(v, L, M ) 12: M [n] ← M [n] + M [v] 13: end for 14: M [n] ← 1 -M [n]
N 1 + N 2 + • • • + N k = O(N k ),

Let us define for any positive integer

f E(α, f )(t) = - ζ: ζ f =1 Res z=ζ z -t-1 F (α)(z) dz.
Proposition 2.7. Let k be a fixed integer, then

E P >N -k (t) = - f ∈G >N -k (α) µ(f )E(α, f )(t).
(2.2)

Thus we have reduced the computation to the fast computation of E(α, f )(t).

3. Polyhedral reinterpretation of the generating function

E(α, f )(t)
To complete the proof of Theorem 2.4 we need only to prove the following proposition.

Proposition 3.1. For any integer f ∈ G >N -k (α), the coefficient functions of the quasipolynomial function E(α, f )(t) and hence E P >N -k (t) are computed in polynomial time as step polynomials of t.

By Proposition 2.7 we know we need to compute the value of E(α, f )(t). Our goal now is to demonstrate that this function can be thought of as the generating function of the lattice points inside a convex cone. This is a key point to guarantee good computational bounds. Before we can do that we review some preliminaries on generating functions of cones. We recall the notion of generating functions of cones; see also [START_REF] Baldoni | Computation of the highest coefficients of weighted Ehrhart quasi-polynomials of rational polyhedra[END_REF].

Let V = R r provided with a lattice Λ, and let V * denote the dual space. A (rational)

simplicial cone c = R ≥0 w 1 + • • • + R ≥0 w r
is a cone generated by r linearly independent vectors w 1 , . . . , w r of Λ. We consider the semi-rational affine cone s + c, s ∈ V . Let ξ ∈ V * be a dual vector such that ξ, w i < 0, 1 ≤ i ≤ r. Then the sum

S(s + c, Λ)(ξ) = n∈(s+c)∩Λ e ξ,n
is summable and defines an analytic function of ξ. It is well known that this function extends to a meromorphic function of ξ ∈ V * C . We still denote this meromorphic extension by S(s + c, Λ)(ξ). There is a very special and important case when the function M (s, c, Λ)(ξ) = e -ξ,s S(s + c, Λ)(ξ) is easy to write down. A unimodular cone, is a cone u whose primitive generators g u i form a basis of the lattice Λ. We introduce the following notation. Definition 3.5. Let u be a unimodular cone with primitive generators g u i and let s ∈ V . Then, write s = i s i g u i , with s i ∈ R, and define

{-s} u = i {-s i }g u i . Thus s + {-s} u = i s i g u i . Note that if t ∈ Λ, then {-(s + t)} u = {-s} u .
Thus, s → {-s} u is a function on V /Λ with value in V . For any ξ ∈ V * , we then find

S(s + u, Λ)(ξ) = e ξ,s e ξ,{-s}u 1 j (1 -e ξ,g u j )
and thus

M (s, u, Λ)(ξ) = e ξ,{-s}u 1 j (1 -e ξ,g u j )
.

For a general cone c, we can decompose its indicator function [c] as a signed sum of indicator functions of unimodular cones, u u [u], modulo indicator functions of cones containing lines. As shown by Barvinok (see [START_REF] Barvinok | Polynomial time algorithm for counting integral points in polyhedra when the dimension is fixed[END_REF] for the original source and [START_REF]Integer points in polyhedra[END_REF] for a great new exposition), if the dimension r of V is fixed, this decomposition can be computed in polynomial time. Then we can write

S(s + c, Λ)(ξ) = u u S(s + u, Λ)(ξ).
Thus we obtain, using Formula (3.3),

M (s, c, Λ)(ξ) = u u e ξ,{-s}u 1 j (1 -e ξ,g u j )
.

(3.4)

Here u runs through all the unimodular cones occurring in the decomposition of c, and the g u j ∈ Λ are the corresponding generators of the unimodular cone u. Remark 3.6. For computing explicit examples, it is convenient to make a change of variables that leads to computations in the standard lattice Z r . Let B be the matrix whose columns are the generators of the lattice Λ; then Λ = BZ r .

M (s, c, Λ)(ξ) = e -ξ,s n∈(s+c)∩BZ r e ξ,n = e -B ξ,B -1 s x∈(B -1 (s+c)∩Z r e B ξ,x = M (B -1 s, B -1 c, Z r )(B ξ).
3.1. Back to the computation of E(α, f )(t). After the preliminaries we will see how to rewrite E(α, f )(t) in terms of lattice points of simplicial cones. This will require some suitable manipulation of the initial form of E(α, f )(t). To start with, define the function

E(α, f )(t, T ) = -res x=0 e -tx ζ: ζ f =1 ζ -T N +1 i=1 (1 -ζ α i e α i x )
.

Writing z = ζe x , changing coordinates in residue and computing dz = z dx we write:

E(α, f )(t, T ) = -res z=ζ z -t-1 ζ t ζ : ζ f =1 ζ -T N +1 i=1 (1 -z α i )
.

By evaluating at T = t, we obtain:

E(α, f )(t) = E(α, f )(t, T ) T =t . (3.5)
We can now define:

Definition 3.7. Let k be fixed. For f ∈ G >N -k (α), define F(α, f, T )(x) := ζ: ζ f =1 ζ -T N +1 i=1 (1 -ζ α i e α i x )
, and

E i (f )(T ) := res x=0 (-x) i i! F(α, f, T )(x). Then E(α, f )(t, T ) = -res x=0 e -tx F(α, f, T )(x).
The dependence in T of F(α, f, T )(x) is through ζ T . As ζ f = 1, the function F(α, f, T )(x) is a periodic function of T modulo f whose values are meromorphic functions of x. Since the pole in x is of order at most N + 1, we can rewrite E(α, f )(t, T ) in terms of E i (f )(T ) and prove: Theorem 3.8. Let k be fixed. Then for f ∈ G >N -k (α) we can write

E(α, f )(t, T ) = N i=0 t i E i (f )(T )
with E i (f )(T ) a step polynomial of degree less than or equal to N -i and periodic of T modulo f . This step polynomial can be computed in polynomial time.

It is now clear that once we have proved Theorem 3.8, then the proof of Theorem 2.4 will follow. Writing everything out, for m such that 0 ≤ m ≤ N , the coefficient of t m in the Ehrhart quasi-polynomial is given by

E m (T ) = -res x=0 (-x) m m! f ∈G>m(α) µ(f ) ζ: ζ f =1 ζ -T i (1 -ζ α i e α i x ) . (3.6)
As an example, we see that E N is indeed independent of T because G >N (α) = {1}; thus E N is a constant. We now concentrate on writing the function F(α, f, T )(x) more explicitly.

Definition 3.9. For a list α and integers f and T , define meromorphic functions of x ∈ C by:

B(α, f )(x) := 1 i : f |α i (1 -e α i x ) , S(α, f, T )(x) := ζ: ζ f =1 ζ -T i : f α i (1 -ζ α i e α i x )
.

Thus F(α, f, T )(x) = B(α, f )(x) S(α, f, T )(x)
. The expression we obtained will allow us to compute F(α, f, T ) by relating S(α, f, T ) to a generating function of a cone. This cone will have fixed dimension when k is fixed.

3.2. E(α, f )(t) as the generating function of a cone in fixed dimension. To this end, let f be an integer from G >N -k (α). By definition, f is the greatest common divisor of a sublist of α. Thus the greatest common divisor of f and the elements of α which are not a multiple of f is still equal to 1. Let J = J(α, f ) be the set of indices i ∈ {1, . . . , N + 1} such that α i is indivisible by f , i.e., f α i . Note that f by definition is the greatest common divisor of all except at most k of the integers α j . Let r denote the cardinality of J; then r ≤ k. Let V J = R J and let V * J denote the dual space. We will use the standard basis of R J , and we denote by R J ≥0 the standard cone of elements in R J having non-negative coordinates. We also define the sublist α J = [α i ] i∈J of elements of α indivisible by f and view it as a vector in V * J via the standard basis. Definition 3.10. For an integer T , define the meromorphic function of ξ ∈ V * J ,

Q(α, f, T )(ξ) := ζ: ζ f =1 ζ -T j∈J(α,f ) (1 -ζ α j e ξ j )
.

Remark 3.11. Observe that Q(α, f, T ) can be restricted at ξ = α J x, for x ∈ C generic, to give S(α, f, T )(x).

We find that Q(α, f, T )(ξ) is the discrete generating function of an affine shift of the standard cone R J ≥0 relative to a certain lattice in V J which we define as:

Λ(α, f ) := y ∈ Z J : α J , y = j∈J y j α j ∈ Zf . (3.7)
Consider the map φ : Z J → Z/Zf , y → α, y +Zf . Its kernel is the lattice Λ(α, f ). Because the greatest common divisor of f and the elements of α J is 1, by Bezout's theorem there exist s 0 ∈ Z and s ∈ Z J such that 1 = i∈J s i α i + s 0 f . Therefore, the map φ is surjective, and therefore the index |Z J : Λ(α, f )| equals f . Theorem 3.12. Let α = [α 1 , . . . , α N +1 ] be a list of positive integers and f be the greatest common divisor of a sublist of α. Let J = J(α, f ) = { i : f α i }. Let s 0 ∈ Z and s ∈ Z J such that 1 = i∈J s i α i + s 0 f using Bezout's theorem. Consider s = (s i ) i∈J as an element of V J = R J . Let T be an integer, and

ξ = (ξ i ) i∈J ∈ V * J with ξ i < 0. Then Q(α, f, T )(ξ) = f e ξ,T s n∈(-T s+R J ≥0 )∩Λ(α,f ) e ξ,n
Remark 3.13. The function Q(α, f, T )(ξ) is a function of T periodic modulo f . Since f Z J is contained in Λ(α, f ), the element f s is in the lattice Λ(α, f ), and we see that the right hand side is also a periodic function of T modulo f .

Proof of Theorem 3.12. Consider ξ ∈ V * J with ξ j < 0. Then we can write the equality 1

j∈J (1 -ζ α j e ξ j ) = j∈J ∞ n j =0 ζ n j α j e n j ξ j . So Q(α, f, T )(ξ) = n∈Z J ≥0 ζ : ζ f =1
ζ j n j α j -T e j∈J n j ξ j .

We note that ζ:

ζ f =1 ζ m is zero except if m ∈ Zf
, when this sum is equal to f . Then we obtain that Q(α, f, T ) is the sum over n ∈ Z J ≥0 such that j n j α j -T ∈ Zf . The equality 1 = j∈J s j α j + s 0 f implies that T ≡ j ts j α j modulo f , and the condition j n j α j -T ∈ Zf is equivalent to the condition j (n j -T s j )α j ∈ Zf . We see that the point n -T s is in the lattice Λ(α, f ) as well as in the cone -T s + R J ≥0 (as n j ≥ 0). Thus the claim.

By definition of the meromorphic functions

S -T s+R J ≥0 , Λ(α, f ) (ξ) and M -T s, R J ≥0 , Λ(α, f ) (ξ)
, we obtain the following equality.

F(α, f, T )(x) = f M -T s, R J ≥0 , Λ(α, f ) (α J x) j : f |α j 1 1 -e α j x .
3.3. Unimodular decomposition in the dual space. The cone R J ≥0 is in general not unimodular with respect to the lattice Λ(α, f ). By decomposing R J ≥0 in cones u that are unimodular with respect to Λ(α, f ), modulo cones containing lines, we can write

M -T s, R J ≥0 , Λ(α, f ) = u u M (-T s, u, Λ),
where u ∈ {±1}. This decomposition can be computed using Barvinok's algorithm in polynomial time for fixed k because the dimension |J| is at most k.

Remark 3.16. For this particular cone and lattice, this decomposition modulo cones containing lines is best done using the "dual" variant of Barvinok's algorithm, as introduced in [START_REF] Barvinok | An algorithmic theory of lattice points in polyhedra[END_REF]. This is in contrast to the "primal" variant described in [START_REF] Brion | Residue formulae, vector partition functions and lattice points in rational polytopes[END_REF][START_REF] Köppe | Computing parametric rational generating functions with a primal Barvinok algorithm[END_REF]; see also [START_REF] Baldoni | Intermediate sums on polyhedra: Computation and real Ehrhart theory[END_REF] for an exposition of Brion-Vergne decomposition and its relation to both decompositions. To explain this, let us determine the index of the cone R J ≥0 in the lattice Λ = Λ(α, f ); the worst-case complexity of the signed cone decomposition is bounded by a polynomial in the logarithm of this index.

Let B be a matrix whose columns form a basis of Λ, so Λ = BZ J . Then |Z J : Λ| = |det B| = f . By Remark 3.6, we find

M -T s, R J ≥0 , Λ (ξ) = M (-T B -1 s, B -1 R J ≥0 , Z J )(B ξ). Let c denote the cone B -1 R J ≥0
, which is generated by the columns of B -1 . Since B -1 is not integer in general, we find generators of c that are primitive vectors of Z J by scaling each of the columns by an integer. Certainly |det B|B -1 is an integer matrix, and thus we find that the index of the cone c is bounded above by f r-1 . We can easily determine the exact index as follows. For each i ∈ J, the generator e i of the original cone R J ≥0 needs to be scaled so as to lie in the lattice Λ. The smallest multiplier y i ∈ Z >0 such that α J , y i e i ∈ Zf is

y i = lcm(α i , f )/α i . Thus the index of R J
≥0 in Z J is the product of the y i , and finally the

index of R J ≥0 in Λ is 1 |Z r : Λ| i∈J lcm(α i , f ) α i = 1 f i∈J lcm(α i , f ) α i .
Instead we consider the dual cone, c

• = { η ∈ V * J : η, y ≥ 0 for y ∈ c }. We have c • = B R J ≥0 . Then the index of the dual cone c • equals |det B | = f , which is much smaller than f r-1 .
Following [START_REF] Dyer | On Barvinok's algorithm for counting lattice points in fixed dimension[END_REF], we now compute a decomposition of c • in cones u • that are unimodular with respect to Z J , modulo lower-dimensional cones,

[c • ] ≡ u u [u • ] (modulo lower-dimensional cones).
Then the desired decomposition follows:

[c] ≡ u u [u] (modulo cones with lines).

Because of the better bound on the index of the cone on the dual side, the worst-case complexity of the signed decomposition algorithm is reduced. This is confirmed by computational experiments.

Remark 3.17. Although we know that the meromorphic function M -T s, R J ≥0 , Λ(α, f ) (ξ) restricts via ξ = α J x to a meromorphic function of a single variable x, it may happen that the individual functions M -T s, u, Λ(α, f ) (ξ) do not restrict. In other words, the line α J x may be entirely contained in the set of poles. If this is the case, we can compute (in polynomial time) a regular vector β ∈ Q J so that, for = 0, the deformed vector (α J + β)x is not a pole of any of the functions M -T s, u, Λ(α, f ) (ξ) occurring. We then consider the meromorphic functions → M -T s, u, Λ(α, f ) ((α J + β)x) and their Laurent expansions at = 0 in the variable . We then add the constant terms of these expansions (multiplied by u ). This is the value of M -T s, R J ≥0 , Λ(α, f ) (ξ) at the point ξ = α J x. 3.4. The periodic dependence in T . Now let us analyze the dependence in T of the functions M (-T s, u, Λ(α, f )), where u is a unimodular cone. Let the generators be g u i , so the elements g u i form a basis of the lattice Λ(α, f ). Recall that the lattice f Z r is contained in Λ(α, f ). Thus as s ∈ Z r , we have s = i s i g u i with f s i ∈ Z and hence {-T s} u = i {-T s i }g u i with {-T s i } a function of T periodic modulo f .

Thus the function T → {-T s} u is a step linear function, modulo f , with value in V . We then write

M (-T s, u, Λ(α, f ))(ξ) = e ξ,{T s}u r j=1 1 1 -e ξ,g j .
Recall that by Corollary 3.15,

F(α, f, T )(x) = f M -T s, R J ≥0 , Λ(α, f ) (α J x) j : f |α j 1 1 -e α j x .
Thus this is a meromorphic function of the variable x of the form:

u e lu(T )x h(x)
x N +1 , where h(x) is holomorphic in x and l u (T ) is a step linear function of T , modulo f . Thus to compute

E i (f )(T ) = res x=0 (-x) i i! F(α, f, T )(x)
we only have to expand the function x → e lu(T )x up to the power x N -i . This expansion can be done in polynomial time. We thus see that, as stated in Theorem 3.8, E i (f )(T ) is a step polynomial of degree less than or equal to N -i, which is periodic of T modulo f . This completes the proof of Theorem 3.8 and thus the proof of Theorem 2.4.

Periodicity of coefficients

Now that we have the main algorithmic result we can prove some consequences to the description of the periodicity of the coefficients. In this section, we determine the largest i with a non-constant coefficient E i (t) and we give a polynomial time algorithm for computing it. This will complete the proof of Theorem 1.3. Theorem 4.1. Given as input a list of integers α = [α 1 , . . . , α N +1 ] with their prime factorization

α i = p a i1 1 p a i2 2 • • • p a in n ,
there is a polynomial time algorithm to find all of the largest sublists where the greatest common divisor is not one. Moreover, if denotes the size of the largest sublists with greatest common divisor different from one, then (1) there are polynomially many such sublists, (2) the poset P> -1 is a fan (a poset with a maximal element and adjacent atoms), and (3) the Möbius function for

P > -1 is µ(f ) = 1 if G(f ) = G(1) and µ(1) = 1 -(|G > -1 (α)| -1).
Proof. Consider the matrix A = [a ij ]. Let c i 1 , . . . , c i k be column indices of A that denote the columns that contain the largest number of non-zero elements among the columns. Let α (c i j ) be the sublist of α that corresponds to the rows of A where column c i j has a non-zero entry. Each α (c i j ) has greatest common divisor different from one. If is the size of the largest sublist of α with greatest common divisor different from one, then there are many α i 's that share a common prime. Hence each column c i 1 of A has many non-zero elements. Then each α (c i j ) is a largest sublist where the greatest common divisor is not one. Note that more than one column index c i might produce the same sublist α (c i j ) . The construction of A, counting the non-zero elements of each column, and forming the sublist indexed by each c i j can be done in polynomial time in the input size.

To show the poset P> -1 is a fan, let G = {1, f 1 , . . . , f m } be the set of greatest common divisors of sublists of size > -1. Each f i corresponds to a greatest common divisor of a sublist α (i) of α with size . We cannot have f i | f j for i = j because if f i | f j , then f i is also the greatest common divisor of α (i) ∪ α (j) , a contradiction to the maximality of . Then the Möbius function is µ(f i ) = 1, and µ(1) = 1 -m.

As an aside, gcd(f i , f j ) = 1 for all f i = f j as if gcd(f i , f j ) = 1, then we can take the union of the sublist that produced f i and f j thereby giving a larger sublist with greatest common divisor not equal to one, a contradiction. where the columns are the powers of the primes indexed by [START_REF] Agnarsson | On the Sylvester denumerants for general restricted partitions[END_REF][START_REF] Barvinok | Short rational generating functions for lattice point problems[END_REF][START_REF] Barvinok | An algorithmic theory of lattice points in polyhedra[END_REF][START_REF] Dyer | On Barvinok's algorithm for counting lattice points in fixed dimension[END_REF]41). We see the largest sublists that have gcd not equal to one are [2 2 7 4 41 1 , 2 1 7 2 11 1 ] and [2

1 7 2 11 1 , 11 4 ]. Then G = {1, 2 1 7 2 , 11}. The poset P >1 is G(1) G(2 1 7 2 ) G (11)
0 and µ(1) = -1, µ(11) = µ(2 1 7 2 ) = 1.
Proof of Theorem 1.3. Let be the greatest integer for which there exists a sublist α J with |J| = , such that its gcd f is not 1. Then for m ≥ the coefficient of degree m, E m (T ), is constant because in Equation (3.6), G >m (α) = {1}. Hence E m (T ) does not depend on T . We now focus on E -1 (T ). To simplify Equation (3.6), we first compute the µ(f ) values.

Lemma 4.3. For as in Theorem 1.3, the poset G > -1 (α) is a fan, with one maximal element 1 and adjacent elements f which are pairwise coprime. In particular, µ(f ) = 1 for f = 1.

Proof. Let α J 1 , α J 2 be two sublists of length with gcd's f 1 = f 2 both not equal to 1. If f 1 and f 2 had a nontrivial common divisor d, then the list α J 1 ∪J 2 would have a gcd not equal to 1, in contradiction with its length being strictly greater than .

Next we recall a fact about Fourier series and use it to show that each term in the summation over f ∈ G > -1 (α) in Equation (3.6) has smallest period equal to f . Lemma 4.4. Let f be a positive integer and let φ(t) be a periodic function on Z/f Z with Fourier expansion

φ(t) = f -1 n=0
c n e 2iπnt/f . If c n = 0 for some n which is coprime to f then φ(t) has smallest period equal to f . Proof. Assume φ(t) has period m with f = qm and q > 1. We write its Fourier series as a function of period m.

φ(t) = m-1 j=0 c j e 2iπjt/m = m-1 j=0
c j e 2iπ(jq)t/f . By uniqueness of the Fourier coefficients, we have c n = 0 if n is not a multiple of q (and c qj = c j ). In particular, c n = 0 if n is coprime to f , a contradiction. Theorem 1.3 is thus the consequence of the following lemma. Proof. For f = 1, the statement is clear. Assume f = 1. We observe that the f -term in (3.6) is a periodic function (of period f ) which is given as the sum of its Fourier expansion and is written as f -1 n=0 c n e -2iπnT /f where

c n = -res x=0 (-x) -1 ( -1)! j 1 -e -2iπnα j /f e α j x .
Consider a coefficient for which n is coprime to f . We decompose the product according to whether f divides α j or not. The crucial observation is that there are exactly indices j such that f divides α j , because of the maximality assumption on . Therefore x = 0 is a simple pole and the residue is readily computed. We obtain

c n = (-1) -1 ( -1)! • 1 j:f α j 1 -e 2iπnα j /f • 1 j:f |α j α j .
Thus c n = 0 for an n coprime with f . By Lemma 4.4, each f -term has minimal period f .

As the various numbers f in G > -1 (α) different from 1 are pairwise coprime and the corresponding terms have minimal period f , E -1 (T ) has minimal period

f ∈G > -1 (α) f > 1.
This completes the proof of Theorem 1.3.

Summary of the algorithm and computational experiments

In this last section we report on experiments using our algorithm. But first, let us review the key steps of the algorithm:

Given a sequence of integers α of length N + 1, we wish to compute the top k+1 coefficients of the quasi-polynomial E(α)(t) of degree N . Recall that

E(α)(t) = N i=0 E i (t)t i
where E i (t) is a periodic function of t modulo some period q i . We assume that greatest common divisor of the list α is 1. (2) Next, we run our algorithms on a few knapsacks that have been studied in the literature. We chose these examples because some of these problems are considered difficult in the literature. We also present a comparison with other available software that can also compute information of the denumerant E α (t): the codes CTEuclid6 [START_REF] Xin | A Euclid style algorithm for MacMahon partition analysis[END_REF] and pSn [START_REF] Sills | Formulae for the number of partitions of n into at most m parts (using the quasi-polynomial ansatz)[END_REF]. 7 These codes use mathematical ideas that are different from those used in this paper.

All computations were performed on a 64-bit Ubuntu machine with 64 GB of RAM and eight Dual Core AMD Opteron 880 processors. 5.2. M-Knapsack vs. LattE Knapsack vs. LattE Top-Ehrhart. Here we compare our two implementations with the LattE Top-Ehrhart algorithm from [START_REF] Baldoni | Computation of the highest coefficients of weighted Ehrhart quasi-polynomials of rational polyhedra[END_REF]. We constructed a test set of 768 knapsacks. For each 3 ≤ d ≤ 50, we constructed four families of knapsacks: random-3: Five random knapsacks in dimension d -1 where a 1 = 1 and the other coefficients, a 2 , . . . , a d , are 3-digit random numbers picked uniformly random-15: Similar to the previous case, but with a 15-digit random number repeat: Five knapsacks in dimension d -1 where α 1 = 1 and all the other α i 's are the same 3-digit random number. These produce few poles and have a simple poset structure. These are among the simplest knapsacks that produce periodic coefficients. partition: One knapsack in the form α i = i for 1 ≤ i ≤ d.

For each knapsack, we successively compute the highest degree terms of the quasi-polynomial, with a time limit of 200 CPU seconds for each coefficient. Once a term takes longer than 200 seconds to compute, we skip the remaining terms, as they are harder to compute than the previous ones. We then count the maximum number of terms of the quasi-polynomial, starting from the highest degree term (which would, of course, be trivial to compute), that can be computed subject to these time limits. Figures 1,2, 3, 4 show these maximum numbers of terms for the random-3, random-15, repeat, and partition knapsacks, respectively. For example, in Figure 1, for each of the five random 3-digit knapsacks in ambient dimension 50, the LattE Knapsack method computed at most 6 terms of an Ehrhart polynomial, the M-Knapsack computed at most four terms, and the LattE Top-Ehrhart method computed at most the trivially computable highest degree term.

In each knapsack family, we see that each algorithm has a "peak" dimension where after it, the number of terms that can be computed subject to the time limit quickly decreases; for the LattE Knapsack method, this is around dimension 25 in each knapsack family. In each family, there is a clear order to which algorithm can compute the most: LattE Knapsack computes the most coefficients, while the LattE Top-Ehrhart method computes the least number of terms. In Figure 3, the simple poset structure helps every method to compute more terms, but the two Maple scripts seem to benefit more than the LattE Knapsack method.

Figure 4 demonstrates the power of the LattE implementation. Note that a knapsack of this particular form in dimension d does not start to have periodic terms until around d/2. Thus even though half of the coefficients are only constants we see that the M-Knapsack code cannot compute past a few periodic term in dimension 10-15 while the LattE Knapsack method is able to compute the entire polynomial. In Figure 5 we plot the average speedup ratio between the M-Knapsack and LattE Top-Ehrhart implementations along with the maximum and minimum speedup ratios (we wrote both algorithms in Maple). The ratios are given by the time it takes LattE Top-Ehrhart to compute a term, divided by the time it takes M-Knapsack to compute the same term, where both times are between 0 and 200 seconds. For example, among all the terms computed in dimension 15 from random 15-digit knapsacks, the average speedup between the two methods was 8000, the maximum ratio was 20000, and the minimum ratio was 200. We see that in dimensions 3-10, there are a few terms for which the LattE Top-Ehrhart method was faster than the M-Knapsack method, but this only occurs for the highest degree terms. Also, after dimension 25, there is little variance in the ratios because the LattE Top-Ehrhart method is only computing the trivial highest term. Similar results hold for the other knapsack families, and so their plots are omitted.

Other examples.

Next we focus on ten problems listed in Table 1. Some of these selected problems have been studied before in the literature [START_REF] Vergne | Hard equality constrained integer knapsacks[END_REF][START_REF] De Loera | Effective lattice point counting in rational convex polytopes[END_REF][START_REF] Xin | A Euclid style algorithm for MacMahon partition analysis[END_REF][START_REF] Xin | A fast algorithm for MacMahon's partition analysis[END_REF]. Table 2 shows the time in seconds to compute the entire denumerant using the M-Knapsack , LattE Knapsack and LattE Top-Ehrhart codes with two other algorithms: CTEuclid6 and pSn.

The CTEuclid6 algorithm [START_REF] Xin | A Euclid style algorithm for MacMahon partition analysis[END_REF] computes the lattice point count of a polytope, and supersedes an earlier algorithm in [START_REF] Xin | A fast algorithm for MacMahon's partition analysis[END_REF]. 8 Instead of using Barvinok's algorithm to construct unimodular cones, the main idea used by the CTEuclid6 algorithm to find the constant term in the generating function F (α)(z) relies on recursively computing partial fraction decompositions to construct the series. Notice that the CTEuclid6 method only computes the number of integer points in one dilation of a polytope and not the full Ehrhart polynomial. We can estimate how long it would take to find the Ehrhart polynomial using an interpolation method by computing the time it takes to find one lattice point count times the periodicity of the polynomial and degree. Hence, in Table 2, column "one point" refers to the running time of finding one lattice point count, while column "estimate" is an estimate for how long it would take to find the Ehrhart polynomial by interpolation. We see that the CTEuclid6 algorithm is fast for finding the number of integer points in a knapsack, but this would lead to a slow method for finding the Ehrhart polynomial.

The pSn algorithm of [START_REF] Sills | Formulae for the number of partitions of n into at most m parts (using the quasi-polynomial ansatz)[END_REF] computes the entire denumerant by using a partial fraction decomposition based method. 9 More precisely the quasi-polynomials are represented as a function f (t) given by q polynomials f [1] (t), f [2] (t), . . . , f [q] (t) such that f (t) = f [i] (t) when t ≡ i (mod q). To find the coefficients of the f [i] their method finds the first few terms of the Maclaurin expansion of the partial fraction decomposition to find enough evaluations of those polynomials and then recovers the coefficients of each the f [i] as a result of solving a linear system. This algorithm goes back to Cayley and it was implemented in Maple. Looking at Table 2, we see that the pSn method is competitive with LattE Knapsack for knapsacks 1, 2, . . . , 6, and beats LattE Knapsack in knapsack 10. However, the pSn method is highly sensitive to the number of digits in the knapsack coefficients, unlike our M-Knapsack and LattE Knapsack methods. For example, the knapsacks [START_REF] Vergne | Hard equality constrained integer knapsacks[END_REF][START_REF] Agnarsson | On the Sylvester denumerants for general restricted partitions[END_REF][START_REF] Baldoni | A user's guide for LattE integrale v1.7.2[END_REF][START_REF] Baldoni | Intermediate sums on polyhedra: Computation and real Ehrhart theory[END_REF][START_REF] Barvinok | Polynomial time algorithm for counting integral points in polyhedra when the dimension is fixed[END_REF] takes 0.320 seconds to find the full Ehrhart polynomial, [START_REF] Vergne | Hard equality constrained integer knapsacks[END_REF][START_REF] Hemmecke | Computing holes in semi-groups and its application to transportation problems[END_REF]40,60, 80] takes 5.520 seconds, and [START_REF] Vergne | Hard equality constrained integer knapsacks[END_REF]200,600,900,400] takes 247.939 seconds. Similar results hold for other three-digit knapsacks in dimension four. However, the partition knapsack [1, 2, 3, . . . , 50] only takes 102.7 seconds. Finally, comparing the two Maple scripts, the LattE Top-Ehrhart method outperforms the M-Knapsack method.

Table 2 ignores one of the main features of our algorithm: that it can compute just the top k terms of the Ehrhart polynomial. In Table 3, we time the computation for finding the top three and four terms of the Ehrhart polynomial on the knapsacks in Table 1. We immediately see that our LattE Knapsack method takes less than one thousandth of a second in each example. Comparing the two Maple scripts, M-Knapsack greatly outperforms LattE Top-Ehrhart. Hence, for a fixed k, the LattE Knapsack is the fastest method.

In summary, the LattE Knapsack is the fastest method for computing the top k terms of the Ehrhart polynomial. The LattE Knapsack method can also compute the full Ehrhart polynomial in a reasonable amount of time up to around dimension 25, and the number of digits in each knapsack coefficient does not significantly alter performance. However, if the coefficients each have one or two digits, the pSn method is faster, even in large dimensions.

  the set of poles of the meromorphic function F (α) and by p(ζ) the order of the pole ζ for ζ ∈ P.

Example 2 . 2 .

 22 (a) Let α = [98, 59, 44, 100], so N = 3, and let k = 1. Then P >N -k consists of poles of order greater than 2. Of course ζ = 1 is a pole of order 4. Note that ζ = -1 is a pole of order 3. So P >N -k = { ζ : ζ 2 = 1 }. (b) Let α = [6, 2, 2, 3, 3], so N = 4, and let k = 2. Let ζ 6 = e 2πi/6 be a primitive 6th root of unity. Then ζ 6 6 = 1 is a pole of order 5, ζ 6 and ζ 5 6 are poles of order 1, and

Example 2 . 5 (

 25 Example 2.2, continued). (a) Here we have I >N -k = I >2 = {, , }, {, , }, {, , }, {, , }, {, , , } and

  an array M of size |V | to keep the value of the Möbius function. Initialize M to hold the Möbius values of infinity: M [v] ← ∞ for all v ∈ V . Then call Algorithm 2.3 below with findMöbius(1, V, M ). Algorithm 1 findMöbius(n, V , M ) Input: n: the label of node G(n) in the poset P>N-k Input: V : list of numbers in the poset P>N-k Input: M : array of current Möbius values computed for P >N -k Output: updates the array M of Möbius values 1: if M [n] < ∞ then 2:

Algorithm 2 .

 2 3 terminates because the number of nodes v with M [v] = ∞ decreases to zero in each iteration. To show correctness, consider a node v in the poset P N -k . If v covers 0, then we must have M [v] = 1 as there is no other G(w) with G(f ) ⊂ G(w). Else if v does not cover 0, we set M [v] to be 1 minus the sum w: v|w M [w] which guarantees that the poles in G(v) are only counted once because w: v|w M [w] is how many times G(v) is a subset of another element that has already been counted.

  which is a polynomial for k fixed. For each sublist, the greatest common divisor of a set of integers is computed in polynomial time. Hence |V | = O(N k ). Notice that lines 4 to 14 of Algorithm 2.3 are executed at most O(|V |) times as once a M [v] value is computed, it is never recomputed. The number of additions on line 12 is O(|V | 2 ) while the number of divisions on line 4 is also O(|V | 2 ). Hence this algorithm finds the Möbius function in O(|V | 2 ) = O(N 2k ) time where k is fixed.

Example 3 . 2 .Theorem 3 . 3 .Definition 3 . 4 .

 323334 Let V = R with lattice Z, c = R ≥0 , and s ∈ R. Then S(s + R ≥0 , Z)(ξ) = n≥s e nξ = e s ξ 1 1 -e ξ . Using the function {x} = x -x , we find s = s + {-s} and can write e -sξ S(s + R ≥0 , Z)(ξ) = e {-s}ξ 1 -e ξ . (3.1) Recall the following result: Consider the semi-rational affine cone s+c and the lattice Λ. The series S(s+ c, Λ)(ξ) is a meromorphic function of ξ such that r i=1 ξ, w i • S(s + c, Λ)(ξ) is holomorphic in a neighborhood of 0. Let t ∈ Λ. Consider the translated cone t+s+c of s+c by t. Then we have the covariance formula S(t + s + c, Λ)(ξ) = e ξ,t S(s + c, Λ)(ξ). (3.2) Because of this formula, it is convenient to introduce the following function. Define the function M (s, c, Λ)(ξ) := e -ξ,s S(s + c, Λ)(ξ). Thus the function s → M (s, c, Λ)(ξ) is a function of s ∈ V /Λ (a periodic function of s) whose values are meromorphic functions of ξ. It is interesting to introduce this modified function since, as seen in Equation (3.1) in Example 3.2, its dependance in s is via step linear functions of s.

Corollary 3. 14 .

 14 Q(α, f, T )(ξ) = f M -T s, R J ≥0 , Λ(α, f ) (ξ).Using Remark 3.11 we thus obtain by restriction to ξ = α J x the following equality. Corollary 3.15.

Example 4 . 2 .

 42 [2 2 7 4 41 1 , 2 1 7 2 11 1 , 11 4 , 17 3] gives the matrix 

Lemma 4 . 5 .

 45 Let f ∈ G > -1 (α).The term in the summation over f in (3.6) has smallest period f as a function of T .

Figure 1 .

 1 Figure 1. Random 3-digit knapsacks: Maximum number of coefficients each algorithm can compute where each coefficient takes less than 200 seconds.

Figure 2 .

 2 Figure 2. Random 15-digit knapsacks: Maximum number of coefficients each algorithm can compute where each coefficient takes less than 200 seconds.

Figure 3 .

 3 Figure 3. Repeat knapsacks: Maximum number of coefficients each algorithm can compute where each coefficient takes less than 200 seconds.

Figure 4 .

 4 Figure 4. Partition knapsacks: Maximum number of coefficients each algorithm can compute where each coefficient takes less than 200 seconds.

Figure 5 .

 5 Figure 5. Average speedup ratio (dots) between the M-Knapsack and LattE Top-Ehrhart codes along with maximum and minimum speedup ratio bounds (vertical lines) for the random 15-digit knapsacks.

  1. We haveE(α)(t) = E P >N -k (t) + E P ≤N -k (t)with E P ≤N -k (t) a periodic polynomial of degree strictly less than N -k. Computing the first k + 1 coefficients means to compute E P >N -k (t).Command line usage:dest/bin/integrate --valuation=top-ehrhart --top-ehrhart-save= output file --num-coefficients= k value LattE style knapsack file .

	Command line usage:	dest/bin/top-ehrhart-knapsack -f knapsack file -o output file -k
	k value .

6

Table 1 .

 1 Ten selected instances Problem Data

Table 2 .

 2 Computation times in seconds for finding the full Ehrhart polynomial using five different methods.

	CTEuclid6

Table 3 .

 3 Computation times in seconds for finding the top three and four terms of the Ehrhart polynomial

			Top 3 coefficients		Top 4 coefficients
		LattE	M-Knapsack LattE	LattE	M-Knapsack LattE
		Knapsack		Top-Ehrhart Knapsack	Top-Ehrhart
	#1	0	0.305	0.128	-	-	-
	#2	0	0.004	0.768	0	0.096	1.356
	#3	0	0.004	0.788	0	0.080	1.308
	#4	0	0.003	0.792	0	0.124	1.368
	#5	0	0.004	0.784	0	0.176	1.424
	#6	0	0.004	1.660	0	0.088	2.976
	#7	0	0.004	0.836	0	0.272	1.652
	#8	0	0.068	1.828	0	0.112	3.544
	#9	0	0.004	18.437	0	0.016	59.527
	#10	0	0.012	142.104	0	0.044	822.187

This notion of degree only induces a filtration, not a grading, on the algebra of step polynomials, because there exist polynomial relations between step linear functions and therefore several step-polynomial formulas with different degrees may represent the same function.

Available under the GNU General Public License at https://www.math.ucdavis.edu/ ~latte/.

Available under the GNU General Public License at https://www.math.ucdavis.edu/ ~latte/. The Maple code M-Knapsack is also available separately at https://www.math.ucdavis.edu/ ~latte/ software/packages/maple/.

Maple usage: coeff Nminusk knapsack( knapsack list , t, k value ).

Both codes can be downloaded from the locations indicated in the respective papers. Maple scripts that correspond to our tests of these codes are available at https://www.math.ucdavis.edu/ ~latte/software/ denumerantSupplemental/.

Maple usage: CTEuclid(F (α)(x)/x b , t, [x]); where b = α 1 + • • • + α N +1 .

Maple usage: QPStoTrunc(pSn( knapsack list ,n,j),n); where j is the smallest value in {100, 200, 500,

1000, 2000, 3000} that produces an answer.

Acknowledgments. We are grateful to Doron Zeilberger and an anonymous referee for suggestions and comments. The work for this article was done in large part during a SQuaRE program at the American Institute of Mathematics, Palo Alto, in March 2012. V. Baldoni was partially supported by the Cofin 40%, MIUR. De Loera was partially supported by NSF grant DMS-0914107, M. Köppe was partially supported by NSF grant DMS-0914873. B. Dutra was supported by the NSF-VIGRE grant DMS-0636297. The support received is gratefully acknowledged.

By writing [P >N

we have

3. Fix f an integer. Write

4. We fix f and let r be the number of elements α i such that α i is not a multiple of f .

We then list such α i in the list α = [α 1 , α 2 , . . . , α r ].

We introduce a lattice Λ := Λ(α, f ) ⊂ Z r and an element s ∈ Z r so that f s ∈ Λ. We decompose the standard cone R r ≥0 as a signed decomposition, modulo cones containing lines, in unimodular cones u for the lattice Λ, obtaining

.

5.

To compute E N -i (f )(T ), we compute the Laurent series of F(α, f, T )(x) at x = 0 and take the coefficient in x -N -1+i of this Laurent series. As the Laurent series of F(α, f, T )(x) starts by x -N -1 , if i is less than k, we just have to compute at most k terms of this Laurent series.

5.1.

Experiments. We first wrote a preliminary implementation of our algorithm in Maple, which we call M-Knapsack in the following. Later we developed a faster implementation in C++, which is referred to as LattE Knapsack in the following (we use the term knapsack to refer to the Diophantine problem

Both implementations are released as part of the software package LattE integrale [START_REF] Baldoni | A user's guide for LattE integrale v1.7.2[END_REF], version 1.7.2. 3 We report on two different benchmarks tests:

(1) We test the performance of the implementations M-Knapsack 4 and LattE Knapsack 5 , and also the implementation of the algorithm from [START_REF] Baldoni | Computation of the highest coefficients of weighted Ehrhart quasi-polynomials of rational polyhedra[END_REF], which refer to as LattE Top-Ehrhart 6 , on a collection of over 750 knapsacks. The latter algorithm can compute the weighted Ehrhart quasi-polynomials for simplicial polytopes, and hence it is more general than the algorithm we present in this paper, but this is the only other available algorithm for computing coefficients directly. Note that the implementations of the M-Knapsack algorithm and the main computational part of the LattE Top-Ehrhart algorithm are in Maple, making comparisons between the two easier.