\section{{\tt LattE}'s implementation of Barvinok's algorithm} \label {software}

In this section, we go through the
steps of Barvinok's algorithm, showing how we implemented them in
{\tt LattE}. Barvinok's algorithm relies on two important new
ideas: the use of rational functions as efficient data
structures and the signed decompositions of cones into unimodular cones.

Let $P \subset \R^d$ be a rational convex polyhedron and let $f(P, z)$ be the multivariate generating function defined in (\ref{generating}). Let $v$ be a vertex of $P$.  
Then, the {\em supporting cone} $K(P, v)$ of $P$ at
$v$ is $K(P, v)= v+ \,\{u \in \real^d : v+\delta u \in P \ {\rm for \ all \ sufficiently \ small} \ \delta>0\}$.
Let $V(P)$ be the vertex set of $P$. One crucial component of
Barvinok's algorithm is the ability to distribute the computation on
the vertices of the polytope. This follows from the seminal theorem of
Brion \citep{brion88}:

\begin{theorem} \citep{brion88}
Let $P$ be a rational polyhedra and let $V(P)$ be the vertex set of $P$.  
Then, $$f(P, z)=\sum_{v \in V(P)}
f(K(P,v), z).$$
\end{theorem}


\begin{example}\label{example_1}
Consider the integral quadrilateral shown in Figure \ref{brionfig}. The
vertices are $V_1=(0,0)$, $V_2=(5,0)$, $V_3=(4,2)$, and $V_4=(0,2)$.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=5 cm]{examplebrion.eps} 
\caption{A quadrilateral in Example \ref{example_1}.}\label{brionfig}
\end{center} 
\end{figure}

We obtain four rational generation functions whose formulas are

$f(K_{V_1}, z)={\frac {1}{\left (1-z_1\right )\left (1-z_2\right )}}, \quad 
f(K_{V_2}, z)=\frac{({z_1}^{5}+{z_1}^{4}z_2)}{ (1-{z_1}^{-1}) (1-z_2^2z_1^{-1})}, $


$f(K_{V_3}, z)=\frac{({z_1}^{4}{z_2}+{z_1}^{4}{z_2}^{2})}{ (1-{z_1}^{-1})
(1-z_1z_2^{-2})},
 \quad  f(K_{V_4}, z)=\frac{z_2^{2}}{(1-{z_2}^{-1} )(1-z_1 )}.$

Indeed, the result of adding the rational functions is equal to the
polynomial

${z_1}^{5}+{z_1}^{4}z_2+{z_1}^{4}+{z_1}^{4}{z_2}^{2}+z_2{z_1}^{3}+{z_1}^{3}+{z_1}^{3}{z_2}^{2
}+z_2{z_1}^{2}+{z_1}^{2}+{z_1}^{2}{z_2}^{2}+z_1z_2+z_1+z_1{z_2}^{2}+{z_2}^{2}+z_2+1$. \eopf
\end{example}

In order to use Brion's theorem for counting lattice points in convex
polyhedra, we need to know how to compute the rational generating
function of convex rational pointed cones.  For polyhedral cones this
generating function is a rational function whose numerator and
denominator have a well-understood geometric meaning (see in
\citet[Chapter 4]{stanley} and in \citet[Corollary 4.6.8]{stanley2} for a clear explanation). We
already have a ``simple'' formula when the cone is a simple cone: Let
$\{u_1, u_2, \ldots, u_k\}$ be a set of linearly independent integral
vectors of $\real^d$, where $k \leq d$.  Let $K$ be a cone which is
generated by $\{u_1, u_2, \ldots, u_k\}$, in other words, $K=\{
\lambda_1 u_1 + \lambda_2 u_2 + \ldots + \lambda_k u_k,\mbox{ for some
} \lambda_i \geq 0 \mbox{ and } i=1,2, \ldots , k\}.$ Consider the
parallelepiped $S=\{ \lambda_1 u_1+ \lambda_2 u_2+ \ldots +\lambda_k
u_k, \ 0 \leq \lambda_i < 1 \mbox{, } i=1,2,\ldots, k\}.$

It is well-known \citep{stanley} that the generating function 
for the lattice points in $K$ equals

$$\sum_{\beta \in K \cap Z^d} z^{\beta} = (\sum_{\tau \in S \cap Z^d} 
z^{\tau})
\prod^k_{i=1} \frac{1}{1-z^{u_i}}. \quad (*) $$

Thus, to derive a formula for arbitrary pointed cones one could
decompose them into simple cones, via a triangulation, and then
apply the formula above and the inclusion-exclusion principle 
in \citet[Proposition 1.2]{stanley2}. Instead, Barvinok's idea is that
it is more efficient to further decompose each simple cone into
simple unimodular cones. A {\em unimodular} cone is a simple
cone with generators $\, \{u_1,\dots,u_k\} \,$ that form an integral
basis for the lattice $\, \real \{u_1,\dots,u_k\} \cap \Z^d $. Note
that in this case the numerator of the formula has a single monomial,
in other words, the parallelepiped has only one lattice point.

\subsection{Simple signed decompositions} \label{decomposition}

We now focus our attention on how the cone decomposition is done.  To
decompose a cone into simple cones the first step is to do a
triangulation ({\em triangulation} of a cone $K$ in dimension $d$ is a
collection of $d$-dimensional simple cones such that their union
is $K$, their interiors are disjoint, and any pair of them intersect
in a (possibly empty) common face).  There are efficient algorithms,
when the dimension is fixed, to carry a triangulation (see
\citet{auren,lee} for details). In {\tt LattE} we use the well-known
Delaunay triangulation which we compute via a convex hull calculation.
The idea is to ``lift'' the rays of the cone into a higher dimensional
paraboloid by adding a new coordinate which is the sum of the
squares of the other coordinates, take the lower convex hull of the
lifted points, and then ``project'' back those simple facets. We
use Fukuda's implementation in {\tt CDD} \citep{fukuda} of this
lift-and-project algorithm.  This is not the only choice of
triangulation, and definitely not the smallest one. In Section
\ref{theory} we discuss some situations when the choice of
triangulation in fact gives a better rational function.

In principle, one could at this point list the points of the
fundamental parallelepiped, for example, using a fast Hilbert bases
code such as {\tt 4ti2} \citep{Hemmecke:Hilbert} or {\tt NORMALIZ}
\citep{brunskoch}, and then use formula $(*)$ for a general simple
cone.  Theoretically this is bad because the number of lattice points
in the parallelepiped is exponentially large already for fixed
dimension. In practice, this can often be done and in some situations
 is useful.  Barvinok instead decomposes each simple cone as a
(signed) sum of simple {\em unimodular} cones. To be more formal,
for a set $A \subset \real^d$, the indicator function $[ A ]: \real^d
\rightarrow \real$ of $A$ is defined as $$[ A ] (x) = \left \{
\begin{array}{ll} 1 \mbox{ if }x \in A, \\ 0 \mbox{ if }x \not \in
A.\\
\end{array}\right .$$ 

We want to express the indicator function of a simple cone as an
integer linear combination of the indicator functions of unimodular simple
cones.  There is a nice valuation from the algebra of
indicator functions of polyhedra to the field of rational functions \citep{BarviPom}, and many of its properties can be
used in the calculation. For example, the valuation is zero when 
the polyhedron contains a line.

\begin{theorem} \citep[Theorem 3.1]{BarviPom} 
There is a valuation $f$ from the algebra of indicator functions of
rational polyhedra into the field of multivariate rational functions
such that for any polyhedron $P$, $f([P])=\sum_{\alpha \in P \cap
\Z^d} x^\alpha$.
\end{theorem}

Therefore once we have a unimodular cone decomposition, the rational
generating function of the original cone is a signed sum of
``short'' rational functions.  Next we focus on how to decompose
a simple cone into unimodular cones.


 Let $u_1, u_2, \ldots , u_d$ be linearly independent integral vectors
which generate a simple cone $K$.  We denote the \textsl{index} of
$K$ by ind$(K)$ which tells how far $K$ is from being unimodular.
That is, ind$(K) =|\det (u_1| u_2 | \ldots | u_d)|$ which is the
volume of the parallelepiped spanned by $u_1, u_2, \ldots , u_d$. It
is also equal to the number of lattice points inside the half-open
parallelepiped. $K$ is unimodular if and only if the index of $K$ is
$1$. Now we discuss how we implemented the following key result of
Barvinok:

\begin{theorem} \citep[Theorem 4.2]{BarviPom} \label{break}
Fix the dimension $d$.  Then, there exists a polynomial time algorithm
with a given rational polyhedral cone $K \subset \real^d$, which
computes unimodular cones $K_i$, $i \in I = \{1,
2, \ldots , l\}$, and numbers $\epsilon_i \in \{ -1, 1 \}$ such that
$$[K] = \sum_{i \in I} \epsilon_i [K_i].$$
\end{theorem}

Let $K$ be a rational pointed simple cone.  Consider the closed
parallelepiped $$\Gamma = \{ \alpha_1u_1+\alpha_2u_2 + \ldots +
\alpha_du_d : \mbox{ } |\alpha_j| \leq (\ind (K))^{- \frac{1}{d}}
\mbox{, } j = 1, 2, \ldots , d\}.$$ Note that this parallelepiped
$\Gamma$ is centrally symmetric and one can show that the volume of
$\Gamma$ is $2^d$. Minkowski's First Theorem \citep{schrijver}
guarantees that because $\Gamma \subset \real^d$ is a centrally
symmetric convex body with volume $\geq 2^d$, there exists a non-zero
lattice point $w$ inside of $\Gamma$. We will use $w$ to build the
decomposition.


We need to find $w$ explicitly. We take essentially the approach
suggested by \citet{Dyer}. We require a
subroutine that computes the shortest vector in a lattice. For fixed dimension this can be done in polynomial time using lattice basis reduction (this follows trivially from \citet[Corollary 6.4b][page 72]{schrijver}).  It is worth observing that when the dimension is not fixed the problem becomes NP-hard \citep{ajtai}.  
We use the basis reduction algorithm
of Lenstra, Lenstra, and Lov\'asz \citep{grolosch,schrijver} to find a short vector. Given $A$, an
integral $d \times d$ matrix whose columns generate a lattice, LLL's
algorithm outputs $A'$, a new $d \times d$ matrix, spanning the same
lattice generated by $A$. The column vectors of $A'$, $u'_1, u'_2,
..., u'_d$, are short and nearly orthogonal to each other, and each
$u'_i$ is an approximation of the shortest vector in the lattice, in
terms of Euclidean length. It is well-known \citep{schrijver} that
there exists a unique unimodular matrix $U$ such that $A U = A'$.

The method proposed in \citet{Dyer} to find $w$ is
the following: Let $A = (u_1|u_2| \ldots | u_d)$, where the $u_i$
are the rays of the simple cone we wish to decompose. Compute the
reduced basis of $A^{-1}$ using the LLL algorithm. Let $A'= (u_1'|u_2'| \ldots | u_d')$ be the
reduced basis of $A^{-1}$.  Dyer and Kannan observed that we can find
the smallest vector with respect to the $l^{\infty}$ norm by searching
over all linear integral combinations of the column vectors of $A'$
with small coefficients. We call this search the {\em enumeration step}.  
This enumeration step can be performed in polynomial time in fixed dimension.
We will briefly describe the enumeration step. 
First we introduce some notation.  
Let $u_1', \dots , u_d'$ be linearly independent integral vectors in $\Z^d$.
Let $\| \cdot \|_2$ be the $l^2$ norm and let $\| \cdot \|_\infty$ be the 
infinity norm.


We will need to recall the Gram-Schmidt process that computes a set of orthogonal 
vectors $u^*_j$, $1\leq j \leq d$,
from independent vectors $u_j'$, $1\leq j \leq d$.  
In particular we need some values from this process.  The vectors 
$u^*_j$, and real numbers 
$\mu_{j, k}\mbox{, } 1 \leq k < j \leq d$ are computed from $u_j'$
by the recursive process: $$u^*_1=u_1'$$
$$u^*_j=u_j'-\sum_{k=1}^{j-1} \mu_{j, k} u^*_k\mbox{, } \, \, 2 \leq j \leq d$$
$$\mu_{j, k}=\frac{u_j' \cdot u^*_k}{ \| u^*_k\|_2 ^2}.$$
Letting $w_i := u_i^* /\|u_i^*\|_2$, there exists real numbers $u_i(j)$ 
such that
\begin{equation}\label{eqkannan1}
u_i' = \sum_{j=1}^d u_i(j) w_j.
\end{equation}
Note that $u_i(j) = \mu_{i, j} \|u_j^*\|_2$ 
for $1 \leq k < j \leq d$ and $u_i(i) = \|u^*_i\|_2$. 
These $u_i(j)$ will be used below.
Let $L(u_1', \dots , u_d')$ be the lattice generated by 
$u_1', \dots , u_d'$.  Then 
we denote $L_j(u_1', \dots , u_d')$ be the projection of $L(u_1', \dots , u_d')$ 
orthogonal to the vector space $V_j$ spanned by $u_1', \dots , u_j'$. 

Now we are ready to describe the process of the enumeration step.
Let $\lambda$ be a shortest vector in the
lattice spanned by $A'$ with respect to the $l^{\infty}$ norm. 
Then we can write $\lambda$ as an integral linear 
combination of columns of $A'$.  Let $\lambda = \sum_{i=1}^d \alpha_i u'_i$,
where $\alpha = (\alpha_1, \dots , \alpha_d) \in \Z^d$. 
The goal is to find some finite set $T \subset \Z^d$ such that $\alpha \in T$ and
the cardinality of $T$ is polynomial size in fixed dimension. $T$ will be
contained inside a certain parallelepiped.  Then we can
search $\lambda$ by enumerating all lattice points inside $T$.


Suppose $A' = (u_1'| \dots | u_d')$ form the reduced basis obtained by LLL algorithm.
Let $m:= \min \{j: u_j(j) \geq u_1(1) \} - 1$. 
Now we will apply the inequalities 
\begin{equation}\label{eqkannan3}
\| x \|_\infty \leq \| x \|_2,
\end{equation}

\begin{equation}\label{eqkannan4}
\| x \|_2 \leq \sqrt{d} \| x \|_\infty.
\end{equation}
We are going to prove that a shortest vector of $L(u_1', \dots , u_m')$ is a shortest
vector of $L(u_1', \dots , u_d')$ with respect to the $l^\infty$ norm.  Any 
vector in $L(u_1', \dots , u_d') \backslash L(u_1', \dots , u_m')$ must have
$l^2$ norm at least $u_1(1)$.  Since $u_1(1) = \|u_1'\|_2$ it must have $l^\infty$ norm 
at least $u_1(1)$ which is at least the $l^\infty$ norm of $u_1'$ by 
(\ref{eqkannan3}).  

We will show how to construct $T$.
Let $y = \sum_{i=1}^d \alpha_i u_i'$ be a candidate for a shortest vector
with respect to the $l^\infty$ norm.

Applying the fact that $|\alpha_i u_i(i)| \leq \|y\|_2$ 
(using the same trick as on page 423, \cite{kannan}),
we have $|\alpha_i u_i(i)| / \sqrt{d} \leq  \|y\|_2 / \sqrt{d} 
\leq  \|u_1'\|_\infty $
for any candidate vector for a shortest vector with respect to the 
$l^\infty$ norm. Therefore, we have
$$|\alpha_i|u_i(i)  / \sqrt{d} \leq
 \|u_1' \|_\infty \leq u_1(1).$$  From this
$$|\alpha_i| \leq  \sqrt{d} u_1(1)/u_i(i) \mbox{ for } i = 1, \dots, m, $$
which defines a parallelepiped in the variables $\alpha_i$ such that, 
$$Q := \{ \alpha \in \R^d: -\sqrt{d} u_1(1)/u_i(i) \leq \alpha_i \leq \sqrt{d} u_1(1)/u_i(i) \mbox{ for } i = 1, \dots , d\}.$$
Finally we set $T := Q \cap \Z^d$.

Now we are going to show that $Q$ contains polynomially many
lattice points. For each $\alpha_i$, there exist at most 
$1+2 \sqrt{d}u_1(1)/u_i(i)$ candidates.  So the total number of candidates is 
$$\prod_{i=1}^m (1+2 \sqrt{d} u_1(1)/u_i(i)).$$  
With the fact that $u_1(1) \geq u_i(i)$ (by the definition of $m$), we have 
$$\prod_{i=1}^m (1+2 \sqrt{d} u(1)_1/u_i(i)) \leq 3^m d^{m/2} \prod_{i=1}^m(u_1(1)/u_i(i)).$$
We derive the following from Minkowski's theorem 
$$u_1(1)^m \leq (2m)^{m/2} \det(L(u_1', \dots, u_m')),$$
$$\det(L(u_1', \dots, u_m')) = \prod_{i=1}^m u_i(i).$$ 
Therefore, we have $\prod_{i=1}^m(u_1(1)/u_i(i)) 
\leq (2m)^{m/2}$.  This implies that 
$$\prod_{i=1}^m (1+2 \sqrt{d} u(1)_1/u_i(i))
\leq (3d)^{d},$$ which is a constant if we fix $d$. With this method, we 
can compute a shortest vector $\lambda$ with respect to the $l^\infty$ norm
in polynomial time in fixed dimension by the enumeration step.


After we compute $\lambda$ in polynomial time in fixed dimension,
we know that
there exists a unique unimodular matrix $U$ such that $A' = A^{-1}U$.
Minkowski's theorem for the  $l^{\infty}$ norm implies 
that for the non-singular matrix $A'$, there exists a non-zero integral vector
$z$ such that $\|\lambda \|_\infty =\|A'z\|_{\infty} \leq |\det(A')|^{1/d}$.
See statement 23 in page 81 
in \citet{schrijver}. We can set

$$\|\lambda \|_{\infty} \leq |\det(A')|^{1/d} = |\det(A^{-1} U)|^{1/d}
=
|\det(A^{-1}) \det(U)|^{1/d}
$$

$$ = |\det(A^{-1})|^{1/d} = |\det(A)|^{-1/d} = |\ind(K)|^{-1/d}.$$

Since $A^{-1}$ and $A'$ span the same lattice, there exists an
integral vector $w \in \real^d$ such that $\lambda = A^{-1} w$. Then,
we have $$w = A \lambda.$$ Note that $w$ is
a non-zero integral vector which is a linear integer combination of
the generators $u_i$ of the cone $K$ {\em with possibly negative
coefficients}, and with coefficients at most $|\ind(K)|^{-1/d}$.
Therefore, we have found a non-zero integral vector $w \in
\Gamma$. In {\tt LattE}, we try to avoid the enumeration step
because it is very costly. Instead, we choose $\lambda$ to be the
shortest of the columns in $A'$. This may not be the smallest vector,
but for practical purposes, it often decreases the index $| \ind(K)|$
just like for a shortest vector. Experimentally
we have observed that we rarely use the enumeration step.

In the next step of the algorithm, for $i = 1, 2, \ldots , d$, we set
$$K_i = \co\{ u_1, u_2, \ldots , u_{i-1}, w , u_{i+1}, \ldots ,
u_d \}.$$ Now, we have to show that for each $i$, $\ind(K_i)$ is
smaller than $\ind(K)$.  Let $w = \sum_{i = 1}^d \alpha_i u_i$.
Then, we have

\begin{eqnarray*}
\ind(K_i)  &=& |\det((u_1|u_2|\ldots | u_{i-1}| w |u_{i+1}|
\ldots|u_d))| \,  \\
 &=& |\alpha_i||\det((u_1|u_2|\ldots | u_{i-1}| u_i |u_{i+1}| \ldots
|u_d))|  \,  \\
&=&|\alpha_i|\ind(K) \leq (\ind(K))^{\frac{d-1}{d}}  \, . \\ 
\end{eqnarray*}


There is one more technical condition that $w$ needs to satisfy. This
is that $w$ and $u_1,\dots,u_d$ belong to an open half-space \citep[Lemma 5.2]{bar}. This is easy to achieve as either the $w$ we
found or $-w$ satisfy this condition.  We can now decompose the
original cone $K$ into cones $K_i$ for $i= 1, 2, \ldots , d$, of
smaller index, $[K]=\sum \pm [K_i]$. This sum of indicator functions
carries signs which depend on the position of $w$ with respect to the
interior or exterior of $K$. We iterate this process until $K_i$
becomes a unimodular cone for $i = 1, 2, \ldots , d$. For implementing
Barvinok's decomposition of cones, we use the package {\tt NTL} by \citet{shoup} to compute the reduced basis of a cone and
to compute with matrices and determinants. All our calculations were
done in exact long integer arithmetic using the routines integrated in {\tt NTL}.  Here is the
pseudo-code of the algorithm and an example.

\begin{Algorithm}(Barvinok's Decomposition of a Simple Cone)


\noindent
\textbf{Input:}  A simple cone $K = \co \{u_1, u_2, 
\ldots , u_d \}$ given by its generators.

\noindent
\textbf{Output:} A list of unimodular cones and numbers $\epsilon_i$ 
as in Theorem \ref{break}.


\begin{tabbing}

\quad \= \quad \= \quad \= \quad \kill
Set two queues Uni and NonUni.\\
\keyw{if} $K$ is unimodular\\
\> \keyw{then} Uni $= \mbox{Uni } \cup \{K\}$.\\  
\keyw{else} NonUni $= \mbox{NonUni } \cup \{K\}$.\\
\keyw{while} NonUni is not empty \keyw{do}\\
\> Take a cone $K \in $ NonUni and set $A = (u_1, u_2, \ldots 
, u_d)$\\
\> \> to be a matrix  whose columns are the rays of $K$.\\
\> Compute the smallest vector $\lambda$ in the lattice,\\
\> \> with respect to  $l^{\infty}$, which is spanned by  the column vectors of $A^{-1}$.\\
\> Find a non-zero integral vector $z$ such that $\lambda = A^{-1}z$.\\
\> \keyw{if} vectors $z, u_1, u_2, \ldots , u_d$ are in an open half plane\\
\> \> \keyw{then} set $z := z$.\\
\> \keyw{else} set $z := -z$.\\
\> \keyw{for} $i = 1, 2, \ldots , d$ \keyw{do}\\
\> \> set $K_i = \co \{ u_1, \ldots , u_{i-1}, z, u_{i+1}, \ldots,
u_d \}$\\
\> \>  and set $A_i = (u_1, \ldots , u_{i-1}, z, u_{i+1}, \ldots,
u_d )$.\\
\>\keyw{for} $i = 1, 2, \ldots , d$ \keyw{do}\\
\> \> \keyw{if} $\det(A_i)$ and $\det(A)$ have the same sign\\
\> \> \> \keyw{then} assign $\epsilon_{K_i} = \epsilon_K$.\\
\> \> \keyw{else} $\epsilon_{K_i} = -\epsilon_K$.  \\
\> \keyw{for} $i = 1, 2, \ldots , d$ \keyw{do}\\
\> \>  \keyw{if} $K_i$ is unimodular\\
\> \> \> \keyw{then}  Uni = Uni $\cup \{K_i\}$. \\
\> \> \keyw{else} NonUni = NonUni $\cup \{K_i\}$.\\
\keyw{return} all elements in Uni.


\end{tabbing}
\end{Algorithm}

It is very important to remark that, in principle, 
one also needs to keep track of lower dimensional cones present in the
decomposition for the purpose of writing the inclusion-exclusion
formula of the generating function $f(K, z)$. For example in Figure
\ref{contribution} we have counted a ray twice, and thus it needs
to be removed.

\begin{figure}[htpb]
\begin{center}
\includegraphics[width = 10cm]{sumcones.eps}
\caption{Contribution of lower dimensional cones} \label{contribution}
\end{center}
\end{figure}

But this is actually not necessary thanks to a {\em Brion's
polarization trick} \citep[Remark 4.3]{BarviPom}: Let $K^*$ be
the dual cone to $K$. Apply the iterative procedure above to $K^*$
instead of $K$, ignoring the lower dimensional cones. This can be done 
because once we polarize the result back, the
contribution of the lower dimensional cones is zero with respect to
the valuation that assigns to an indicator function its generating
function counting the lattice points \citep[Corollary 2.8]{BarviPom}. In the current implementation of {\tt LattE} we do
the following:

\begin{enumerate}
\item Find the vertices of the polytope and their defining supporting cones. 

\item Compute the polar cone to each of the cones.

\item Apply the Barvinok decomposition to each of the polars.

\item Polarize back the cones to obtain a decomposition, into
full-dimensional unimodular cones, of the original supporting cones.

\item Recover the generating function of each cone and, by Brion's
theorem, of the whole polytope.
\end{enumerate}

Here is an example of how we carry out the decomposition. 

\noindent
\begin{example}
Let $K$ be a cone generated by $(2, 7)^T$ and $(1, 0)^T$.  Let
$$A = \left( \begin{array}{cc}
2& 1\\
7& 0\\
\end{array}\right).$$  Then, we have $\det(A) = -7$ and $$A^{-1} = \left( \begin{array}{cc}
0 & \frac{1}{7}\\
1 & \frac{-2}{7}\\
\end{array}\right).$$
The reduced basis $A'$ of $A^{-1}$ and the unimodular matrix $U$ for 
the
transformation from $A^{-1} $ to $A'$ are:
$A' = \left( \begin{array}{cc}
\frac{1}{7} &\frac{3}{7}\\
\frac{-2}{7} & \frac{1}{7}\\
\end{array}\right),$
and 
$U = \left( \begin{array}{cc} 
0 & 1\\
1 & 3\\
\end{array}\right).$
By enumerating the column vectors, we can verify that $(\frac{-2}{7},
\frac{1}{7})^T$ is the smallest vector with respect to $l^{\infty}$ in
the lattice generated by the column vectors of $A^{-1}$.  So, we have $z =
(1, 0)^T$.  Then, we have two cones:
$$\left( \begin{array}{cc}
2 & 0\\
7 & 1\\
\end{array}\right) \mbox{ and } \left( \begin{array}{cc} 
0 & 1\\
1 & 0\\
\end{array}\right).$$
The second cone is unimodular of index $-1$ which is the 
same sign as the determinant of $A$.  Thus, Uni $= \mbox{ Uni } \cup \{\left( 
\begin{array}{cc} 
0 & 1\\
1 & 0\\
\end{array} \right)\},$ and assign to it $\epsilon = 1$.  The first cone 
has 
determinant $2$.  So, we assign $\epsilon = -1$.  Since the first cone is 
not unimodular, we have NonUni = NonUni $\cup \{\left( 
\begin{array}{cc} 
2 & 0\\
7 & 1\\
\end{array} \right)\}.$  Set $$A =\left( \begin{array}{cc}
2 & 0\\
7 & 1\\
\end{array}\right).$$
Then, we have $\det(A) = 2$ and $$A^{-1} = \left( \begin{array}{cc}
\frac{1}{2} & 0\\
\frac{-7}{2} & 1\\
\end{array}\right) \mbox{, } A' = \left( \begin{array}{cc}  
\frac{1}{2} & \frac{1}{2}\\
\frac{-1}{2} & \frac{1}{2}\\
\end{array}\right) \mbox{ and } U = \left( \begin{array}{cc}
1 & 1\\
3 & 4\\
\end{array}\right).$$ Since $\lambda = (\frac{1}{2}, \frac{-1}{2})^T$ is
the smallest vector with respect to $l^{\infty}$, we have $z = (1,
3)^T$.  So, we get two cones: $$\left( \begin{array}{cc}
2 & 1\\
7 & 3\\
\end{array}\right) \mbox{ and } \left( \begin{array}{cc} 
1 & 0\\
3 & 1\\
\end{array}\right).$$ 
The first matrix has negative determinant which is not the same sign
as the determinant of its parent matrix $A$.  Since $\epsilon_{A} =
-1$, we assign to the first cone $\epsilon = 1$ and the second one has
positive determinant, so we assign to it $\epsilon = 1$.  Since both
of them are unimodular, we take them into Uni and since NonUni is
empty, we end while loop and print all elements in Uni.

\noindent
This gives a full decomposition:
$$ \co \{ \left( \begin{array}{c}
2\\
7\\
\end{array}\right) ,
\left( \begin{array}{c}
1\\
0\\
\end{array}\right) \}$$
$$ = \ominus
\co\{ \left( \begin{array}{c}
1\\
3\\
\end{array}\right),  
\left( \begin{array}{c}
0\\
1\\
\end{array}\right) \} \oplus \co\{ \left( \begin{array}{c}
0\\
1\\
\end{array}\right), 
\left( \begin{array}{c}
1\\
0\\
\end{array}\right) \} \oplus \co\{ \left( \begin{array}{c}
2\\
7\\
\end{array}\right), 
\left( \begin{array}{c}
1\\
3\\
\end{array}\right) \}.$$ \eopf

\begin{figure}[htpb]
\begin{center}
\includegraphics[width = 10cm]{barvinok_example.eps}
\caption{Example of Barvinok's decomposition}
\end{center}
\end{figure}

\end{example}

From the previous example, we notice that the determinant of each cone
gets much smaller in each step. This is not an accident as Theorem
\ref{break} guarantees that the cardinality of the index set $I$ of
cones in the decomposition is bounded polynomially in terms of the
determinant of the input matrix.  We have looked experimentally at how
many levels of iteration are necessary to carry out the
decomposition. We observed experimentally that it often grows linearly with
the dimension. We tested two kinds of instances. We used random square
matrices whose entries are between 0 and 9, thinking of their columns
as the generators of a cone centered at the origin. We tested from $2
\times 2$ matrices all the way to $8 \times 8$ matrices, and we tested
fifteen random square matrices for each dimension. We show the results
in Table \ref{randomcones}. For computation, we used a $1$ GHz Pentium
PC machine running Red Hat Linux.

\begin{table} 
\begin{center}

\begin{tabular}{|c|c|c|c|c|} \hline
Dimension & Height of tree  & \# of cones & $|$ determinant$| $& Time 
(seconds)\\ \hline 
2 &  1.33 & 2.53 & 11.53 & 0 \\ \hline
3 & 2.87 & 12.47 & 55.73 & 0.005\\ \hline
4 & 3.87 & 65.67  & 274.667& 0.153  \\ \hline
5 & 5.87 & 859.4 & 3875.87 & 0.25 \\ \hline
6 & 7.47 & 10308 & 19310.4 & 3.67 \\ \hline
7 & 8.53 & 91029.4 & 72986.3 & 41.61 \\ \hline 
8 & 10.67  & 2482647.533  & 1133094.733  & 2554.478  \\ \hline

\end{tabular}
\caption{Averages of 15 random matrices for computational experiences} 
\label{randomcones}
\end{center}
\end{table}

The second set of examples comes from the Birkhoff polytope $B_n$ of
doubly stochastic matrices \citep{schrijver}. Each vertex of
the polytope is a permutation matrix which is a $0/1$ matrix whose 
column sums and row sums are all $1$ \citep{schrijver}. We decompose
the cone with vertex at the origin and whose rays are the $n!$ 
permutation matrices. The results are reported in Table \ref{birkhoff}.


\begin{table}
\begin{center}

\begin{tabular}{|c|c|c|c|} \hline
Dimension & \# of vertices & \# of unimodular cones at a vertex cone &
Time (seconds)\\ \hline
$B_3=4$ & 6 & 3 & 0.05\\ \hline
$B_4=9$ & 24 & 16 & 0.15  \\ \hline
$B_5=16$ & 120 & 125 & 0.5 \\ \hline
$B_6=25$ & 720 & 1296 & 7.8 \\ \hline


\end{tabular}
\caption{The numbers of unimodular cones for the Birkhoff polytopes} 
\label{birkhoff}
\end{center}
\end{table} 

\subsection{From cones to rational functions and counting}
\label{limit}

Once we decompose all cones into simple unimodular cones, it is
easy to find the generating function attached to the $i$th cone
$K_i$. In the denominator there is a product of binomials of the form
$(1-z^{B_{ij}})$ where $B_{ij}$ is the $j$th ray of the cone $K_i$.
Thus the denominator is the polynomial $\prod (1-z^{B_{ij}})$. How
about the numerator? The cone $K_i$ is unimodular, thus it must have a
single monomial $z^{A_i}$, corresponding to the unique lattice point
inside the fundamental parallelepiped of $K_i$. Remember that the
vertex of $K_i$ is one of the vertices of our input polytope. If that
vertex $v$ has all integer coordinates then $A_i=v$, or else $v$ can
be written as a linear combination $\sum \lambda_j B_{ij}$ where all
the $\lambda_i$ are rational numbers and can be found by solving a
system of equations (remember the $B_{ij}$ form a vector space basis
for $\real^d$). The unique lattice point inside the parallelepiped of
the cone $K_i$ is simply $\sum \lceil \lambda_j \rceil B_{ij}$ \citep[Lemma 4.1]{BarviPom}.


Brion's theorem says the sum of the rational functions coming from the
unimodular cones at the vertices is a polynomial with one monomial per
lattice point inside the input polytope. One might think that to
compute the number of lattice points inside of a given convex
polyhedron, one could directly substitute the value of $1$ at each of
the variables. Unfortunately, $(1,1,\dots,1)$ is a singularity of all
the rational functions. Instead we discuss the method used in {\tt
LattE} to compute this value, which is different from that presented
by Barvinok \citep{BarviPom}. The typical generating function of
lattice points inside a unimodular cone forms:

$$ E[i] \frac{z^{A_i}}{\prod (1-z^{B_{ij}})},$$

where $z^a$ is monomial in $d$ variables, each $A_i$ (cone vertex)
and $B_{ij}$ (a generator of cone $i$) are integer vectors of length
$d, i$ ranges over all cones given, $j$ ranges over the generators of
cone $i$, and $E[i]$ is 1 or -1. Adding these rational functions and
simplifying would yield the polynomial function of the lattice point
of the polytope. Now this is practically impossible as the number of
monomials is too large. But calculating the number of monomials in
this polynomial is equivalent to evaluating the limit as $z_i$ goes to
$1$ for all $i$. We begin by finding an integer vector $\lambda$ and making the
substitution $z_i \rightarrow t^{\lambda_i}$. This is with the
intention of obtaining a univariate polynomial. To do this, $\lambda$
must be picked such that there is no zero denominator in any cone
expression, i.e. no dot product of $\lambda$ with a $B_{ij}$ can be
zero. Barvinok showed that such a $\lambda$ can be picked in polynomial time
by choosing points on the moment curve. Unfortunately, this method
yields large values in the entries of $\lambda$. Instead we try random
vectors with small integer entries, allowing small increments if
necessary, until we find  $\lambda$. Since we are essentially trying
to avoid a measure zero set, this process terminates very quickly in
practice.

After substitution, we have expressions of the form $\pm
t^{N_i}/\prod (1-t^{D_{ij}})$, where $N_i$ and $D_{ij}$ are
integers. Notice that this substitution followed by summing these
expressions yields the same polynomial as would result from first
summing and then substituting. This follows from the fact that we can
take Laurent series expansions, and the sum of Laurent series is equal
to the Laurent series of the sum of the original expressions.

Also, note that we have the following identity:
$$\sum_{\alpha \in P \cap \Z^d} z^{\alpha} = \sum_{i=1}^{\# \ of \
cones} E[i] \frac{z^{A_i}}{\prod(1-z^{B_{ij}})}.$$ After substitution
we have the following univariate (Laurent) polynomial such that:
$$\sum_{\alpha \in P \cap \Z^d} t^{\sum_{i=1}^d \lambda_i \alpha_i}=\sum^{\# \ of \ cones}_{i=1} E[i] \frac{t^{N_i}}{\prod (1-t^{D_{ij}})}.$$

With the purpose of avoiding large exponents in the numerators, we
factor out a power of $t$, say $t^c$. Now we need to evaluate the sum
of these expressions at $t=1$, but we cannot evaluate these
expressions directly at $t=1$ because each has a pole there. Consider
the Laurent expansion of the sum of these expressions about $t=1$.
The expansion must evaluate at $t=1$ to the finite number
$\sum_{\alpha \in P \cap \Z^d} 1$. It is a Taylor expansion and its
value at $t=1$ is simply the constant coefficient.  If we expand each
expression about $t=1$ individually and add them up, it will yield the
same result as adding the expressions and then expanding (again the
sum of Laurent expansions is the Laurent expansion of the sum of the
expressions). Thus, to obtain the constant coefficient of the sum, we
add up the constant coefficients of the expansions about $t=1$ of each
summand.  Computationally, this is accomplished by substituting
$t=s+1$ and expanding about $s=0$ via a polynomial division.  Summing
up the constant coefficients with proper accounting for $E[i]$ and
proper decimal accuracy yields the desired result: the number of
lattice points in the polytope. Before the substitution $t=s+1$ we
rewrite each rational function in the sum (recall $t^c$ was factored
to keep exponents small);

$$\sum E[i] \frac{t^{N_i-c}}{\prod (1-t^{D_{ij}})}=\sum E'[i]
\frac{t^{N'_i}}{\prod (t^{D'_{ij}}-1)},$$ involves in such a way that
$D'_{ij}>0$ for all $i,j$. This requires that the powers of $t$ at each
numerator to be modified, and the sign $E[i]$ is also adjusted to $E'[i]$. Then
the substitution $t=s+1$ yields
$$ \sum E'[i] \frac{(1+s)^{N'_i}}{\prod ((1+s)^{D'_{ij}}-1)},$$ 
where it is evident that, in each summand, the pole $s=0$ has an order
equal to the number of factors in the denominator. This is the same as
the number of rays in the corresponding cone and we denote this number
by $d$.


Thus the summand for cone $i$ can be rewritten as
$E'[i]s^{-d}P_i(s)/Q_i(s)$ where $P_i(s)=(1+s)^{N_i}$ and
$Q_i(s)=\prod^d ((1+s)^{D'_{ij}}-1)/s)$.  $P_i(s)/Q_i(s)$ is a
Taylor polynomial whose $s^d$ coefficient is the contribution we are
looking for (after accounting for the sign $E'[i]$ of course). The
coefficients of the quotient $P_i(s)/Q_i(s)$ can be obtained
recursively as follows: Let $Q_i(s) =b_0 + b_1 s + b_2 s^2 + \ldots $
and $P_i(s) = a_0 + a_1 s+ a_2 s^2 + \ldots $ and let
$\frac{P_i(s)}{Q_i(s)} = c_0 + c_1 s + c_2 s^2 +\ldots $. Therefore,
we want to obtain $c_d$ which is the coefficient of the constant term
of $P_i/Q_i$.  So, how do we obtain $c_d$ from $Q_i(s)$ and $P_i(s)$?
We obtain this by the following recurrence relation:

$$ c_0 = \frac{a_0}{b_0}, $$
$$ c_k = \frac{1}{b_0}(a_k - b_1 c_{k-1} - b_2c_{k-2} - \ldots - b_k
c_0) \mbox{ for } k = 1, 2, \ldots.$$

In order to obtain $c_d$, only the coefficients $a_0$, $a_1, \ldots ,
a_d$ and $b_0$, $b_1, \ldots , b_d$ are required.

\begin{example}
(A triangle). Let us consider three points in $2$
dimensions such that $V_1 = (0, 1)$, $V_2 = (1, 0)$, and $V_3 = (0,
0)$.  Then, the convex hull of $V_1$, $V_2$, and $V_3$ is a triangle
in $2$ dimensions.  We want to compute the number of lattice points by
using the residue theorem.  Let $K_i$ be the vertex cone at $V_i$ for
$i = 1, 2, 3$.  Then, we have the rational functions:

$$f(K_1, (x, y)) = \frac{y}{(1-y^{-1})(1 - x y^{-1})}, f(K_2, (x, y)) = \frac{x}{(1-x^{-1})(1 - x^{-1} y)},$$ $$ f(K_3, (x, y)) = \frac{1}{(1-x)(1-y)}.$$

We choose a vector $\lambda$ such that the inner products of $\lambda$
and the generators of $K_i$ are not equal to zero.  We choose $\lambda
= (1, -1)$ in this example. Then, reduce multivariate to univariate
with $\lambda$, so that we have:

$$f(K_1, t) = \frac{t^{-1}}{(1-t)(1 - t^{2})}, f(K_2, t) = \frac{t}{(1-t^{-1})(1 - t^{-2})}, f(K_3, t) = \frac{1}{(1-t)(1-t^{-1})}.$$ 

We want to have all the denominators to have
positive exponents.  We simplify them in order to eliminate negative
exponents in the denominators with simple algebra.  Then, we have:

$$f(K_1, t) = \frac{t^{-1}}{(1-t)(1 - t^{2})}, f(K_2, t) =
\frac{t^4}{(1-t)(1 - t^{2})}, f(K_3, t) = \frac{-t}{(1-t)(1-t)}.$$ 

We factor out $t^{-1}$ from each rational function, so that we obtain:

$$f(K_1, t) = \frac{1}{(1-t)(1 - t^{2})}, f(K_2, t) =
\frac{t^5}{(1-t)(1 - t^{2})}, f(K_3, t) = \frac{-t^2}{(1-t)(1-t)}.$$  

We substitute $t = s+1$ and simplify them to the form 
$\frac{P(s)}{s^d Q(s)}$:

$$f(K_1, s) = \frac{1}{s^2(2 + s)}, f(K_2, s) = \frac{1+5s+10s^2+10s^3+5s^4+s^5}{s^2(2 + s)},$$ $$ f(K_3, s) = \frac{-(1+2s+s^2)}{s^2}.$$

Now we use the recurrence relation to obtain the coefficient of
the constant terms.  Then, for $f(K_1, s)$, we have $c_2 = \frac{1}{8}$.
For $f(K_2, s)$, we have $c_2 = \frac{31}{8}$.  For $f(K_3, s)$, we have
$c_2 = - 1$.  Thus, if we sum up all these coefficients, we have $3$,
which is the number of lattice points in this triangle. \eopf

\end{example}

{\tt LattE} produces the sum of rational functions which converges to the generating function of the lattice points of
an input polytope. This generating function is a multivariate
polynomial of finite degree. As we saw in Subsection \ref{limit}
it is possible to count the number of lattice points without expanding
the rational functions into the sum of monomials. Suppose that instead
of wanting to know the number of lattice points we simply wish to {\em
decide} whether there is one lattice point inside the polytope or
not. The integer feasibility problem is an important and difficult
problem \citep{aardaletal1,schrijver}. Obviously, one can simply
compute the residues and then if the number of lattice points is
non-zero, clearly, the polytope has lattice points. 

Before we end our description of {\tt LattE}, we must comment on how
we deal with polytopes that are not full-dimensional (e.g.
transportation polytopes).  Given the lower-dimensional polytope
$P=\{x\in\R^n:Ax=a, Bx\leq b\}$ with the $d\times n$ matrix $A$ of
full row-rank, we will use the equations to transform $P$ into a
polytope $Q=\{x\in\R^{n-d}:Cx\leq c\}$ in fewer variables, whose
integer points are in one-to-one correspondence to the integer points
of $P$. This second polytope will be the input to the main part of
{\tt LattE}. The main idea of this transformation is to find the
general integer solution $x=x_0+\sum_{i=1}^{n-d} \lambda_i g_i$ to
$Ax=a$ and to substitute it into the inequalities $Bx\leq b$, giving a
new system $Cx\leq c$ in $n-d$ variables
$\lambda_1,\ldots,\lambda_{n-d}$.

It is known that the general integer solution $Ax=a$ can be found
via the Hermite normal form $H=(R|0)$ of $A$ \citep{schrijver}.
Here, $R$ is a lower-triangular matrix and $H=AU$ for some
unimodular matrix $U$. Moreover, as $A$ is supposed to have full
row-rank, $R$ is a non-singular $d\times d$ matrix. Let $U_1$ be
the matrix consist of the first $d$ columns of $U$ and $U_2$
consisting of the remaining $n-d$ columns of $U$. Now we have
$AU_1=R$ and $AU_2=0$ and the columns of $U_2$ give the generators
$\{g_1,\ldots,g_{n-d}\}$ of the integer null-space of $A$. Thus,
it remains to determine a special integer solution $x_0$ to
$Ax=a$.

To do this, first find an integer solution $y_0$ to $Hy=(R|0)y=a$,
which is easy due to the triangular structure of $R$. With
$x_0=Uy_0$, we get $Ax_0=AUy_0=Hy_0=a$ and have found all pieces
of the general integer solution $x=x_0+\sum_{i=1}^{n-d} \lambda_i
g_i$ to $\{x\in\Z^n:Ax=a\}$. 


\section{Computational experience and performance for counting} \label{experiments}

One can download {\tt LattE} via a web page
\url{www.math.ucdavis.edu/~latte}. You can also find there the files of all the
experiments presented in this section.  At the moment we
have been able to handle polytopes of dimension 30 and several
thousands vertices.  It is known that the theoretical upper bound of
the number of unimodular cones is $2^{dh}$, where $h = \lfloor
\frac{\log \log 1.9 - \log \log D }{\log(d-1/d)} \rfloor$ and where
$D$ is the volume of the fundamental parallelepiped of the input cone
\citep{bar}. If we fix the dimension this upper bound becomes
polynomial time. Unfortunately, if we do not fix the dimension, this
upper bound becomes exponential.  In practice this might be costly and
some families of polytopes have large numbers of unimodular cones.  The
cross polytope family, for instance, has many unimodular cones and behaves badly.  For example, for the cross polytope in $6$
dimensions, with cross6.ine input file \citep{fukuda}, {\tt LattE} took
$147.63$ seconds to finish computing.  The number of lattice points of
this polytope is obviously $13$. Also, for the cross polytope in $8$
dimensions, with cross8.ine input file \citep{fukuda}, {\tt LattE} took
85311.3 seconds to finish computing, even though this polytope has
only $16$ vertices and the number of lattice points of this polytope
is $17$.  For all computations, we used a $1$ GHz Pentium PC machine
running Red Hat Linux.  


Here is a short description of how to use {\tt
LattE}. 

For computations involving a polytope $P$ described by a
system of inequalities $Ax\leq b$, where $A\in\Z^{m\times d}$, 
$A=(a_{ij})$, and $b\in\Z^m$, the {\tt LattE} readable input file
would be as follows: 
\begin{verbatim}
m d+1
b  -A
\end{verbatim}

\noindent
\textbf{EXAMPLE.}  
Let $P=\{(x,y): x\leq 1, y\leq 1, x+y\leq 1, x\geq 0, y\geq 0\}$.
Thus
\[
\begin{array}{ccc}
A=\left(
\begin{array}{rr} 
 1 &  0 \\ 
 0 &  1 \\ 
 1 &  1 \\
-1 &  0 \\ 
 0 & -1 \\ 
\end{array} 
\right) 
& , &
b = \left( 
\begin{array}{r} 
1 \\ 
1 \\ 
1 \\ 
0 \\
0 \\ 
\end{array} 
\right)
\end{array}
\]
and the {\tt LattE} input file would be as such:
\begin{verbatim}
5 3
1 -1  0
1  0 -1
1 -1 -1
0  1  0
0  0  1
\end{verbatim}

In {\tt LattE}, polytopes are represented by {\bf linear constraints},
i.e. equalities or inequalities. By default a constraint is an
inequality of type $ax\leq b$ unless we specify, by using a single
additional line, the line numbers of constraints that are linear
equalities. 

\noindent
\textbf{EXAMPLE.}  
Let $P$ be as in the previous example, but require $x+y=1$ instead of
$x+y\leq 1$, thus, 
$P=\{(x,y): x\leq 1, y\leq 1, x+y=1, x\geq 0, y\geq 0\}$.
Then the {\tt LattE} input file that describes $P$ would be as such:
\begin{verbatim}
5  3
1 -1  0
1  0 -1
1 -1 -1
0  1  0
0  0  1
linearity 1 3
\end{verbatim}
The last line states that among the $5$ inequalities one is to be
considered an equality, the third one.

For bigger examples it quickly becomes cumbersome to state all
nonnegativity constraints for the variables one by one. Instead, you
may use another short-hand.

\noindent
\textbf{EXAMPLE.}  
Let $P$ be as in the previous example, then the {\tt LattE} input file
that describes $P$ could also be described as such: 
\begin{verbatim}
3  3
1 -1  0
1  0 -1
1 -1 -1
linearity 1 3
nonnegative 2 1 2
\end{verbatim}
The last line states that there are two nonnegativity constraints and
that the first and second variables are required to be nonnegative. 
{\bf NOTE} that the first line reads ``3 3'' and not ``5 3'' as above! 


We now report on computations with convex rational polytopes. We used
a $1$ GHz Pentium PC machine running Red Hat Linux. We begin with the
class of {\em multiway contingency tables}. A {\em $d$-table of size
$(n_1,\dots,n_d)$} is an array of non-negative integers
$v=(v_{i_1,\dots,i_d})$, $1\leq i_j\leq n_j$. For $0\leq m <d$, an
{\em $m$-marginal} of $v$ is any of the $d\choose m$ possible
$m$-tables obtained by summing the entries over all but $m$
indices. For instance, if $(v_{i,j,k})$ is a $3$-table then its
$0$-marginal is
$v_{+,+,+}=\sum_{i=1}^{n_1}\sum_{j=1}^{n_2}\sum_{k=1}^{n_3}
v_{i,j,k}$, its $1$-marginals are
$(v_{i,+,+})=(\sum_{j=1}^{n_2}\sum_{k=1}^{n_3}v_{i,j,k})$ and likewise
$(v_{+,j,+})$, $(v_{+,+,k})$, and its $2$-marginals are
$(v_{i,j,+})=(\sum_{k=1}^{n_3} v_{i,j,k})$ and likewise $(v_{i,+,k})$,
$(v_{+,j,k})$.

Such tables appear naturally in statistics and operations research
under various names such as {\em multi-way contingency tables}, or
{\em tabular data}.  We consider the {\em table counting problem}:
{\em given a prescribed collection of marginals, how many $d$-tables
are there that share these marginals?}  Table counting has several
applications in statistical analysis, in particular for independence
testing, and has been the focus of much research (see \citep{DG} and
the extensive list of references therein).  Given a specified
collection of marginals for $d$-tables of size $(n_1,\dots,n_d)$
(possibly together with specified lower and upper bounds on some of
the table entries) the associated {\em multi-index transportation
polytope} is the set of all non-negative {\em real valued} arrays
satisfying these marginals and entry bounds. The counting problem can
be formulated as that of counting the number of integer points in the
associated multi-index transportation polytope.  We begin with a small
example of a three-dimensional table of format $ 2 \times 3 \times 3$
given below. The data displayed in Table \ref{census} have been
extracted from the 1990 decennial census and is used in
\citep{FMMS}. For the $2$-marginals implied by these data we get the
answer of 441 in less than a second.

\begin{table}[ht]
\begin{center}
\begin{tabular}{|l|c|c|c||c|}
\multicolumn{5}{c}{Gender $=$ Male} \\ 
\multicolumn{1}{c}{} & \multicolumn{4}{c}{Income Level} \\ \cline{2-5}
\multicolumn{1}{c|}{Race} & $\le \$10,000$ & $> \$10000 \mbox{ and} \le
\$25000$ & $> \$25000 $ & Total \\ \hline
White & 96 & 72 & 161 & 329 \\ \hline
Black & 10 & 7 & 6 & 23 \\ \hline
Chinese & 1 & 1 & 2 & 4 \\ \hline\hline
Total & 107 & 80 & 169 & 356 \\ \hline
\multicolumn{5}{c}{} \\[11pt] 
\multicolumn{5}{c}{Gender $=$ Female} \\ 
\multicolumn{1}{c}{} & \multicolumn{4}{c}{Income Level} \\ \cline{2-5}
\multicolumn{1}{c|}{Race} & $\le \$10,000$ & $> \$10000 \mbox{ and} \le
\$25000$ & $> \$25000 $ & Total \\ \hline
White & 186 & 127 & 51 & 364 \\ \hline
Black & 11 & 7 & 3 & 21 \\ \hline
Chinese & 0 & 1 & 0 & 1 \\ \hline\hline
Total & 197 & 135 & 54 & 386 \\ \hline
\end{tabular}
\end{center}
\caption{Three-way cross-classification of gender, race, and income for a
selected U.S. census tract. {\em Source}: 1990 Census Public Use Microdata
Files.}\label{census}
\end{table}

 We present now an example of a $3 \times 3 \times 3$ table with
fairly large $2$-marginals. They are displayed in Table \ref{3by3by3}.
{\tt LattE} took only $19.67$ seconds of CPU time.  The number of
lattice points inside of this polytope is
$ 2249847900174017152559270967589010977293 $.


\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|} \hline
      164424  &    324745  &    127239  \\ \hline
      262784  &    601074  &   9369116  \\ \hline
      149654  &   7618489  &   1736281  \\ \hline
\end{tabular}
\end{center}


\begin{center}
\begin{tabular}{|c|c|c|} \hline
      163445   &    49395 &     403568  \\ \hline
     1151824   &   767866 &    8313284  \\ \hline
     1609500   &  6331023 &    1563901  \\ \hline


\end{tabular}
\end{center}

\begin{center}
\begin{tabular}{|c|c|c|} \hline
      184032  &    123585 &     269245  \\ \hline
      886393  &   6722333 &     935582  \\ \hline
     1854344  &    302366 &    9075926  \\ \hline

\end{tabular}
\end{center}
\caption{2-Marginals for the $3 \times 3 \times 3$ example.} \label{3by3by3}
\end{table}

 Next we present an example of a $3 \times 3 \times 4$ table with
large 2-marginals. The $2$-marginals are displayed in Table
\ref{3by3by4}.  The CPU time for this example was $44$ minutes $42.22$
seconds. The number of lattice points inside of this polytope is

$ 4091700129572445106288079361219676736812805058988286839062994.$


\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|c|} \hline
      273510  &    273510  &    273510  &    191457  \\ \hline
      273510  &    273510  &    547020  &    191457  \\ \hline
      273510  &    547020  &    273510  &    191457  \\ \hline
\end{tabular}
\end{center}


\begin{center}
\begin{tabular}{|c|c|c|} \hline

      464967  &    273510  &    273510   \\ \hline 
      547020  &    273510  &    464967  \\ \hline
      410265  &    601722  &    273510  \\ \hline

\end{tabular}
\end{center}


\begin{center}
\begin{tabular}{|c|c|c|} \hline

      273510  &    273510  &    273510   \\ \hline
      410265  &    547020  &    136755   \\ \hline
      547020  &    136755  &    410265   \\ \hline
      191457  &    191457  &    191457   \\ \hline

\end{tabular}
\end{center}
\caption{2-Marginals for the $3 \times 3 \times 4$ example.} \label{3by3by4}
\end{table} 


The next family of examples are some hard knapsack-type
problems. Suppose we have a set of positive relatively prime integers
$\{a_1, a_2, \ldots , a_d \}$. Denote by $a$ the vector $(a_1, a_2,
\ldots , a_d)$.  Consider the following problem: does there exist a
non-negative integral vector $x$ satisfying $ax = a_0$ for some
positive integer $a_0$?  We take several examples from
\citep{aardaletal3} which have been found to be extremely hard to solve
by commercial quality branch-and-bound software. This is very
surprising since the number of variables is at most 10. It is not very
difficult to see that if the right-hand-side value $a_0$ is large
enough, the equation will surely have a non-negative integer solution.
The {\em Frobenius number} for a knapsack problem is the largest value
$a_0$ such that the knapsack problem is infeasible. Aardal and Lenstra
\citep{aardaletal3} solved them using the reformulation in
\citep{aardaletal1}.  Their method works significantly better than
branch-and-bound using {\tt CPLEX} 6.5. Here we demonstrate that our
implementation of Barvinok's algorithm is fairly fast and, on the
order of seconds, we resolved the first 15 problems in Table 1 of
\citep{aardaletal3} and verified all are infeasible except {\bf prob9},
where there is a mistake. The vector $(3480,1,4,4,1,0,0,0,0,0)$ is a
solution to the right-hand-side $13385099.$ In fact, using {\tt LattE}
we know that the exact number of solutions is $838908602000$.  For
comparison we named the problems exactly as in Table 1 of
\citep{aardaletal3}. We present our results in Table \ref{knapsack}.
It is very interesting to know the
number of lattice points if we add $1$ to the Frobenius number for
each problem.  In Table \ref{knap2}, we find the number of solutions
if we add $1$ to the Frobenius number on each of the (infeasible)
problems. The speed is practically the same as in the previous
case. In fact the speed is the same regardless of the right-hand-side
value $a_0$.


\begin{sidewaystable}
\centering
\begin{tabular}{|c|cccccccccc|c|c|} \hline

&&&&&&&&& &&Frobenius \# & Time (m, s)\\ \hline
cuww1&12223 &12224 &36674 &61119 &85569& & & & & & 89643481 & 0.55s\\ 
\hline
cuww2&12228 &36679 &36682 &48908 &61139 &73365 & & & & &89716838& 1.78s  
\\ \hline
cuww3&12137 &24269 &36405 &36407 &48545 &60683 & & & & &58925134 & 1.27s\\ 
\hline
cuww4&13211 &13212 &39638 &52844 &66060 &79268 &92482 & & & &104723595& 
2.042s \\ \hline
cuww5&13429 &26850 &26855 &40280 &40281 &53711 &53714 &67141 & & & 
45094583& 16.05s\\ \hline
pro1& 25067 &49300 &49717 &62124 &87608 &88025 &113673 &119169 & & & 
33367335& 47.07s  \\ \hline 
prob2& 11948 &23330 &30635 &44197 &92754 &123389 &136951 &140745 & & 
&14215206 & 1m0.58s \\ \hline
prob3  &39559 &61679 &79625 &99658 &133404 &137071 &159757 &173977 & & 
&58424799 & 1m28.3s \\ \hline
prob4& 48709 &55893 &62177 &65919 &86271 &87692 &102881 &109765 & & 
&60575665 & 59.04s \\ \hline
prob5& 28637 &48198 &80330 &91980 &102221 &135518 &165564 &176049 & & & 
62442884& 1m41.78s \\ \hline
prob6 &20601 &40429 &40429 &45415 &53725 &61919 &64470 &69340 &78539 
&95043 &
22382774 & 3m45.86s  \\ \hline
prob7 &18902 &26720 &34538 &34868 &49201 &49531 &65167 &66800 &84069 
&137179 &
27267751 &  2m57.64s\\  \hline
prob8& 17035 &45529 &48317 &48506 &86120 &100178 &112464 &115819 &125128 
&129688
& 21733990 & 8m29.78s\\ \hline
prob10 & 45276 &70778 &86911 &92634 &97839 &125941 &134269 &141033 &147279 
&153525 &106925261 & 4m24.67s \\ \hline
\end{tabular}
\centering
\caption{Infeasible knapsack problems.} \label{knapsack}
\end{sidewaystable}


\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|}\hline
problem & RHS & \# of lattice points. \\ \hline
cuww1& 89643482 & 1 \\ \hline
cuww2 & 89716839 & 1 \\ \hline
cuww3 & 58925135 & 2 \\ \hline
cuww4 & 104723596 & 1 \\ \hline
cuww5 & 45094584  & 1 \\ \hline
prob1 & 33367336 & 859202692 \\ \hline
prob2 & 14215207  & 2047107 \\ \hline
prob3 & 58424800   &  35534465752\\ \hline
pro4 & 60575666   & 63192351 \\ \hline
pro5 & 62442885   & 21789552314 \\ \hline
pro6 & 22382775   &  218842 \\ \hline
pro7 & 27267752   & 4198350819898 \\ \hline
pro8 & 21733991   & 6743959 \\ \hline
pro10 & 106925262  & 102401413506276371 \\ \hline
\end{tabular}
\caption{The number of lattice points if we add 1 to the Frobenius
number.}\label{knap2}
\end{center}
\end{table}

Already counting the lattice points of large width convex polygons is
a non-trivial task if one uses brute-force enumeration (e.g. list one
by one the points in a bounding box of the polygon and see whether it
is inside the polygon). Here we experiment with very large convex {\em
almost} regular $n$-gons. Regular $n$-gons cannot have rational
coordinates, but we can approximate them to any desired accuracy by
rational polygons. In the following experiment we take regular
$n$-gons, from $n=5$ to $n=12$ centered at the origin (these have only
a handful of lattice points). We take a truncation of the coordinates
up to 3, 9, and 15 decimal digits, then we multiply by a large enough
power of 10 to make those vertex coordinates integral and we count the
number of lattice points in the dilation. All experiments take less
than a second.

{\tiny
\begin{table} 
\begin{center}
{\small

\begin{tabular}{|c|c|c|c|} \hline
 & $10^3$ (seconds) & $10^9$ (seconds)& $10^{15}$ (seconds) \\ \hline
5gon & 2371673(0.136) & 2377641287748905186(0.191) & 2377641290737895844565559026875(0.289) \\ \hline
6gon & 2596011(0.153) & 2598076216000000011(0.193) & 2598076211353321000000000000081(0.267) \\ \hline
7gon & 2737110(0.175) & 2736410188781217941(0.318) & 2736410188638105174143840143912(0.584s) \\ \hline
8gon & 2820021(0.202) & 2828427120000000081(0.331) & 2828427124746200000000000000201(0.761) \\ \hline
9gon & 2892811(0.212) &  2892544245156317460(0.461) & 2892544243589428566861745742966(0.813) \\ \hline
10gon & 2931453(0.221) &  2938926257659276211(0.380) & 2938926261462380264188126524437(0.702) \\ \hline
11gon & 2974213(0.236) & 2973524496796366173(0.745) &  2973524496005786351949189500315(1.858) \\ \hline
12gon & 2997201(0.255) & 3000000004942878881(0.466) & 3000000000000005419798779796241(0.696) \\ \hline
\end{tabular}
\caption{The numbers of the approximated regular polygons. We show
the number of lattice points in different dilation factors (powers
of ten) and time of computation.}
}
\end{center}
\end{table}
}

The next two sets of examples are families that have been studied quite
extensively in the literature and provide us with a test for speed.
In the first case we deal with {\em two-way contingency tables.} The
polytope defined by a two-way contingency table is called the {\em
transportation polytope}.  We present the results in Table
\ref{transport}.  The second family consists of flow polytopes for the
complete $4$-vertex and the complete $5$-vertex tournaments (directed complete  graphs).  Consider the
directed complete graph $K_l$ for $l \in \N$ and $l \geq 3$.  We
assign a number to each node of the graph.  Then, we orient the arcs
 from the node of smaller index to the node of bigger index. 
Let $N$ be the node set of the
complete graph $K_l$, let $w_i$ be a weight assigned to node $i$ for $i =
1, 2, \ldots , l$, and let $A$ be the arc set of $K_l$.  Then, we have
the following constraints, with as many variables as arcs:

$$\sum_{(j,i) \mbox{arc enters} i} x_{ji} - \sum_{(i,j) \mbox{arc has tail} i} x_{ij} = w_i,$$
$$ x_{ij} \geq 0 \mbox{ } \forall  (i,j) \in
A.$$ These equalities and inequalities define a polytope and this
polytope is the special case of a {\em flow polytope}. The results for the complete
graphs $K_4$ and $K_5$, with different weight vectors, are shown in Tables 
\ref{k4} and \ref{k5} respectively.


These two families of polytopes have been studied by several authors\\
\citep{baldonideloeravergne, Beck2, deloerasturmfels, mount} and thus
are good for testing the performance of {\tt LattE}. We used several
examples of transportation polytopes, as presented in the table
below. In general, {\tt LattE} runs at comparable performance to the
software of \citep{baldonideloeravergne, Beck2} for generic vectors
$(a,b)$ but is slower for degenerate inputs (those that do not give a
simple polytope). The reason seems to be that at each non-simplex vertex {\tt
LattE} needs to triangulate each cone which takes considerable time in
problems of high dimension.


\begin{sidewaystable}

\centering
{\tiny

\begin{tabular}{|p{2.7in}|r|c|} \hline 
Margins  & \# of lattice points & Time (seconds)\\ \hline 
[220, 215, 93, 64], \newline 
[108, 286, 71, 127] &  1225914276768514 & 1.048 \\ \hline
[109, 127, 69, 109], \newline [119, 86, 108, 101] & 993810896945891 & 1.785 \\ \hline
[72, 67, 47, 96], \newline [70, 70, 51, 91] &  25387360604030 &  1.648\\ \hline
[179909, 258827, 224919, 61909], \newline [190019, 90636, 276208, 168701] &  13571026063401838164668296635065899923152079 &  1.954  \\ \hline
[229623, 259723, 132135, 310952],\newline  [279858, 170568, 297181, 184826] &646911395459296645200004000804003243371154862 &  1.765 \\ \hline
[249961, 232006, 150459, 200438],\newline  [222515, 130701, 278288, 201360] & 319720249690111437887229255487847845310463475 &  1.854  \\ \hline 
[140648, 296472, 130724, 309173], \newline [240223, 223149, 218763, 194882]&322773560821008856417270275950599107061263625 &  1.903 \\ \hline
[65205, 189726, 233525, 170004],\newline  [137007, 87762, 274082, 159609]  & 6977523720740024241056075121611021139576919 &1.541 \\ \hline
[251746, 282451, 184389, 194442], \newline [146933, 239421, 267665, 259009] & 861316343280649049593236132155039190682027614 & 1.880 \\ \hline
[138498, 166344, 187928, 186942], \newline [228834, 138788, 189477, 122613] & 63313191414342827754566531364533378588986467 & 1.973 \\ \hline
[20812723, 17301709, 21133745, 27679151], \newline [28343568, 18410455,
19751834, 20421471] & 
665711555567792389878908993624629379187969880179721169068827951
& 2.917\\ 
\hline

[15663004, 19519372, 14722354, 22325971], \newline [17617837, 25267522,
20146447, 9198895] &
63292704423941655080293971395348848807454253204720526472462015 &
3.161
\\
\hline

[13070380, 18156451, 13365203, 20567424], \newline [12268303, 20733257,
17743591, 14414307] &
43075357146173570492117291685601604830544643769252831337342557 &
2.990
\\
\hline

\end{tabular}
\caption{Testing for $4 \times 4$ transportation polytopes.}

}
\label{transport}
\end{sidewaystable} 


\begin{table}
\begin{center}

\begin{tabular}{|l|r|c|} \hline 
Weights on nodes  & \# of lattice points & Time (seconds)\\ \hline 
[-6, -8, 5, 9] &  223 & 0.288 \\ \hline
[-9, -11, 12, 8] & 330 & 0.286 \\ \hline
[-1000, -1, 1000, 1] &  3002 & 0.287 \\ \hline
[-4383, 886, 2777, 720] &  785528058 &  0.287  \\ \hline
[-4907, -2218, 3812, 3313] & 20673947895 &  0.288 \\ \hline
[-2569, -3820, 1108, 5281] & 14100406254 &  0.282  \\ \hline 
[-3842, -3945, 6744, 1043]& 1906669380  &  0.281 \\ \hline
[-47896, -30744, 46242, 32398]  &  19470466783680  & 0.282 \\ \hline
[-54915, -97874, 64165, 88624] & 106036300535520 & 0.281 \\ \hline
[-69295, -62008, 28678, 102625] & 179777378508547 & 0.282 \\ \hline
[-3125352, -6257694, 926385, 8456661] & 34441480172695101274 & 0.509 \\
\hline

[-2738090, -6701290, 190120, 9249260] & 28493245103068590026 & 0.463 \\
\hline

[-6860556, -1727289, 934435, 7653410] & 91608082255943644656 & 0.503 \\
\hline

\end{tabular}
\caption{Testing for the complete graph $K_4$.} \label{k4}
\end{center}
\end{table}

{\tiny


\begin{table}
{\small


\begin{tabular}{|l|r|c|} \hline 
Weights on nodes  & \# of lattice points & secs\\ \hline 
[-12, -8, 9, 7, 4]&  14805 & 0.319 \\ \hline
[-125, -50, 75, 33, 67] & 6950747024 & 0.325 \\ \hline
[-763, -41, 227, 89, 488] &  222850218035543 & 0.325 \\ \hline
[-11675, -88765, 25610, 64072, 10758] &  563408416219655157542748  &  0.319 \\ \hline
[-78301, -24083, 22274, 19326, 60784] & 1108629405144880240444547243 &  0.336 \\ \hline
[-52541, -88985, 1112, 55665, 84749] & 3997121684242603301444265332 &  0.331  \\ \hline 
[-71799, -80011, 86060, 39543, 26207]& 160949617742851302259767600 & 0.316  \\ \hline
[-45617, -46855, 24133, 54922, 13417]  & 15711217216898158096466094 & 0.285 \\ \hline
[-54915, -97874, 64165, 86807, 1817] & 102815492358112722152328 & 0.277 \\ \hline
[-69295, -62008, 28678, 88725, 13900] & 65348330279808617817420057 & 0.288 \\ \hline
[-8959393, -2901013, 85873, 533630, 11240903] &
6817997013081449330251623043931489475270 & 0.555 \\ \hline
[-2738090, -6701290, 190120, 347397, 8901863] &
277145720781272784955528774814729345461 & 0.599 \\ \hline
[-6860556, -1727289, 934435, 818368, 6835042] &
710305971948234346520365668331191134724 & 0.478 \\ \hline


\end{tabular}
\centering
\caption{Testing for the complete graph $K_5$. Time is given in seconds} \label{k5}

}
\end{table}

}

The following experiment is from \cite{mcmc2}.  
This is from an actual data in German survey.  For 2,262 German citizens they asked the following questions: if you order the following items, how you order them? 
\begin{enumerate}
\item Maintain order.
\item Give people more say.
\item Fight rising prices.
\item Protect freedom of speech.
\end{enumerate}

Then we have the data in Table \ref{S4:data}.

\begin{table}
\begin{center}
\begin{tabular}{||c|c||c|c||c|c||c|c||}\hline
1234 & 137 & 2134 & 48 & 3124 & 330 & 4123 & 21 \\ \hline
1243 & 29  & 2143 & 23 & 3142 & 294 & 4132 & 30 \\ \hline
1324 & 309 & 2314 & 61 & 3214 & 117 & 4213 & 29 \\ \hline
1342 & 255 & 2341 & 55 & 3241 & 69  & 4231 & 52 \\ \hline
1423 & 52  & 2413 & 33 & 3412 & 70  & 4312 & 35 \\ \hline
1432 & 93  & 2431 & 39 & 3421 & 34  & 4321 & 27 \\ \hline
Total & 875 & & 279 && 914 && 194 \\ \hline
\end{tabular}
\caption{Permutation $S_4$ problem from \cite{mcmc2}.}\label{S4:data}
\end{center}
\end{table}

We take sum for the number of people who picked $j \in \{1, 2, 3, 4\}$ for the $i$th order for all $j = 1, 2, 3, 4$ and for all $i = 1, 2, 3, 4$.  Then we have the following condition given in Table \ref{S4}.  

\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|c|}\hline
875 & 279 & 914 & 194 \\ \hline
746 & 433 & 742 & 341 \\ \hline
345 & 773 & 419 & 725 \\ \hline
296 & 777 & 187 & 1002 \\ \hline
\end{tabular}
\caption{Marginal conditions for the permutation problem from \cite{mcmc2}.}\label{S4}
\end{center}
\end{table}

Then we want to compute how many functions $f$ such that $f: S_4 \to \N$ and $\sum_{\sigma \in S_4} f(\sigma) = 2262$.

The solution to this problem is: Total number of functions is\\ $11606690287805167142987310121$ and 
CPU Time is 523.12 sec.

The last experiment in this section
is from Ian Dinwoodie.
Ian Dinwoodie 
communicated to us the problem of counting all $7 \times 7$
contingency tables whose entries are nonnegative integers $x_i$, with
diagonal entries multiplied by a constant as presented in Table
\ref{dinwoodie}. The row sums and column sums of the entries are given
there too. Using {\tt LattE} we obtained the exact answer {\em
{8813835312287964978894}}.

\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c||c|} \hline
  2{\color{blue}{$x_1$}}  &   $x_2$   &   $x_3$ &  $x_4$ &  $x_5$ & $x_6$ & $x_7$ & 205\\ \hline
     $x_2$   &  2{\color{blue}{$x_8$}}   &   $x_{9}$ &  $x_{10}$ &  $x_{11}$ & $x_{12}$ & $x_{13}$ & 600 \\ \hline
     $x_3$  & $x_9$   &  2{\color{blue}{$x_{14}$}}  &  $x_{15}$ &  $x_{16}$ & $x_{17}$ & $x_{18}$ & 61\\ \hline
     $x_4$   & $x_{10}$ & $x_{15}$      &  2{\color{blue}{$x_{19}$}} &  $x_{20}$ & $x_{21}$ & $x_{22}$ & 17 \\ \hline
     $x_5$  & $x_{11}$ & $x_{16}$     &$x_{20}$  & 2{\color{blue}{$x_{23}$}}  & $x_{24}$ & $x_{25}$ & 11 \\ \hline
     $x_6$  &$x_{12}$  & $x_{17}$    & $x_{21}$  & $x_{24}$  & 2{\color{blue}{$x_{26}$}} & $x_{27}$ & 152 \\ \hline
     $x_7$ &  $x_{13}$ & $x_{18}$    & $x_{22}$  & $x_{25}$  & $x_{27}$  & 2{\color{blue}{$x_{28}$}} & 36 \\ \hline
    205 & 600 & 61 & 17 & 11 & 152 &36 & 1082 \\ \hline
\end{tabular}
\caption{The conditions for retinoblastoma RB1-VNTR genotype data from
 the Ceph database.} \label{dinwoodie} 
\end{center}
\end{table}

\section{New Ehrhart (quasi-)polynomials}\label{formulas}
 
Given a rational polytope $P \subset \R^d$, the function \[ i_P(t) :=
  \# \left( t P \cap \Z^d \right),\] for a positive integer $t$, was
  first studied by E. Ehrhart \citep{Ehrhart} and has received a lot of
  attention in combinatorics.  It is known to be a polynomial when all
  vertices of $P$ are integral and it is a quasi-polynomial for
  arbitrary rational polytopes. It is called the {\em Ehrhart
  quasi-polynomial} in honor of its discoverer 
  \citep[Chapter 4]{stanley}.  A function $f: \N \rightarrow \C$ is a
  quasi-polynomial if there is an integer $N > 0$ and polynomials
  $f_0, \ldots , f_{N-1}$ such that $f(s) = f_i(s) \mbox{ if } s
  \equiv i \mod N$.  The integer $N$ is called a {\em quasi-period} of
  $f$.  Therefore, by counting the number of lattice points for
  sufficiently many dilations of a rational polytope, we can
  interpolate its Ehrhart quasi-polynomial.


Using {\tt LattE}, {\tt Maple}, and interpolation, we have calculated the Ehrhart polynomials and
  quasi-polynomials for polytopes that are slices or nice truncations
  of the unit $d$-cube.  To the best of our knowledge these values
  were not known before. For example, the $24$-cell polytope centered
  at the origin with smallest integer coordinates has Ehrhart
  polynomial $i_{24\_cell}(s) = 8s^4 + \frac{32s^3}{3} + 8s^2 +
  \frac{16s}{3} + 1.$ In Table \ref{hypersi}, we see the Ehrhart
  polynomials for the hypersimplices $\Delta(n,k)$. They are defined
  as the slice of the $n$-cube by the hyperplane of equation $\sum x_i=k$ 
  with $k \leq n$. Note that $\Delta(n,k)=\Delta(n,n-k)$ because of the
symmetries of the regular cube. The hypersimplices form one of the most famous families of $0/1$-polytopes.
It is known that hypersimplices are {\em compressed polytopes} \citep{hibi}.
This means that their Ehrhart polynomials can be recovered from the $f$-vectors
of any of their reverse lexicographic triangulations. Instead, we recovered 
them explicitly for the first time  using {\tt LattE} and interpolation.


\begin{sidewaystable}
\centering


{\tiny


\begin{tabular}{|c|c|c|} \hline
n & k  & the Ehrhart polynomial $P(s)$ \\ \hline 
4 & 1 & $\frac{s^3}{6} + s^2 + \frac{11s}{6} + 1$  \\ \hline
4 & 2 & $\frac{2x^3}{3} + 2s^2 + \frac{7s}{3} + 1$ \\ \hline
5 & 1 & $\frac{s^4}{24} + \frac{5s^3}{12} + \frac{35s^2}{24} +\frac{25s}{12} + 1 $\\ \hline
5 & 2 & $\frac{11s^4}{24} + \frac{25s^3}{12} + \frac{85s^2}{24} + 
\frac{35s}{12} + 1$ \\ \hline
6 & 1 & $\frac{s^5}{120} + \frac{s^4}{8} + \frac{17s^3}{24} + \frac{15 
s^2}{8} + \frac{137s}{60} + 1$  \\ \hline
6 & 2 & $\frac{13s^5}{60} + \frac{3s^4}{2} + \frac{47s^3}{12} + 5s^2 + \frac{101 s}{30} + 1$ \\ \hline
6 & 3 & $\frac{11s^5}{20}+\frac{11s^4}{4} + \frac{23s^3}{4} + \frac{25s^2}{4} +\frac{37s}{10} + 1$   \\ \hline
7 & 1 & $\frac{s^6}{720} + \frac{7s^5}{240} + \frac{35s^4}{144} + \frac{49s^3}{48} + \frac{203s^2}{90} + \frac{49s}{20} + 1$ \\ \hline 
7 & 2  & $\frac{19s^6}{240}+ \frac{63s^5}{80} + \frac{49s^4}{16} + \frac{287s^3}{48} + \frac{763s^2}{120} + \frac{56s}{15} + 1$ \\ \hline
7 & 3 & $\frac{151s^6}{360} + \frac{161s^5}{60} + \frac{256s^4}{36} + \frac{21s^3}{2} + \frac{3199s^2}{360} + \frac{259s}{60} + 1$ \\ \hline
8 & 1 &
$\frac{s^7}{5040}+\frac{s^6}{180}+\frac{23s^5}{360}+\frac{7s^4}{18} +
\frac{967s^3}{720} + \frac{469s^2}{180}+\frac{363s}{140} + 1$ \\ \hline
8 & 2 &
$\frac{s^7}{42}+\frac{29s^6}{90}+\frac{53s^5}{30}+\frac{91s^4}{18}+
\frac{49s^3}{6}+\frac{343s^2}{45}+\frac{283s}{70}+1$ \\ \hline
8 & 3 &
$\frac{397s^7}{1680}+\frac{359s^6}{180}+\frac{281s^5}{40}+\frac{245s^4}{18}
+ \frac{1273s^3}{80}+\frac{2051s^2}{180}+\frac{2027s}{420}+1$ \\ \hline
8 & 4 & $\frac{151s^7}{315}+\frac{151s^6}{45}+\frac{463s^5}{45} +
\frac{161s^4}{9}+\frac{862s^3}{45}+\frac{574s^2}{45}+\frac{533s}{105}+1$ 
\\ \hline
9 & 1 &$\frac{s^8}{40320}+\frac{s^7}{1120}+\frac{13s^6}{960}+\frac{9s^5}{80}+\frac{1069s^4}{1920}+\frac{267s^3}{160}+ \frac{29531s^2}{10080}+\frac{761s}{280}+1$ \\ \hline
9&2& $\frac{247s^8}{40320}+\frac{121s^7}{1120} + \frac{763s^6}{960}+\frac{253s^5}{80}+\frac{14203s^4}{1920}+\frac{1667s^3}{160}+\frac{88721s^2}{10080}+\frac{1207s}{280}+1$ \\ \hline
9 & 3 & $\frac{477s^8}{4480}+\frac{1311s^7}{1120}+\frac{1731s^6}{320}+\frac{1107s^5}{80}+\frac{13899s^4}{640}+\frac{3477s^3}{160}+\frac{15419s^2}{1120}+\frac{1473s}{280}+1$ \\ \hline
9 & 4 & $\frac{15619s^8}{40320}+\frac{3607s^7}{1120}+\frac{11311s^6}{960}+\frac{1991s^5}{80}+\frac{63991s^4}{1920}+\frac{4669s^3}{160}+\frac{166337s^2}{10080}+\frac{1599s}{280}+1$ \\ \hline

10 & 1 & $\frac{s^9}{362880}+\frac{s^8}{8064}+\frac{29s^7}{12096}+\frac{5s^6}{192}+\frac{3013s^5}{17280}+\frac{95s^4}{128}+\frac{4523s^3}{2268}+\frac{6515s^2}{2016}+\frac{7129s}{2520}+1$ \\ \hline

10 & 2 & $\frac{251s^9}{181440}+\frac{31s^8}{1008}+\frac{1765s^7}{6048}+\frac{37s^6}{24}+\frac{42863s^5}{8640}+\frac{481s^4}{48}+\frac{115205s^3}{9072}+\frac{4993s^2}{504}+\frac{5729s}{1260}+1$ \\ \hline

10 & 3 & $\frac{913s^9}{22680}+\frac{1135s^8}{2016}+\frac{5071s^7}{1512}+\frac{179s^6}{16}+\frac{3128s^5}{135} +\frac{2999s^4}{96}+\frac{63041s^3}{2268}+\frac{ 8069s^2}{504} +\frac{3553s}{630}+1$ \\ \hline

10 & 4 & $\frac{44117s^9}{181440}+\frac{2489s^8}{1008}+\frac{66547s^7}{6048}+\frac{683s^6}{24}+\frac{409361s^5}{ 8640}+\frac{2543s^4}{48}+\frac{363947s^3}{9072}+\frac{10127s^2}{504}+\frac{7883s}{1260}+1$ \\ \hline

10 & 5 & $ \frac{15619s^9}{36288}+\frac{15619s^8}{4032}+\frac{94939s^7}{6048}+\frac{3607s^6}{96}+\frac{101311s^5}{1728}+\frac{11911s^4}{192}+\frac{25394s^3}{567}+\frac{21689s^2}{1008}+\frac{1627s}{252}+1$ \\ \hline

11 & 1 & $\frac{s^{10}}{3628800}+\frac{11s^9}{725760}+\frac{11s^8}{30240}+\frac{121s^7}{24192}+\frac{7513s^6}{172800}+\frac{8591s^5}{34560}+\frac{341693s^4}{362880}+\frac{84095s^3}{36288}+\frac{177133s^2}{50400}+\frac{7381s}{2520}+1$ \\ \hline

11 & 2 & $\frac{1013s^{10}}{3628800}+\frac{ 5533s^9}{725760}+\frac{2189s^8}{24192}+\frac{14795s^7}{24192}+\frac{447689s^6}{172800}+\frac{246697s^5}{34560}+\frac{14597s^4}{1134}+\frac{543763s^3}{36288}+\frac{91949s^2}{8400}+\frac{1199s}{252}+1$ \\ \hline

11 & 3 & $\frac{299s^{10}}{22680}+\frac{16621s^9}{72576}+\frac{41591s^8}{24192}+\frac{88693s^7}{12096}+\frac{170137s^6}{8640}+\frac{604109s^5}{17280}+\frac{3043997s^4}{72576}+\frac{308473s^3}{9072}+\frac{60929s^2}{3360}+\frac{15059s}{2520}+1$ \\ \hline

11 & 4 & $\frac{56899s^{10}}{453600}+\frac{565631s^9}{362880}+\frac{205733s^8}{24192}+\frac{326491s^7}{12096}+\frac{2400629s^6}{43200}+\frac{1348787s^5}{17280}+\frac{5535695s^4}{72576}+\frac{468655s^3}{9072}+\frac{1185701s^2}{ 50400}+\frac{16973s}{2520}+1$ \\ \hline

11 & 5 & $\frac{655177s^{10}}{1814400}+\frac{336083s^9}{90720}+\frac{2078791s^8}{120960}+\frac{287639s^7}{6048}+\frac{7525771s^6}{86400}+\frac{95557s^5}{864}+\frac{35914087s^4}{362880}+\frac{1125575s^3}{18144}+\frac{443179s^2}{16800}+\frac{17897s}{2520}+1$ \\ \hline

12 & 1 & $\frac{s^{11}}{39916800}+\frac{s^{10}}{604800}+\frac{s^9}{20736}+\frac{11s^8}{13440}+\frac{10831s^7}{1209600}+\frac{1903s^6}{28800}+\frac{242537s^5}{725760}+\frac{139381s^4}{120960}+\frac{341747s^3}{129600}+\frac{190553s^2}{50400 }+\frac{83711s}{ 27720}+1 $ \\ \hline

12 & 2 & $\frac{509s^{11}}{9979200}+\frac{169s^{10}}{100800}+\frac{551s^9}{22680}+\frac{2057s^8}{10080}+\frac{332249s^7}{302400}+\frac{18997s^6}{4800}+\frac{876959s^5}{90720}+\frac{80179s^4}{5040 }+\frac{244681s^3}{14175}+\frac{150293s^2}{12600}+\frac{68591s}{13860}+1$ \\ \hline

12 & 3 & $\frac{50879s^{11}}{13305600}+\frac{6979s^{10}}{86400}+\frac{60271s^9}{80640}+\frac{32153s^8}{8064}+\frac{5483809s^7}{403200}+\frac{897259s^6}{28800}+\frac{11875111s^5}{241920}+\frac{185339s^4}{3456}+\frac{451173s^3}{11200}+\frac{338503s^2}{16800}+\frac{58007s}{9240}+1$ \\ \hline

12 & 4 & $\frac{1093s^{11}}{19800}+\frac{62879s^{10}}{75600}+\frac{20893s^9}{3780}+\frac{10813s^8}{504}+\frac{684323s^7}{12600}+\frac{340967s^6}{3600}+\frac{5258s^5}{45}+\frac{38819s^4}{378}+\frac{1202029s^3}{18900}+\frac{42218s^2}{1575}+\frac{1103s}{154}+1$ \\ \hline

12 & 5 & $\frac{1623019s^{11}}{6652800}+\frac{882773s^{10}}{302400}+\frac{1908073s^9}{120960}+\frac{1028401s^8}{20160}+\frac{7395023s^7}{67200}+\frac{2401619s^6}{14400}+\frac{4398559s^5}{24192}+\frac{8661917s^4}{60480}+\frac{12163441s^3}{151200}+\frac{782969s^2}{25200}+\frac{8861s}{1155}+1 $ \\ \hline

12 & 6 & $\frac{655177s^{11}}{1663200}+\frac{655177s^{10}}{151200}+\frac{5507s^9}{252}+\frac{336083s^8}{ 5040}+\frac{6898277s^7}{50400}+\frac{1430341s^6}{7200}+\frac{3152491s^5}{15120}+\frac{1200463s^4}{7560}+\frac{30291s^3}{350}+\frac{68321s^2}{2100}+\frac{18107s}{2310}+1$ \\ \hline

\end{tabular}
\caption{The Ehrhart polynomials for the hypersimplices $\Delta(n, k)$ } 
\label{hypersi}

}

\end{sidewaystable}


We also have the Ehrhart quasi-polynomials of some truncated unit cubes.

\begin{proposition}
The Ehrhart quasi-polynomial for the truncated unit cube in Figure 
\ref{tru1}, where its vertices are at $1/3$ and $2/3$ of the way along
edges of the cube, is given by:

$$i_{tru\_cube1}(s) = \left \{ \begin{array}{ll} \frac{77s^3}{81} + 
\frac{23s^2}{9} + \frac{19s}{9} + 1 \mbox{ if } s \equiv 0 \mod 3, \\ 
 \frac{77s^3}{81} + \frac{61s^2}{27} - \frac{7s}{27}
- \frac{239}{81} \mbox{ if } s \equiv 1 \mod 3, \\
 \frac{77s^3}{81} + \frac{65s^2}{27} +
\frac{29s}{27} - \frac{31}{81} \mbox{ if } s \equiv 2 \mod 3.\\
\end{array}\right .$$


\begin{figure}[ht]
\begin{center}
\includegraphics[width=5 cm]{tru_cube1.eps}  
\caption{The truncated cube.} \label{tru1}
\end{center} 
\end{figure}
\end{proposition}

\begin{proposition}
The Ehrhart quasi-polynomial for the cuboctahedron (Figure \ref{tru2}) is:

$$i_{tru\_cube2}(s) = \left \{ \begin{array}{ll} \frac{5s^3}{6} + 2s^2 + 
\frac{5s}{3} + 1 \mbox{ if } s\equiv 0 \mod 2, \\
 \frac{5s^3}{6} + \frac{3s^2}{2} - \frac{5s}{6} - \frac{3}{2} \mbox{ if } 
s \equiv 1 \mod 2.\\ 
\end{array}\right .$$


\begin{figure}[ht]
\begin{center}
\includegraphics[width=5 cm]{tru_cube2.eps} 
\caption{The cuboctahedron.} \label{tru2}
\end{center} 
\end{figure}

\end{proposition}

\begin{proposition}
The Ehrhart quasi-polynomial for the truncated regular simplex, where
the vertices are at $1/3$ and $2/3$ of the way along the simplex edges
(see Figure \ref{tru3}), is given by:

$$i_{tru\_simplex}(s) = \left \{ \begin{array}{ll}
\frac{23s^3}{81} + \frac{7s^2}{9} + \frac{13s}{9}
+ 1 \mbox{ if } s \equiv 0 \mod 3, \\
 \frac{23s^3}{81} + \frac{19s^2}{27} +
\frac{5s}{27} - \frac{95}{81} \mbox{ if } s \equiv 1 \mod 3, \\
 \frac{23s^3}{81} + \frac{17s^2}{27} +
\frac{23s}{27} + \frac{41}{81} \mbox{ if } s \equiv 2 \mod 3.\\
\end{array}\right .$$


\begin{figure}[ht]
\begin{center}
\includegraphics[width=5 cm]{tru_simplex.eps} 
\caption{The truncated simplex.} \label{tru3}
\end{center} 
\end{figure}
\end{proposition}


\section{Computational results via Homogenized Barvinok's algorithm} \label{theory}

We have demonstrated the practical relevance of Barvinok's cone
decomposition approach for counting lattice points and deriving
formulas. Several other algorithms are available to carry out the same
kind of enumeration. It is important to implement them all in the same
computer system for comparison of performance and to corroborate that
the answers are correct. Some problems are solvable by some
methods but not by others.

We would like to remind a reader that one of the bottle necks of the 
original Barvinok algorithm in \citep{bar} is the fact that a polytope
may have too many vertices.  Originally, we visit all vertices of the input
polytope to compute Barvinok's short rational function and this can be
costly in terms of computational time.  For example, the well-known polytope of semi-magic cubes in
the $ 4\times 4 \times 4$ case has over two million vertices, but only
64 linear inequalities describe the polytope. 
Algorithm \ref{homogbarvinok} shows that Homogenized Barvinok's algorithm
works with only a single cone.  In this section, we will show some
practical results with Homogenized Barvinok's algorithm.


 


A \emph{normal semigroup} $S$ is the intersection of the lattice
$\jdlZ^d$ with a rational convex polyhedral cone in $\R^d$. 
Each pointed affine semigroup $S \subset \Z^d$ can be {\em graded}. This means that there is a linear map $\,deg: S \rightarrow \N$ with
$deg(x)=0$ if and only if $x=0$.  Given a pointed graded affine
semigroup, we define $S_r$ to be the set of elements with degree $r$, i.e.
$S_r=\{x \in S : deg(x)=r\}$.  The {\em Hilbert series} of $S$ is the
formal power series $H_S(t)=\sum_{k=0}^\infty \# (S_r)
t^r$, where $\# (S_r)$ is the cardinality of $S_r$. Algebraically, this is just the Hilbert series of the semigroup
ring $\C[S]$. It is a well-known property that $H_S$ is
represented by a rational function of the form

$$\frac{Q(t)}{(1-t^{s_1})(1-t^{s_2})\dots (1-t^{s_d})}$$

where $Q(t)$ is a polynomial of degree less than $s_1+\dots+s_d$ (see
Chapter 4 \citep{stanley}).  
The first challenge is to compute the Hilbert series of magic cubes.
Several other methods had been tried to
compute the Hilbert series explicitly (see \citep{ahmed} for
references).  One of the most well-known challenges was that of
counting the $5 \times 5$ magic squares of magic sum $n$. Similarly
several authors had tried before to compute the Hilbert series of the
$3 \times 3 \times 3 \times 3$ magic cubes. It is not difficult
to see this is equivalent to determining an Ehrhart series. Using
Algorithm \ref{homogbarvinok} we finally present the solution, which had been
inaccessible using Gr\"obner bases methods. For comparison, the reader
familiar with Gr\"obner bases computations should be aware that the
$5\times 5$ magic squares problem required a computation of a
Gr\"obner bases of a toric ideal of a matrix $A$ with 25 rows and over
4828 columns. Our attempts to handle this problem with {\tt CoCoA} and
{\tt Macaulay2} were unsuccessful.  We now give the numerator and then
the denominator of the rational functions computed with the software
{\tt LattE}:

\begin{theorem} 
\label{magic}

\item The generating function $\,\sum_{n \geq 0} f(n) t^n \,$ for the number $f(n)$ of $5\times
5$ magic squares of magic sum $n$ is given by the rational function
$p(t)/q(t)$ with numerator

{\footnotesize

$p(t)=
{t}^{76}+28\,{t}^{75}+639\,{t}^{74}+11050\,{t}^{73}+136266\,{t}^{72}+
1255833\,{t}^{71}+9120009\,{t}^{70}+54389347\,{t}^{69}+\newline
274778754\,{t}^{68}+1204206107\,{t}^{67}+4663304831\,{t}^{66}+
16193751710\,{t}^{65}+51030919095\,{t}^{64}+ \newline
147368813970\,{t}^{63}+393197605792\,{t}^{62}+
975980866856\,{t}^{61}+
2266977091533\,{t}^{60}+
4952467350549\,{t}^{59}+\newline
10220353765317\,{t}^{58}+
20000425620982\,{t}^{57}+37238997469701\,{t}^{56}+
66164771134709\,{t}^{55}+112476891429452\,{t}^{54}+\newline
183365550921732\,{t}^{53}+
287269293973236\,{t}^{52}+433289919534912\,{t}^{51}+
630230390692834\,{t}^{50}\newline +885291593024017\,{t}^{49}+
1202550133880678\,{t}^{48}+1581424159799051\,{t}^{47}+
2015395674628040\,{t}^{46}+\newline 2491275358809867\,{t}^{45}+
2989255690350053\,{t}^{44}+3483898479782320\,{t}^{43}+
3946056312532923\,{t}^{42}+\newline
4345559454316341\,{t}^{41}+
4654344257066635\,{t}^{40}+4849590327731195\,{t}^{39}+
4916398325176454\,{t}^{38}+\newline 4849590327731195\,{t}^{37}+
4654344257066635\,{t}^{36}+
4345559454316341\,{t}^{35}+
3946056312532923\,{t}^{34}+\newline 3483898479782320\,{t}^{33}+
2989255690350053\,{t}^{32}+2491275358809867\,{t}^{31}+
2015395674628040\,{t}^{30}+\newline
1581424159799051\,{t}^{29}+
1202550133880678\,{t}^{28}+885291593024017\,{t}^{27}+
630230390692834\,{t}^{26}+433289919534912\,{t}^{25}+
287269293973236\,{t}^{24}+183365550921732\,{t}^{23}+
112476891429452\,{t}^{22}+66164771134709\,{t}^{21}+37238997469701
\,{t}^{20}+\newline  20000425620982\,{t}^{19}+
10220353765317\,{t}^{18}+
4952467350549\,{t}^{17}+2266977091533\,{t}^{16}+975980866856\,{t}^
{15}+\newline
393197605792\,{t}^{14}+147368813970\,{t}^{13}+51030919095\,{t
}^{12}+16193751710\,{t}^{11}+4663304831\,{t}^{10}+1204206107\,{t}^
{9}+\newline 274778754\,{t}^{8}+54389347\,{t}^{7}+9120009\,{t}^{6}+1255833
\,{t}^{5}+136266\,{t}^{4}+11050\,{t}^{3}+639\,{t}^{2}+28\,t+1
$ \ \ \ \ \ and denominator

$ q(t)=
\left ({t}^{2}-1\right )^{10}\left ({t}^{2}+t+1\right )
^{7}\left ({t}^{7}-1\right )^{2}
\left ({t}^{6}+{t}^{3}+1\right )\left ({t}^{5}+{t}^{3}+{t}^{2}+t+1
\right )^{4}\left (1-t\right )^{3}\left ({t
}^{2}+1\right )^{4}
$.

}

The generating function $\,\sum_{n \geq 0} f(n) t^n \,$ for the number
$f(n)$ of $3\times 3 \times 3 \times 3$ magic cubes with magic sum $n$
is given the rational function $r(t)/s(t)$ where

{\footnotesize
$
r(t)={t}^{54}+150\,{t}^{51}+5837\,{t}^{48}+63127\,{t}^{45}+331124\,{t}^{42}
+1056374\,{t}^{39}+2326380\,{t}^{36}+3842273\,{t}^{33}+5055138\,{t}^{
30}+5512456\,{t}^{27}+5055138\,{t}^{24}+3842273\,{t}^{21}+2326380\,{t}
^{18}+1056374\,{t}^{15}+331124\,{t}^{12}+63127\,{t}^{9}+5837\,{t}^{6}+
150\,{t}^{3}+1
$ \ \ \ \ \ and

$
s(t)=\left ({t}^{3}+1\right )^{4}\left ({t}^{
12}+{t}^{9}+{t}^{6}+{t}^{3}+1\right )\left (1-{t}^{3}\right )^{9}
\left ({t}^{6}+{t}^{3}+1\right )
$.
}
\end{theorem}


\input{IPThesis}
{"config": "arxiv", "file": "math0406284/softwareThesis.tex"}