short_rat::print: correctly print out terms with a zero coefficient

[barvinok.git] / doc / implementation.tex

\section{Implementation details}

\subsection{An interior point of a polyhedron}
\label{s:interior}

We often need a point that lies in the interior of a polyhedron.
The function \ai[\tt]{inner\_point} implements the following algorithm.
Each polyhedron $P$ can be written as the sum of a polytope $P'$ and a cone $C$
(the \ai{recession cone} or \ai{characteristic cone} of $P$).
Adding a positive multiple of the sum of the extremal rays of $C$ to
the \ai{barycenter}
$$
\frac 1 N \sum_i \vec v_i(\vec p)
$$
of $P'$, where $N$ is the number of vertices, results in a point
in the interior of $P$.

\subsection{The integer points in the fundamental parallelepiped of a simple cone}

\label{s:fundamental}

This section is based on \shortciteN[Lemma 5.1]{Barvinok1992volume} and
\shortciteN{Koeppe2006experiments}.

\sindex{simple}{cone}
\sindex{open}{facet}
\sindex{open}{ray}
\sindex{explicit}{representation}
In this section we will deal exclusively with \ai{simple cone}s,
i.e. $d$-dimensional cones with $d$ extremal rays and $d$ facets.
\index{open facet}%
Some of the facets of these cones may be open.
Since we will mostly be dealing with cones in their
\ai{explicit representation}, we will have occasion to speak of
``\ai{open ray}s'', by which we will mean that the facet not
containing the ray is open.  (There is only one such facet because the cone
is simple.)

\sindex{fundamental}{parallelepiped}
\begin{definition}[Fundamental parallelepiped]
Let $K = \vec v + \poshull \lb\, \vec u_i \,\rb$ be 
a closed (shifted) cone, then the \defindex{fundamental parallelepiped} $\Pi$
of $K$ is
$$
\Pi = \vec v +
\lb\, \sum_i \alpha_i \vec u_i \mid 0 \leq \alpha_i < 1 \,\rb
.
$$
If some of the rays $\vec u_i$ of $K$ are open, then the constraints on
the corresponding coefficient $\alpha_i$ are such that $0 < \alpha_i \le 1$.
\end{definition}

\begin{lemma}[Integer points in the fundamental parallelepiped of a simple cone]
\label{l:fundamental}
Let $K = \vec v + \poshull \lb\, \vec u_i \,\rb$ be a closed simple cone
and let $A$ be the matrix with the generators $\vec u_i$ of $K$
as rows.
Furthermore let $V A W^{-1} = S = \diag \vec s$ be the \indac{SNF} of $A$.
Then the integer points in the fundamental parallelepiped of $K$ are given
by
\begin{eqnarray}
\label{eq:parallelepiped}
\vec w^\T & = & \vec v^\T + \fractional{(\vec k^\T W - \vec v^\T) A^{-1}} A
\\
\nonumber
& = &
\vec v^\T +
\sum_{i=1}^d
    \fractional{\sps{\sum_{j=1}^d k_j \vec w^\T_j - \vec v^\T}{\vec u^*_i}}
        \vec u_i^\T,
\end{eqnarray}
where $\vec u^*_i$ are the columns of $A^{-1}$ and $k_j \in \ZZ$ ranges
over $0 \le k_j < s_j$.
\end{lemma}

\begin{proof}
Since $0 \le \fractional{x} < 1$, it is clear that each such $\vec w$
lies inside the fundamental parallelepiped.
Furthermore,
\begin{eqnarray*}
\vec w^\T & = & \vec v^\T + \fractional{(\vec k^\T W - \vec v^\T) A^{-1}} A
\\
& = &
\vec v^T + 
\left(
(\vec k^\T W - \vec v^\T) A^{-1} - \floor{(\vec k^\T W - \vec v^\T) A^{-1}}
\right) A
\\
& = &
\underbrace{\vec k^\T W\mathstrut}_{\in \ZZ^{1\times d}}
-
\underbrace{\floor{(\vec k^\T W - \vec v^\T) A^{-1}}}_{\in \ZZ^{1\times d}}
\underbrace{A\mathstrut}_{\in \ZZ^{d\times d}} \in \ZZ^{1\times d}.
\end{eqnarray*}
Finally, if two such $\vec w$ are equal, i.e., $\vec w_1 = \vec w_2$,
then
\begin{eqnarray*}
\vec 0^\T = \vec w_1^\T - \vec w_2^\T
& = & \vec k_1^\T W - \vec k_2^\T W + \vec p^\T A
\\
& = & \left(\vec k_1^\T - \vec k_2^\T \right) W + \vec p^\T V^{-1} S W,
\end{eqnarray*}
with $\vec p \in \ZZ^d$,
or $\vec k_1 \equiv \vec k_2 \mod \vec s$, i.e., $\vec k_1 = \vec k_2$.
Since $\det S = \det A$, we obtain all points in the fundamental parallelepiped
by taking all $\vec k \in \ZZ^d$ satisfying $0 \le k_j < s_j$.
\end{proof}

If the cone $K$ is not closed then the coefficients of the open rays
should be in $(0,1]$ rather than in $[0,1)$.
In (\ref{eq:parallelepiped}),
we therefore need to replace the fractional part $\fractional{x} = x - \floor{x}$
by $\cractional{x} = x - \ceil{x-1}$ for the open rays.

\begin{figure}
\intercol=1.2cm
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS@i@={(0,-3),(0,0),(4,2),(4,-3)},{0*[grey]\xypolyline{*}}
\POS@i@={(0,-3),(0,0),(4,2)},{0*[|(2)]\xypolyline{}}
\POS(-1,0)\ar(4.5,0)
\POS(0,-3)\ar(0,3)
\POS(0,0)\ar@[|(3)](0,-1)
\POS(0,0)\ar@[|(3)](2,1)
\POS(0,-1)\ar@{--}@[|(2)](2,0)
\POS(2,1)\ar@{--}@[|(2)](2,0)
\POS(0,0)*{\bullet}
\POS(1,0)*{\bullet}
\end{xy}
\caption{The integer points in the fundamental parallelepiped of $K$}
\label{f:parallelepiped}
\end{figure}

\begin{example}
Let $K$ be the cone
$$
K = \sm{0 \\ 0} + \poshull \lb\, \sm{2 \\ 1}, \sm{0 \\ -1} \,\rb
,
$$
shown in Figure~\ref{f:parallelepiped}.
Then
$$
A = \sm{2 & 1\\0 & -1} \qquad A^{-1} = \sm{1/2 & 1/2 \\ 0 & -1 }
$$
and
$$
\sm{1 & 0 \\ 1 & 1 } \sm{2 & 1\\0 & -1} = \sm{1 & 0 \\ 0 & 2} \sm{2 & 1 \\ 1 & 0}.
$$
We have $\det A = \det S = 2$ and
$\vec k_1^\T = \sm{0 & 0}$ and $\vec k_2^\T = \sm{0 & 1}$.
Therefore,
$$
\vec w_1^\T = \fractional{\sm{0 & 0} \sm{2 & 1 \\ 1 & 0} \sm{1/2 & 1/2 \\ 0 & -1 }}
\sm{2 & 1\\0 & -1} = \sm{0 & 0}
$$
and
\begin{eqnarray*}
\vec w_2^\T & = & 
\fractional{\sm{0 & 1} \sm{2 & 1 \\ 1 & 0} \sm{1/2 & 1/2 \\ 0 & -1 }}
\sm{2 & 1\\0 & -1}
\\
& = &
\sm{1/2 & 1/2} \sm{2 & 1\\0 & -1} = \sm{1 & 0}.
\end{eqnarray*}
\end{example}




\subsection{Barvinok's decomposition of simple cones in primal space}
\label{s:decomposition}

As described by \shortciteN{DeLoera2003effective}, the first
implementation of Barvinok's counting algorithm applied
\ai{Barvinok's decomposition} \shortcite{Barvinok1994} in the \ai{dual space}.
\ai{Brion's polarization trick} \shortcite{Brion88} then ensures that you
do not need to worry about lower-dimensional faces in the decomposition.
Another way of avoiding the lower-dimensional faces, in the \ai{primal space},
is to perturb the vertex of the cone such that none of the lower-dimensional
face encountered contain any integer points \shortcite{Koeppe2006primal}.
In this section, we describe another technique that is based on allowing
some of the facets of the cone to be open.

The basic step in Barvinok's decomposition is to replace a
$d$-dimensional simple cone 
$K = \poshull \lb\, \vec u_i \,\rb_{i=1}^d \subset \QQ^d$
by a signed sum of (at most) $d$ cones $K_j$
with a smaller determinant (in absolute value).
The cones are obtained by successively replacing each generator
of $K$ by an appropriately chosen
$\vec w = \sum_{i=1}^d \alpha_i \vec u_i$, i.e.,
\begin{equation}
\label{eq:K_j}
K_j = 
\poshull \left(\lb\, \vec u_i \,\rb_{i=1}^d 
\setminus \{\, \vec u_j \,\} \cup \{\, \vec w \,\}\right)
.
\end{equation}
To see that we can use these $K_j$ to perform a decomposition,
rearrange the $\vec u_i$ such that for all $1 \le i \le k$ we have
$\alpha_i < 0$ and for all $k+1 \le i \le d'$ we have $\alpha_i > 0$,
with $d - d'$ the number of zero $\alpha_i$.
We may assume $k < d'$; otherwise replace $\vec w \in B$ by
$-\vec w \in B$.  We have
$$
\vec w + \sum_{i=1}^k (-\alpha_i) \vec u_i =
\sum_{i=k+1}^{d'} \alpha_i \vec u_i
$$
or
\begin{equation}
\label{eq:sub}
\sum_{i=0}^k \beta_i \vec u_i =
\sum_{i=k+1}^{d'} \alpha_i \vec u_i
,
\end{equation}
with $\vec u_0 = \vec w$, $\beta_0 = 1$ and $\beta_i = -\alpha_i > 0$
for $1 \le i \le k$.  Any two $\vec u_j$ and $\vec u_l$ on the same side
of the equality are on opposite sides of the linear hull $H$ of
the other $\vec u_i$s since there exists a convex combination
of $\vec u_j$ and $\vec u_l$ on this hyperplane.
In particular, since $\alpha_j$ and $\alpha_l$ have the same sign,
we have
\begin{equation}
\label{eq:opposite}
\frac {\alpha_j}{\alpha_j+\alpha_l} \vec u_j
+
\frac {\alpha_l}{\alpha_j+\alpha_l} \vec u_l
\in H
\qquad\text{for $\alpha_i \alpha_l > 0$}
.
\end{equation}
The corresponding cones $K_j$ and $K_l$ (with $K_0 = K$)
therefore intersect in a common face $F \subset H$.
Let 
$$
K' := 
\poshull \left(\lb\, \vec u_i \,\rb_{i=1}^d \cup \{\, \vec w \,\}\right)
,
$$
then any $\vec x \in K'$ lies both in some cone $K_i$ with
$0 \le i \le k$ and in some cone $K_i$ with $k+1 \le i \le d'$.
(Just subtract an appropriate multiple of Equation~(\ref{eq:sub}).)
The cones 
$\{\, K_i \,\}_{i=0}^k$
and
$\{\, K_i \,\}_{i=k+1}^{d'}$
therefore both form a triangulation of $K'$ and hence
\begin{equation}
\label{eq:triangulations}
\indf{K'}
=
\indf{K} + \sum_{i=1}^k \indf{K_i} - \sum_{j\in J_1} \indf{F_j}
=
\sum_{i=k+1}^{d'} \indf{K_i} - \sum_{j\in J_2} \indf{F_j}
\end{equation}
or
\begin{equation}
\label{eq:decomposition}
\indf{K} = \sum_{i=1}^{d'} \varepsilon_i \indf{K_i} + \sum_j \delta_j \indf{F_j}
,
\end{equation}
with $\varepsilon_i = -1$ for $1 \le i \le k$,
$\varepsilon_i = 1$ for $k+1 \le i \le d'$,
$\delta_j \in \{ -1, 1 \}$ and $F_j$ some lower-dimensional faces.
Figure~\ref{fig:w} shows the possible configurations
in the case of a $3$-dimensional cone.

\begin{figure}
\intercol=0.48cm
\begin{center}
\begin{minipage}{0cm}
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
*
\xybox{
\POS(-2,-1)="a"*+!U{+}
\POS(2,0)="b"*+!L{+}
\POS(0,2)="c"*+!D{+}
\POS(0,0)="w"*+!DR{\vec w}
\POS"a"\ar@{-}"b"
\POS"b"\ar@{-}"c"
\POS"c"\ar@{-}"a"
\POS"a"\ar@{--}"w"
\POS"b"\ar@{--}"w"
\POS"c"\ar@{--}"w"
}="a"
+R+(2,0)*!L
\xybox{
\POS(-2,-1)="a"*+!U{+}
\POS(2,0)="b"*+!L{-}
\POS(0,2)="c"*+!D{+}
\POS(-3,1)="w"*+!DR{\vec w}
\POS"a"\ar@{-}"b"
\POS"b"\ar@{-}"c"
\POS"c"\ar@{-}"a"
\POS"a"\ar@{--}"w"
\POS"b"\ar@{--}"w"
\POS"c"\ar@{--}"w"
}="b"
+R+(2,0)*!L
\xybox{
\POS(-2,-1)="a"*+!U{-}
\POS(2,0)="b"*+!U{+}
\POS(0,2)="c"*+!D{-}
\POS(5,-1)="w"*+!L{\vec w}
\POS"a"\ar@{-}"b"
\POS"b"\ar@{-}"c"
\POS"c"\ar@{-}"a"
\POS"a"\ar@{--}"w"
\POS"b"\ar@{--}"w"
\POS"c"\ar@{--}"w"
}
\POS"a"
+D-(0,1)*!U
\xybox{
\POS(-2,-1)="a"*+!U{0}
\POS(2,0)="b"*+!L{+}
\POS(0,2)="c"*+!D{+}
\POS(1,1)="w"*+!DL{\vec w}
\POS"a"\ar@{-}"b"
\POS"b"\ar@{-}"c"
\POS"c"\ar@{-}"a"
\POS"a"\ar@{--}"w"
}
\POS"b"
+DL-(0,1)*!UL
\xybox{
\POS(-2,-1)="a"*+!U{0}
\POS(2,0)="b"*+!U{+}
\POS(0,2)="c"*+!D{-}
\POS(4,-2)="w"*+!L{\vec w}
\POS"a"\ar@{-}"b"
\POS"b"\ar@{-}"c"
\POS"c"\ar@{-}"a"
\POS"a"\ar@{--}"w"
\POS"b"\ar@{--}"w"
}
\end{xy}
\end{minipage}
\end{center}
\caption[Possible locations of the vector $\vec w$ with respect to the rays
of a $3$-dimensional cone.]
{Possible locations of $\vec w$ with respect to the rays
of a $3$-dimensional cone.  The figure shows a section of the cones.}
\label{fig:w}
\end{figure}

As explained above there are several ways of avoiding the lower-dimensional
faces in (\ref{eq:decomposition}).  Here we will apply the following proposition.
\begin{proposition}[\shortciteN{Koeppe2007parametric}]
\label{p:inclusion-exclusion}
  Let 
  \begin{equation}
    \label{eq:full-source-identity}
    \sum_{i\in {I_1}} \epsilon_i [P_i] + \sum_{i\in {I_2}} \delta_k [P_i] = 0
  \end{equation}
  be a (finite) linear identity of indicator functions of closed
  polyhedra~$P_i\subseteq\QQ^d$, where the
  polyhedra~$P_i$ with $i \in I_1$ are full-dimensional and those with $i \in I_2$
  lower-dimensional.  Let each closed polyhedron be given as 
$$
    P_i = \left\{\, \vec x \mid \sp{b^*_{i,j}}{x} \ge \beta_{i,j} \text{
      for $j\in J_i$}\,\right\}
  .
$$
  Let $\vec y\in\QQ^d$ be a vector such that $\langle \vec b^*_{i,j}, \vec
  y\rangle \neq 0$ for all $i\in I_1\cup I_2$, $j\in J_i$.
  For each $i\in I_1$, we define the half-open polyhedron
  \begin{equation}
    \label{eq:half-open-by-y}
    \begin{aligned}
      \tilde P_i = \Bigl\{\, \vec x\in\QQ^d \mid {}&
            \sp{b^*_{i,j}}{x} \ge \beta_{i,j}
      \text{ for $j\in J_i$ with $\sp{b^*_{i,j}}{y} > 0$,} \\
      & \sp{b^*_{i,j}}{x} > \beta_{i,j}
      \text{ for $j\in J_i$ with $\sp{b^*_{i,j}}{y} < 0$} \,\Bigr\}.
    \end{aligned}
  \end{equation}
  Then 
  \begin{equation}
    \label{eq:target-identity}
    \sum_{i\in I_1} \epsilon_i [\tilde P_i] = 0.
  \end{equation}
\end{proposition}
When applying this proposition to (\ref{eq:decomposition}), we obtain
\begin{equation}
\label{eq:decomposition:2}
\indf{\tilde K} = \sum_{i=1}^{d'} \varepsilon_i \indf{\tilde K_i}
,
\end{equation}
where we start out
from a given $\tilde K$, which may be $K$ itself, i.e., a fully closed cone,
or the result of a previous application of the proposition, either through
a triangulation (Section~\ref{s:triangulation}) or a previous decomposition.
In either case, a suitable $\vec y$ is available, either as an interior
point of the cone or as the vector used in the previous application
(which may require a slight perturbation if it happens to lie on one of
the new facets of the cones $K_i$).
We are, however, free to construct a new $\vec y$ on each application
of the proposition.
In fact, we will not even construct such a vector explicitly, but
rather apply a set of rules that is equivalent to a valid choice of $\vec y$.
Below, we will present an ``intuitive'' motivation for these rules.
For a more algebraic, shorter, and arguably simpler motivation we
refer to \shortciteN{Koeppe2007parametric}.

The vector $\vec y$ has to satisfy $\sp{b^*_j}y > 0$ for normals $\vec b^*_j$
of closed facets and $\sp{b^*_j}y < 0$ for normals $\vec b^*_j$ of open facets of
$\tilde K$.
These constraints delineate a non-empty open cone $R$ from which
$\vec y$ should be selected.  For some of the new facets of the cones
$\tilde K_j$, the cone $R$ will not be cut by the affine hull of the facet.
The closedness of these facets is therefore predetermined by $\tilde K$.
For the other facets, a choice will have to be made.
To be able to make the choice based on local information and without
computing an explicit vector $\vec y$, we use the following convention.
We first assign an arbitrary total order to the rays.
If (the affine hull of) a facet separates the two rays not on the facet $\vec u_i$
and $\vec u_j$, i.e., $\alpha_i \alpha_j > 0$ (\ref{eq:opposite}), then
we choose $\vec y$ to lie on the side of the smallest ray, according
to the chosen order.
That is, $\sp{{\tilde n}_{ij}}y > 0$, for
$\vec {\tilde n}_{ij}$ the normal of the facet pointing towards this smallest ray.
Otherwise, i.e., if $\alpha_i \alpha_j < 0$,
the interior of $K$ will lie on one side
of the facet and then we choose $\vec y$ to lie on the other side.
That is, $\sp{{\tilde n}_{ij}}y > 0$, for
$\vec {\tilde n}_{ij}$ the normal of the facet pointing away from the cone $K$.
Figure~\ref{fig:primal:examples} shows some example decompositions with
an explicitly marked $\vec y$.

\begin{figure}
\begin{align*}
\intercol=0.48cm
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-2,-1)="a"*+!U{+}
\POS(2,0)="b"*+!L{+}
\POS(0,2)="c"*+!D{+}
\POS"a"\ar@{-}@[|(3)]"b"
\POS"b"\ar@{-}@[|(3)]"c"
\POS"c"\ar@{-}@[|(3)]"a"
\POS(0.3,0.6)*{\bullet},*+!L{\vec y}
\end{xy}
& =
\intercol=0.48cm
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(2,0)="b"*+!L{+}
\POS(0,2)="c"*+!D{+}
\POS(0,0)="w"*+!DR{\vec w}
\POS"b"\ar@{-}@[|(3)]"c"
\POS"b"\ar@{-}@[|(3)]"w"
\POS"c"\ar@{-}@[|(3)]"w"
\POS(0.3,0.6)*{\bullet},*+!L{\vec y}
\end{xy}
+
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-2,-1)="a"*+!U{+}
\POS(0,2)="c"*+!D{+}
\POS(0,0)="w"*+!DR{\vec w}
\POS"c"\ar@{-}@[|(3)]"a"
\POS"a"\ar@{-}@[|(3)]"w"
\POS"c"\ar@{--}@[|(3)]"w"
\POS(0.3,0.6)*{\bullet},*+!L{\vec y}
\end{xy}
+
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-2,-1)="a"*+!U{+}
\POS(2,0)="b"*+!L{+}
\POS(0,0)="w"*+!DR{\vec w}
\POS"a"\ar@{-}@[|(3)]"b"
\POS"a"\ar@{--}@[|(3)]"w"
\POS"b"\ar@{--}@[|(3)]"w"
\POS(0.3,0.6)*{\bullet},*+!L{\vec y}
\end{xy}
\\
\intercol=0.48cm
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-2,-1)="a"*+!U{+}
\POS(2,0)="b"*+!L{+}
\POS(0,2)="c"*+!D{+}
\POS"a"\ar@{--}@[|(3)]"b"
\POS"b"\ar@{-}@[|(3)]"c"
\POS"c"\ar@{--}@[|(3)]"a"
\POS(-2.5,-1.5)*{\bullet},*+!U{\vec y}
\end{xy}
& =
\intercol=0.48cm
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(2,0)="b"*+!L{+}
\POS(0,2)="c"*+!D{+}
\POS(0,0)="w"*+!DR{\vec w}
\POS"b"\ar@{-}@[|(3)]"c"
\POS"b"\ar@{--}@[|(3)]"w"
\POS"c"\ar@{--}@[|(3)]"w"
\POS(-2.5,-1.5)*{\bullet},*+!U{\vec y}
\end{xy}
+
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-2,-1)="a"*+!U{+}
\POS(0,2)="c"*+!D{+}
\POS(0,0)="w"*+!DR{\vec w}
\POS"c"\ar@{--}@[|(3)]"a"
\POS"a"\ar@{--}@[|(3)]"w"
\POS"c"\ar@{-}@[|(3)]"w"
\POS(-2.5,-1.5)*{\bullet},*+!U{\vec y}
\end{xy}
+
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-2,-1)="a"*+!U{+}
\POS(2,0)="b"*+!L{+}
\POS(0,0)="w"*+!DR{\vec w}
\POS"a"\ar@{--}@[|(3)]"b"
\POS"a"\ar@{-}@[|(3)]"w"
\POS"b"\ar@{-}@[|(3)]"w"
\POS(-2.5,-1.5)*{\bullet},*+!U{\vec y}
\end{xy}
\\
\intercol=0.48cm
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-2,-1)="a"*+!U{+}
\POS(2,0)="b"*+!L{+}
\POS(0,2)="c"*+!D{+}
\POS"a"\ar@{--}@[|(3)]"b"
\POS"b"\ar@{-}@[|(3)]"c"
\POS"c"\ar@{-}@[|(3)]"a"
\POS(1,-1.5)*{\bullet},*+!L{\vec y}
\end{xy}
& =
\intercol=0.48cm
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(2,0)="b"*+!L{+}
\POS(0,2)="c"*+!D{+}
\POS(0,0)="w"*+!DR{\vec w}
\POS"b"\ar@{-}@[|(3)]"c"
\POS"b"\ar@{--}@[|(3)]"w"
\POS"c"\ar@{-}@[|(3)]"w"
\POS(1,-1.5)*{\bullet},*+!L{\vec y}
\end{xy}
+
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-2,-1)="a"*+!U{+}
\POS(0,2)="c"*+!D{+}
\POS(0,0)="w"*+!DR{\vec w}
\POS"c"\ar@{-}@[|(3)]"a"
\POS"a"\ar@{--}@[|(3)]"w"
\POS"c"\ar@{--}@[|(3)]"w"
\POS(1,-1.5)*{\bullet},*+!L{\vec y}
\end{xy}
+
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-2,-1)="a"*+!U{+}
\POS(2,0)="b"*+!L{+}
\POS(0,0)="w"*+!DR{\vec w}
\POS"a"\ar@{--}@[|(3)]"b"
\POS"a"\ar@{-}@[|(3)]"w"
\POS"b"\ar@{-}@[|(3)]"w"
\POS(1,-1.5)*{\bullet},*+!L{\vec y}
\end{xy}
\\
\intercol=0.48cm
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-2,-1)="a"*+!U{0}
\POS(2,0)="b"*+!L{+}
\POS(0,2)="c"*+!D{+}
\POS"a"\ar@{-}@[|(3)]"b"
\POS"b"\ar@{-}@[|(3)]"c"
\POS"c"\ar@{-}@[|(3)]"a"
\POS(1,0.2)*{\bullet},*+!R{\vec y}
\end{xy}
& =
\intercol=0.48cm
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-2,-1)="a"*+!U{0}
\POS(2,0)="b"*+!L{+}
\POS(1,1)="w"*+!DL{\vec w}
\POS"a"\ar@{-}@[|(3)]"b"
\POS"a"\ar@{-}@[|(3)]"w"
\POS"b"\ar@{-}@[|(3)]"w"
\POS(1,0.2)*{\bullet},*+!R{\vec y}
\end{xy}
+
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-2,-1)="a"*+!U{0}
\POS(0,2)="c"*+!D{+}
\POS(1,1)="w"*+!DL{\vec w}
\POS"c"\ar@{-}@[|(3)]"a"
\POS"a"\ar@{--}@[|(3)]"w"
\POS"c"\ar@{-}@[|(3)]"w"
\POS(1,0.2)*{\bullet},*+!R{\vec y}
\end{xy}
\\
\intercol=0.48cm
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-2,-1)="a"*+!U{0}
\POS(2,0)="b"*+!U{+}
\POS(0,2)="c"*+!D{-}
\POS"a"\ar@{-}@[|(3)]"b"
\POS"b"\ar@{--}@[|(3)]"c"
\POS"c"\ar@{-}@[|(3)]"a"
\POS(1.5,1.5)*{\bullet},*+!D{\vec y}
\end{xy}
& =
\intercol=0.48cm
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-2,-1)="a"*+!U{0}
\POS(0,2)="c"*+!D{-}
\POS(4,-2)="w"*+!L{\vec w}
\POS"c"\ar@{-}@[|(3)]"a"
\POS"a"\ar@{-}@[|(3)]"w"
\POS"c"\ar@{--}@[|(3)]"w"
\POS(1.5,1.5)*{\bullet},*+!D{\vec y}
\end{xy}
-
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-2,-1)="a"*+!U{0}
\POS(2,0)="b"*+!U{+}
\POS(4,-2)="w"*+!L{\vec w}
\POS"a"\ar@{--}@[|(3)]"b"
\POS"a"\ar@{-}@[|(3)]"w"
\POS"b"\ar@{--}@[|(3)]"w"
\POS(1.5,1.5)*{\bullet},*+!D{\vec y}
\end{xy}
\\
\intercol=0.48cm
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-2,-1)="a"*+!U{0}
\POS(2,0)="b"*+!U{+}
\POS(0,2)="c"*+!D{-}
\POS"a"\ar@{--}@[|(3)]"b"
\POS"b"\ar@{--}@[|(3)]"c"
\POS"c"\ar@{-}@[|(3)]"a"
\POS(4.7,-2.5)*{\bullet},*+!R{\vec y}
\end{xy}
& =
\intercol=0.48cm
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-2,-1)="a"*+!U{0}
\POS(0,2)="c"*+!D{-}
\POS(4,-2)="w"*+!L{\vec w}
\POS"c"\ar@{-}@[|(3)]"a"
\POS"a"\ar@{--}@[|(3)]"w"
\POS"c"\ar@{--}@[|(3)]"w"
\POS(4.7,-2.5)*{\bullet},*+!R{\vec y}
\end{xy}
-
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-2,-1)="a"*+!U{0}
\POS(2,0)="b"*+!U{+}
\POS(4,-2)="w"*+!L{\vec w}
\POS"a"\ar@{-}@[|(3)]"b"
\POS"a"\ar@{--}@[|(3)]"w"
\POS"b"\ar@{--}@[|(3)]"w"
\POS(4.7,-2.5)*{\bullet},*+!R{\vec y}
\end{xy}
\end{align*}
\caption{Examples of decompositions in primal space.}
\label{fig:primal:examples}
\end{figure}

To see that there is a $\vec y$ satisfying the above constraints,
we need to show that $R \cap S$ is non-empty, with
$S = \{ \vec y \mid \sp{{\tilde n}_{i_kj_k}}y > 0 \text{ for all $k$}\}$.
It will be easier to show this set is non-empty when the $\vec u_i$ form
an orthogonal basis.  Applying a non-singular linear transformation $T$
does not change the decomposition of $\vec w$ in terms of the $\vec u_i$ (i.e., the
$\alpha_i$ remain unchanged), nor does this change
any of the scalar products in the constraints that define $R \cap S$
(the normals are transformed by $\left(T^{-1}\right)^\T$).
Finding a vector $\vec y \in T(R \cap S)$ ensures that
$T^{-1}(\vec y) \in R \cap S$.
Without loss of generality, we can therefore assume for the purpose of
showing that $R \cap S$ is non-empty that
the $\vec u_i$ indeed form an orthogonal basis.

In the orthogonal basis, we have $\vec b_i^* = \vec u_i$
and the corresponding inward normal $\vec N_i$ is either
$\vec u_i$ or $-\vec u_i$.
Furthermore, each normal of a facet of $S$ of the first type is of the
form $\vec {\tilde n}_{i_kj_k} = a_k \vec u_{i_k} - b_k \vec u_{j_k}$, with 
$a_k, b_k > 0$ and ${i_k} < {j_k}$,
while for the second type each normal is of the form
$\vec {\tilde n}_{i_kj_k} = -a_k \vec u_{i_k} - b_k \vec u_{j_k}$, with 
$a_k, b_k > 0$.
If $\vec {\tilde n}_{i_kj_k} = a_k \vec u_{i_k} - b_k \vec u_{j_k}$
is the normal of a facet of $S$
then either
$(\vec N_{i_k}, \vec N_{j_k}) = (\vec u_{i_k}, \vec u_{j_k})$
or
$(\vec N_{i_k}, \vec N_{j_k}) = (-\vec u_{i_k}, -\vec u_{j_k})$.
Otherwise, the facet would not cut $R$.
Similarly,
if $\vec {\tilde n}_{i_kj_k} = -a_k \vec u_{i_k} - b_k \vec u_{j_k}$
is the normal of a facet of $S$
then either
$(\vec N_{i_k}, \vec N_{j_k}) = (\vec u_{i_k}, -\vec u_{j_k})$
or
$(\vec N_{i_k}, \vec N_{j_k}) = (-\vec u_{i_k}, \vec u_{j_k})$.
Assume now that $R \cap S$ is empty, then there exist
$\lambda_k, \mu_i \ge 0$ not all zero
such that
$\sum_k \lambda_k \vec {\tilde n}_{i_kj_k} + \sum_l \mu_i \vec N_i = \vec 0$.
Assume $\lambda_k > 0$ for some facet of the first type.
If $\vec N_{j_k} = -\vec u_{j_k}$, then $-b_k$ can only be canceled
by another facet $k'$ of the first type with $j_k = i_{k'}$, but then
also $\vec N_{j_{k'}} = -\vec u_{j_{k'}}$.  Since the $j_k$ are strictly
increasing, this sequence has to stop with a strictly positive coefficient
for the largest $\vec u_{j_k}$ in this sequence.
If, on the other hand, $\vec N_{i_k} = \vec u_{i_k}$, then $a_k$ can only
be canceled by the normal of a facet $k'$ of the second kind
with $i_k = j_{k'}$, but then
$\vec N_{i_{k'}} = -\vec u_{i_{k'}}$ and we return to the first case.
Finally, if $\lambda_k > 0$ only for normals of facets of the second type,
then either $\vec N_{i_k} = -\vec u_{i_k}$ or $\vec N_{j_k} = -\vec u_{j_k}$
and so the coefficient of one of these basis vectors will be strictly
negative.
That is, the sum of the normals will never be zero and
the set $R \cap S$ is non-empty.

For each ray $\vec u_j$ of cone $K_i$, i.e., the cone with $\vec u_i$ replaced
by $\vec w$, we now need to determine whether the facet not containing this
ray is closed or not.  We denote the (inward) normal of this cone by
$\vec n_{ij}$.  Note that cone $K_j$ (if it appears in (\ref{eq:triangulations}),
i.e., $\alpha_j \ne 0$) has the same facet opposite $\vec u_i$
and its normal $\vec n_{ji}$ will be equal to either $\vec n_{ij}$ or
$-\vec n_{ij}$, depending on whether we are dealing with an ``external'' facet,
i.e., a facet of $K'$, or an ``internal'' facet.
If, on the other hand, $\alpha_j = 0$, then $\vec n_{ij} = \vec n_{0j}$.
If $\sp{n_{ij}}y > 0$, then the facet is closed.
Otherwise it is open.
It follows that the two (or more) occurrences of external facets are either all open
or all closed, while for internal facets, exactly one is closed.

First consider the facet not containing $\vec u_0 = \vec w$.
If $\alpha_i > 0$, then $\vec u_i$ and $\vec w$ are on the same side of the facet
and so $\vec n_{i0} = \vec n_{0i}$.  Otherwise, $\vec n_{i0} = -\vec n_{i0}$.
Second, if $\alpha_j = 0$, then replacing $\vec u_i$ by $\vec w$ does not
change the affine hull of the facet and so $\vec n_{ij} = \vec n_{0j}$.
Now consider the case that $\alpha_i \alpha_j < 0$, i.e., $\vec u_i$
and $\vec u_j$ are on the same side of the hyperplane through the other rays.
If we project $\vec u_i$, $\vec u_j$ and $\vec w$ onto a plane orthogonal
to the ridge through the other rays, then the possible locations of $\vec w$
with respect to $\vec u_i$ and $\vec u_j$ are shown in Figure~\ref{fig:w:same}.
If both $\vec n_{0i}$ and $\vec n_{0j}$ are closed then $\vec y$ lies in region~1
and therefore $\vec n_{ij}$ (as well as $\vec n_{ji}$) is closed too.
Similarly, if both $\vec n_{0i}$ and $\vec n_{0j}$ are open then so is
$\vec n_{ij}$.  If only one of the facets is closed, then, as explained above,
we choose $\vec n_{ij}$ to be open, i.e., we take $\vec y$ to lie in region~3
or~5.
Figure~\ref{fig:w:opposite} shows the possible configurations
for the case that $\alpha_i \alpha_j > 0$.
If exactly one of $\vec n_{0i}$ and $\vec n_{0j}$ is closed, then
$\vec y$ lies in region~3 or region~5 and therefore $\vec n_{ij}$ is closed iff
$\vec n_{0j}$ is closed.
Otherwise, as explained above, we choose $\vec n_{ij}$ to be closed if $i < j$.

\begin{figure}
\intercol=0.7cm
\begin{minipage}{0cm}
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-4,0)="a"\ar@[|(2)]@{-}(4,0)
\POS?(0)/\intercol/="b"\POS(2,0)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0i}}(0,0.75)}
\POS(2,0)*{\bullet},*+!U{\vec u_j}
\POS(-2,-3)="a"\ar@[|(2)]@{-}(2,3)
\POS?(0)/\intercol/="b"\POS(1,1.5)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0j}}(0,-0.75)}
\POS(1,1.5)*{\bullet},*+!R{\vec u_i}
\POS(-2,3)="a"\ar@{-}(2,-3)
\POS?(0)/\intercol/="b"\POS(1.5,-2.25)
*\xybox{"b"-"a":(0,0)\ar_{\vec n_{ji}}^{\vec n_{ij}}(0,+0.75)}
\POS(1.5,-2.25)*{\bullet},*+!R{\vec u_0 = \vec w}
\POS(0,0)*{\bullet}
\POS(3,1.5)*+[o][F]{\scriptstyle 1}
\POS(0,2.5)*+[o][F]{\scriptstyle 2}
\POS(-3,1.5)*+[o][F]{\scriptstyle 3}
\POS(-3,-1.5)*+[o][F]{\scriptstyle 4}
\POS(0,-3)*+[o][F]{\scriptstyle 5}
\POS(3,-1.5)*+[o][F]{\scriptstyle 6}
\end{xy}
\end{minipage}
\begin{minipage}{0cm}
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-4,0)="a"\ar@[|(2)]@{-}(4,0)
\POS?(0)/\intercol/="b"\POS(2,0)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0i}}(0,0.75)}
\POS(2,0)*{\bullet},*+!U{\vec u_j}
\POS(-2,-3)="a"\ar@[|(2)]@{-}(2,3)
\POS?(0)/\intercol/="b"\POS(1,1.5)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0j}}(0,-0.75)}
\POS(1,1.5)*{\bullet},*+!R{\vec u_i}
\POS(-2,3)="a"\ar@{-}(2,-3)
\POS?(0)/\intercol/="b"\POS(-1.5,2.25)
*\xybox{"b"-"a":(0,0)\ar_{\vec n_{ji}}^{\vec n_{ij}}(0,+0.75)}
\POS(-1.5,2.25)*{\bullet},*+!R{\vec u_0 = \vec w}
\POS(0,0)*{\bullet}
\POS(3,1.5)*+[o][F]{\scriptstyle 1}
\POS(0,2.5)*+[o][F]{\scriptstyle 2}
\POS(-3,1.5)*+[o][F]{\scriptstyle 3}
\POS(-3,-1.5)*+[o][F]{\scriptstyle 4}
\POS(0,-3)*+[o][F]{\scriptstyle 5}
\POS(3,-1.5)*+[o][F]{\scriptstyle 6}
\end{xy}
\end{minipage}
\caption{Possible locations of $\vec w$ with respect to $\vec u_i$ and
$\vec u_j$, projected onto a plane orthogonal to the other rays, when
$\alpha_i \alpha_j < 0$.}
\label{fig:w:same}
\end{figure}

\begin{figure}
\intercol=0.7cm
\begin{minipage}{0cm}
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-4,0)="a"\ar@[|(2)]@{-}(4,0)
\POS?(0)/\intercol/="b"\POS(2,0)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0i}}(0,0.75)}
\POS(2,0)*{\bullet},*+!U{\vec u_j}
\POS(-2,3)="a"\ar@[|(2)]@{-}(2,-3)
\POS?(0)/\intercol/="b"\POS(-1,1.5)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0j}}(0,+0.75)}
\POS(-1,1.5)*{\bullet},*+!R{\vec u_i}
\POS(-2,-3)="a"\ar@{-}(2,3)
\POS?(0)/\intercol/="b"\POS(1.5,2.25)
*\xybox{"b"-"a":(0,0)\ar^{\vec n_{ji}}(0,+0.75)
\POS(0,0)\ar_{\vec n_{ij}}(0,-0.75)}
\POS(1.5,2.25)*{\bullet},*+!L{\vec u_0 = \vec w}
\POS(0,0)*{\bullet}
\POS(3,1.5)*+[o][F]{\scriptstyle 1}
\POS(0,2.5)*+[o][F]{\scriptstyle 2}
\POS(-3,1.5)*+[o][F]{\scriptstyle 3}
\POS(-3,-1.5)*+[o][F]{\scriptstyle 4}
\POS(0,-3)*+[o][F]{\scriptstyle 5}
\POS(3,-1.5)*+[o][F]{\scriptstyle 6}
\end{xy}
\end{minipage}
\begin{minipage}{0cm}
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-4,0)="a"\ar@[|(2)]@{-}(4,0)
\POS?(0)/\intercol/="b"\POS(2,0)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0i}}(0,0.75)}
\POS(2,0)*{\bullet},*+!U{\vec u_j}
\POS(-2,3)="a"\ar@[|(2)]@{-}(2,-3)
\POS?(0)/\intercol/="b"\POS(-1,1.5)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0j}}(0,+0.75)}
\POS(-1,1.5)*{\bullet},*+!R{\vec u_i}
\POS(-2,-3)="a"\ar@{-}(2,3)
\POS?(0)/\intercol/="b"\POS(-1.5,-2.25)
*\xybox{"b"-"a":(0,0)\ar^{\vec n_{ji}}(0,+0.75)
\POS(0,0)\ar_{\vec n_{ij}}(0,-0.75)}
\POS(-1.5,-2.25)*{\bullet},*+!L{\vec u_0 = \vec w}
\POS(0,0)*{\bullet}
\POS(3,1.5)*+[o][F]{\scriptstyle 1}
\POS(0,2.5)*+[o][F]{\scriptstyle 2}
\POS(-3,1.5)*+[o][F]{\scriptstyle 3}
\POS(-3,-1.5)*+[o][F]{\scriptstyle 4}
\POS(0,-3)*+[o][F]{\scriptstyle 5}
\POS(3,-1.5)*+[o][F]{\scriptstyle 6}
\end{xy}
\end{minipage}
\caption{Possible locations of $\vec w$ with respect to $\vec u_i$ and
$\vec u_j$, projected onto a plane orthogonal to the other rays, when
$\alpha_i \alpha_j > 0$.}
\label{fig:w:opposite}
\end{figure}

The algorithm is summarized in Algorithm~\ref{alg:closed}, where
we use the convention that in cone $K_i$, $\vec u_i$ refers to
$\vec u_0 = \vec w$.
Note that we do not need any of the rays or normals in this code.
The only information we need is the closedness of the facets in the
original cone and the signs of the $\alpha_i$.

\begin{algorithm}
\begin{tabbing}
next \= next \= next \= \kill
if $\alpha_j = 0$ \\
\> closed[$K_i$][$\vec u_j$] := closed[$\tilde K$][$\vec u_j$] \\
else if $i = j$ \\
\> if $\alpha_j > 0$ \\
\> \> closed[$K_i$][$\vec u_j$] := closed[$\tilde K$][$\vec u_j$] \\
\> else \\
\> \> closed[$K_i$][$\vec u_j$] := $\lnot$closed[$\tilde K$][$\vec u_j$] \\
else if $\alpha_i \alpha_j > 0$ \\
\> if closed[$\tilde K$][$\vec u_i$] = closed[$\tilde K$][$\vec u_j$] \\
\> \> closed[$K_i$][$\vec u_j$] := $i < j$ \\
\> else \\
\> \> closed[$K_i$][$\vec u_j$] := closed[$\tilde K$][$\vec u_j$] \\
else \\
\> closed[$K_i$][$\vec u_j$] := closed[$\tilde K$][$\vec u_i$] and
closed[$\tilde K$][$\vec u_j$]
\end{tabbing}
\caption{Determine whether the facet opposite $\vec u_j$ is closed in $K_i$.}
\label{alg:closed}
\end{algorithm}

\subsection{Triangulation in primal space}
\label{s:triangulation}

As in the case for Barvinok's decomposition (Section~\ref{s:decomposition}),
we can transform a triangulation of a (closed) cone into closed simple cones
into a triangulation of half-open simple cones that fully partitions the
original cone, i.e., such that the half-open simple cones do not intersect at their
facets.
Again, we apply Proposition~\ref{p:inclusion-exclusion} with $\vec y$
an interior point of the cone (Section~\ref{s:interior}).

\subsection{Multivariate quasi-polynomials as lists of polynomials}

There are many definitions for a (univariate) \ai{quasi-polynomial}.
\shortciteN{Ehrhart1977} uses a definition based on {\em periodic number}s.

\begin{definition}
\label{d:periodic:1}
A rational \defindex{periodic number} $U(p)$
is a function $\ZZ \to \QQ$,
such that there exists a \defindex{period} $q$
such that $U(p) = U(p')$ whenever $p \equiv p' \mod q$.
\end{definition}

\begin{definition}
\label{d:qp:1}
A (univariate)
\defindex{quasi-polynomial}\/ $f$ of degree $d$ is
a function
$$
f(n) = c_d(n) \, n^d + \cdots + c_1(n) \, n + c_0
,
$$
where $c_i(n)$ are rational periodic numbers.
I.e., it is a polynomial expression of degree $d$ 
with rational periodic numbers for coefficients.
The \defindex{period} of a quasi-polynomial is the \ac{lcm}
of the periods of its coefficients.
\end{definition}

Other authors (e.g., \shortciteNP{Stanley1986})
use the following definition of a quasi-polynomial.
\begin{definition}
\label{d:qp:1:list}
A function $f : \ZZ \to \QQ$ is
a (univariate) \defindex{quasi-polynomial} of period $q$ if there
exists a list of $q$ polynomials $g_i \in \QQ[T]$ for $0 \le i < q$ such
that
\[
f (s) = g_i(s) \qquad \hbox{if $s \equiv i \mod {q}$}
.
\]
The functions $g_i$ are called the {\em constituents}.
\index{constituent}
\end{definition}

In our implementation, we use Definition~\ref{d:qp:1},
but whereas
\shortciteN{Ehrhart1977} uses a list of $q$ rational
numbers enclosed in square brackets to represent periodic
numbers, our periodic numbers are polynomial expressions
in fractional parts (Section~\ref{a:data}).
These fractional parts naturally extend to multivariate
quasi-polynomials.
The bracketed (``explicit'') periodic numbers can
be extended to multiple variables by nesting them
(e.g., \shortciteNP{Loechner1999}).

Definition~\ref{d:qp:1:list} could be extended in a similar way
by having a constituent for each residue modulo a vector period $\vec q$.
However, as pointed out by \citeN{Woods2006personal}, this may not result
in the minimum number of constituents.
A vector period can be considered as a lattice with orthogonal generators and
the number of constituents is equal to the index or determinant of that lattice.
By considering more general lattices, we can potentially reduce the number
of constituents.
\begin{definition}
\label{d:qp}
A function $f : \ZZ^n \to \QQ$ is
a (multivariate) \defindex{quasi-polynomial} of period $L$ if there
exists a list of $\det L$ polynomials $g_{\vec i} \in \QQ[T_1,\ldots,T_n]$
for $\vec i$ in the fundamental parallelepiped of $L$ such
that
\[
f (\vec s) = g_{\vec i}(\vec s) \qquad \hbox{if $\vec s \equiv \vec i \mod L$}
.
\]
\end{definition}

To compute the period lattice from a fractional representation, we compute
the appropriate lattice for each fractional part and then take their intersection.
Recall that the argument of each fractional part is an affine expression
in the parameters $(\sp a p + c)/m$,
with $\vec a \in \ZZ^n$ and $c, m \in \ZZ$.
Such a fractional part is translation invariant over
any (integer) value of $\vec p$
such that $\sp a p + m t = 0$, for some $\vec t \in \ZZ$.
Solving this homogeneous equation over the integers (in our implementation,
we use \PolyLib/'s \ai[\tt]{SolveDiophantine}) gives the general solution
$$
\begin{bmatrix}
\vec p \\ t
\end{bmatrix}
=
\begin{bmatrix}
U_1 \\ U_2
\end{bmatrix}
\vec x
\qquad\text{for $\vec x \in \ZZ^n$}
.
$$
The matrix $U_1 \in \ZZ^{n \times n}$ then has the generators of
the required lattice as columns.
The constituents are computed by plugging in each integer point
in the fundamental parallelepiped of the lattice.
These points themselves are computed as explained in Section~\ref{s:fundamental}.
Note that for computing the constituents, it is sufficient to take any
representative of the residue class.  For example, we could take
$\vec w^\T = \vec k^\T W$ in the notations of Lemma~\ref{l:fundamental}.

\begin{example}[\shortciteN{Woods2006personal}]
Consider the parametric polytope
$$
P_{s,t}=\{\, x \mid 0 \le x \le (s+t)/2 \,\}
.
$$
The enumerator of $P_{s,t}$ is
$$
\begin{cases}
\frac s 2 + \frac t 2 + 1 &
\text{if $\begin{bmatrix}s \\ t \end{bmatrix} \in
\begin{bmatrix}
-1 & -2 \\ 1 & 0
\end{bmatrix}
\ZZ^2 +
\begin{bmatrix}
0 \\ 0
\end{bmatrix}
$}
\\
\frac s 2 + \frac t 2 + \frac 1 2 &
\text{if $\begin{bmatrix}s \\ t \end{bmatrix} \in
\begin{bmatrix}
-1 & -2 \\ 1 & 0
\end{bmatrix}
\ZZ^2 +
\begin{bmatrix}
-1 \\ 0
\end{bmatrix}
$}
.
\end{cases}
$$
The corresponding output of \ai[\tt]{barvinok\_enumerate} is
\begin{verbatim}
         s + t  >= 0
          1 >= 0

Lattice:
[[-1 1]
[-2 0]
]
[0 0]
( 1/2 * s + ( 1/2 * t + 1 )
 )
[-1 0]
( 1/2 * s + ( 1/2 * t + 1/2 )
 )
\end{verbatim}
\end{example}

\subsection{Left inverse of an affine embedding}
\label{s:inverse}

We often map a polytope onto a lower dimensional space to remove possible
equalities in the polytope.  These maps are typically represented
by the inverse, mapping the coordinates $\vec x'$ of the lower-dimensional
space to the coordinates $\vec x$ of (an affine subspace of) the original space,
i.e.,
$$
\begin{bmatrix}
\vec x \\ 1
\end{bmatrix}
=
\begin{bmatrix}
T & \vec v \\ \vec 0^\T & 1
\end{bmatrix}
\begin{bmatrix}
\vec x' \\ 1
\end{bmatrix}
,
$$
where, as usual in \PolyLib/, we work with homogeneous coordinates.
To obtain the transformation that maps the coordinates of the original
space to the coordinates of the lower dimensional space,
we need to compute the \ai{left inverse} of the above \ai{affine embedding},
i.e., an $A$, $\vec b$ and $d$ such that
$$
d
\begin{bmatrix}
\vec x' \\ 1
\end{bmatrix}
=
\begin{bmatrix}
A & \vec b \\ \vec 0^\T & d
\end{bmatrix}
\begin{bmatrix}
\vec x \\ 1
\end{bmatrix}
$$

To compute this left inverse, we first compute the
(right) \indac{HNF} of T,
$$
\begin{bmatrix}
U_1 \\ U_2
\end{bmatrix}
T
=
\begin{bmatrix}
H \\ 0
\end{bmatrix}
.
$$
The left inverse is then simply
$$
\begin{bmatrix}
d H^{-1}U_1 & -d H^{-1} \vec v \\ \vec 0^\T & d
\end{bmatrix}
.
$$
We often also want a description of the affine subspace that is the range
of the affine embedding and this is given by
$$
\begin{bmatrix}
U_2 & - U_2 \vec v \\ \vec 0^T & 1
\end{bmatrix}
\begin{bmatrix}
\vec x \\ 1
\end{bmatrix}
=
\vec 0
.
$$
This computation is implemented in \ai[\tt]{left\_inverse}.

\subsection{Integral basis of the orthogonal complement of a linear subspace}
\label{s:completion}

Let $M_1 \in \ZZ^{m \times n}$ be a basis of a linear subspace.
We first extend $M_1$ with zero rows to obtain a square matrix $M'$
and then compute the (left) \indac{HNF} of $M'$,
$$
\begin{bmatrix}
M_1 \\ 0
\end{bmatrix}
=
\begin{bmatrix}
H & 0 \\ 0 & 0
\end{bmatrix}
\begin{bmatrix}
Q_1 \\ Q_2
\end{bmatrix}
.
$$
The rows of $Q_2$ span the orthogonal complement of the given subspace.
Since $Q_2$ can be extended to a unimodular matrix, these rows form
an integral basis.

If the entries on the diagonal of $H$ are all $1$ then $M_1$
can be extended to a unimodular matrix, by concatenating $M_1$ and $Q_2$.
The resulting matrix is unimodular, since
$$
\begin{bmatrix}
M_1 \\ Q_2
\end{bmatrix}
=
\begin{bmatrix}
H & 0 \\ 0 & I_{n-m,n-m}
\end{bmatrix}
\begin{bmatrix}
Q_1 \\ Q_2
\end{bmatrix}
.
$$
This method for extending a matrix of which
only a few lines are known to a \ai{unimodular matrix}
is more general than the method described by \shortciteN{Bik1996PhD},
which only considers extending a matrix given by a single row.

\subsection{Ensuring a polyhedron has only revlex-positive rays}
\label{s:revlexpos}

The \ai[\tt]{barvinok\_series\_with\_options} function and all
further \ai[\tt]{gen\_fun} manipulations assume that the effective
parameter domain has only \ai{revlex-positive} rays.
When used to computer \rgf/s, the \ai[\tt]{barvinok\_enumerate}
application will therefore transform the effective parameter domain
of a problem if it has revlex-negative rays.
It will then not compute the generating function
$$
f(\vec x) = \sum_{\vec p \in \ZZ^m} \#(P_{\vec p} \cap \ZZ^d) \, x^{\vec p}
,
$$
but
$$
g(\vec z) = \sum_{\vec p' \in \ZZ^n}
                \#(P_{T \vec p' + \vec t} \cap \ZZ^d) \, x^{\vec p'}
$$
instead, where $\vec p = T \vec p' + \vec t$,
with $T \in \ZZ^{m \times n}$ and $\vec t \in \ZZ^m$, is an affine transformation
that maps the transformed parameter space back to the original parameter space.

First assume that the parameter domain does not contain any lines and
that there are no equalities in the description of $P_{\vec p}$ that force
the values of $\vec p$ for which $P_{\vec p}$ contains integer points
to lie on a non-standard lattice.
Let the effective parameter domain be given as
$
\{\, \vec p \mid A \vec p + \vec c \ge \vec 0 \,\}
$,
where $A \in \ZZ^{s \times d}$ of row rank $d$;
otherwise the effective parameter domain would contain a line.
Let $H$ be the (left) \indac{HNF} of $A$, i.e.,
$$
A = H Q
,
$$
with $H$ lower-triangular with positive diagonal elements and
$Q$ unimodular.
Let $\tilde Q$ be the matrix obtained from $Q$ by reversing its rows,
and, similarly, $\tilde H$ from $H$ by reversing the columns.
After performing the transformation
$\vec p' = \tilde Q \vec p$, i.e.,
$\vec p = \tilde Q^{-1} \vec p'$, the transformed parameter domain
is given by
$$
\{\, \vec p' \mid A \tilde Q^{-1} \vec p' + \vec c \ge \vec 0 \,\}
$$
or
$$
\{\, \vec p' \mid \tilde H \vec p' + \vec c \ge \vec 0 \,\}
.
$$
The first constraint of this domain is
$h_{11} p'_m + c_1 \ge 0$.  A ray with non-zero final coordinate
therefore has a positive final coordinate.
Similarly, the second constraint is
$h_{22} p'_{m-1} + h_{21} p'_m + c_2 \ge 0$.
A ray with zero $n$th coordinate, but non-zero $n-1$st coordinate,
will therefore have a positive $n-1$st coordinate.
Continuing this reasoning, we see that all rays in the transformed
domain are revlex-positive.

If the parameter domain does contains lines, but is not restricted
to a non-standard lattice, then the number of points in the parametric
polytope is invariant over a translation along the lines.
It is therefore sufficient to compute the number of points in the
orthogonal complement of the linear subspace spanned by the lines.
That is, we apply a prior transformation that maps a reduced parameter
domain to this subspace,
$$
\vec p = L^\perp \vec p' =
\begin{bmatrix}
L & L^\perp
\end{bmatrix}
\begin{bmatrix}
0 \\ I
\end{bmatrix}
\vec p'
,
$$
where $L$ has the lines as columns, and $L^\perp$ an integral basis
for the orthogonal complement (Section~\ref{s:completion}).
Note that the inverse transformation
$$
\vec p' =
L^{-\perp}
\vec p =
\begin{bmatrix}
0 & I
\end{bmatrix}
\begin{bmatrix}
L & L^\perp
\end{bmatrix}^{-1}
\vec p
$$
has integral coefficients since $L^\perp$ can be extended to a unimodular matrix.

If the parameter values $\vec p$ for which $P_{\vec p}$ contains integer points
are restricted to a non-standard lattice, we first replace the parameters
by a different set of parameters that lie on the standard lattice
through ``\ai{parameter compression}''\shortcite{Meister2004PhD},
$$
\vec p = C \vec p'
.
$$
The (left) inverse of $C$ can be computes as explained in
Section~\ref{s:inverse}, giving
$$
\vec p' = C^{-L} \vec p
.
$$
We have to be careful to only apply this transformation when
both the equalities computed in Section~\ref{s:inverse} are satisfied
and some additional divisibility constraints.
In particular if $\vec a^\T/d$ is a row of $C^{-L}$, with $\vec a \in \ZZ^{n'}$
and $d \in \ZZ$, the transformation can only be applied to parameter values
$\vec p$ such that $d$ divides $\sp a p$.

The complete transformation is given by
$$
\vec p = C L^\perp \hat Q^{-1} \vec p'
$$
or
$$
\vec p' = \hat Q L^{-\perp} C^{-L} \vec p
.
$$

\subsection{Parametric Volume Computation}

The \ai{volume} of a (parametric) polytope can serve as an approximation
for the number of integer points in the polytope.
We basically follow the description of~\shortciteN{Rabl2006} here, except that we
focus on volume computation for {\em linearly}
parametrized polytopes, which we exploit to determine the sign
of the determinants we compute, as explained below.

Note first that
the vertices of a linearly parametrized polytope are affine expressions
in the parameters that may be valid only in parts (chambers)
of the parameter domain.
Since the volume computation is based on the (active) vertices, we perform
the computation in each chamber separately.
Also note that since the vertices are affine expressions, it is
easy to check whether they belong to a facet.

The volume of a $d$-simplex, i.e., a $d$-dimensional polytope with
$d+1$ vertices, is relatively easy to compute.
In particular, if $\vec v_i(\vec p)$, for $0 \le i \le d$,
are the (parametric) vertices
of the simplex $P$ then
\begin{equation}
\label{eq:volume}
\vol P = 
\frac 1 {d!}
\left|
\det
\begin{bmatrix}
v_{11}(\vec p) - v_{01}(\vec p) &
v_{12}(\vec p) - v_{02}(\vec p) &
\ldots &
v_{1d}(\vec p) - v_{0d}(\vec p)
\\
v_{21}(\vec p) - v_{01}(\vec p) &
v_{22}(\vec p) - v_{02}(\vec p) &
\ldots &
v_{2d}(\vec p) - v_{0d}(\vec p)
\\
\vdots & \vdots & \ddots & \vdots
\\
v_{d1}(\vec p) - v_{01}(\vec p) &
v_{d2}(\vec p) - v_{02}(\vec p) &
\ldots &
v_{dd}(\vec p) - v_{0d}(\vec p)
\end{bmatrix}
\right|.
\end{equation}

If $P$ is not a simplex, i.e., $N > d+1$, with $N$ the number of
vertices of $P$, then the standard way of computing the volume of $P$
is to first {\em triangulate} $P$, i.e., subdivide $P$ into simplices,
and then to compute and sum the volumes of the resulting simplices.
One way of computing a triangulation is to
compute the \ai{barycenter}
$$
\frac 1 N \sum_i \vec v_i(\vec p)
$$
of $P$
and to perform a subdivision by computing the convex hulls
of the barycenter with each of the facets of $P$.
If a given facet of $P$ is itself a simplex, then this convex hull
is also a simplex.  Otherwise the facet is further subdivided.
This recursive process terminates as every $1$-dimensional polytope
is a simplex.

The triangulation described above is known as
the boundary triangulation~\shortcite{Bueler2000exact} and is used
by \shortciteN{Rabl2006} in his implementation.
The Cohen-Hickey triangulation~\shortcite{Cohen1979volumes,Bueler2000exact}
is a much more efficient variation and uses one of the vertices
instead of the barycenter.  The facets incident on the vertex
do not have to be considered in this case because the resulting subpolytopes
would have zero volume.
Another possibility is to use a
``lifting'' triangulation~\shortcite{Lee1991,DeLoera1995}.
In this triangulation, each vertex is assigned a (random) ``height'' in
an extra dimension.
The projection of the ``lower envelope'' of the resulting polytope onto
the original space results in a subdivision, which is a triangulation
with very high probability.

A complication with the lifting triangulation is that the constraint system
of the lifted polytope will in general not be linearly parameterized,
even if the original polytope is.
It is, however, sufficient to perform the triangulation for a particular
value of the parameters inside the chamber since the parametric polytope
has the same combinatorial structure throughout the chamber.
The triangulation obtained for the instantiated vertices can then
be carried over to the corresponding parametric vertices.
We only need to be careful to select a value for the parameters that
does not lie on any facet of the chambers.  On these chambers, some
of the vertices may coincide.
For linearly parametrized polytopes, it is easy to find a parameter
point in the interior of a chamber, as explained in Section~\ref{s:interior}.
Note that this point need not be integer.

A direct application of the above algorithm, using any of the triangulations,
would yield for each chamber
a volume expressed as the sum of the absolute values of polynomials in
the parameters.  To remove the absolute value, we plug in a particular
value of the parameters (not necessarily integer)
belonging to the given chamber for which we know that the volume is non-zero.
Again, it is sufficient to take any point in the interior of the chamber.
The sign of the resulting value then determines the sign of the whole
polynomial since polynomials are continuous functions and will not change
sign without passing through zero.

\subsection{Maclaurin series division}
\label{s:division}

If $P(t)$ and $Q(t)$ are two Maclaurin series
\begin{align*}
P(t) & =  a_0 + a_1 t + a_2 t^2 + \cdots \\
Q(t) & =  b_0 + b_1 t + b_2 t^2 + \cdots
,
\end{align*}
then, as outlined by \shortciteN[241--247]{Henrici1974},
we can compute the coefficients $c_l$ in
$$
\frac{P(t)}{Q(t)} =: c_0 + c_1 t + c_2 t^2 + \cdots
$$
by applying the recurrence relation
$$
c_l = \frac 1 {b_0} \left( a_l - \sum_{i=1}^l b_i c_{l-i} \right)
.
$$
To avoid dealing with denominators, we can also compute
$d_l = b_0^{l+1} c_l$ instead as
$$
d_l = b_0^l a_l - \sum_{i=1}^l b_0^{i-1} b_i c_{l-i}
.
$$
The coefficients $c_l$ can then be directly read off as
$$
c_l = \frac{d_l}{b_0^{l+1}}
.
$$

\subsection{Specialization through exponential substitution}

This section draws heavily from \shortciteN{Koeppe2006experiments}.

After computing the \rgf/ of a polytope,
\begin{equation}
\label{eq:rgf}
f(\vec x)=
\sum_{i\in I}\alpha_i
    \frac{\sum_{k=1}^{r} \vec x^{\vec w_{ik} }}
         {\prod_{j=1}^{d}\left(1-\vec x^{\vec b_{ij}}\right)}
,
\end{equation}
the number of lattice points in the polytope can be obtained
by evaluating $f(\vec 1)$.  Since $\vec 1$ is a pole of each
term, we need to compute the constant term in the Laurent expansions
of each term in~\eqref{eq:rgf} about $\vec 1$.
Since it is easier to work with univariate series, a substitution is usually
applied, either a \ai{polynomial substitution}
$$
\vec x = (1+t)^{\vec \lambda}
,
$$
as implemented in \LattE/ \shortcite{latte1.1},
or an \ai{exponential substitution} (see, e.g., \shortciteNP{Barvinok1999}),
$$
\vec x = e^{t \vec \lambda}
,
$$
as implemented in \LattEmk/ \shortcite{latte-macchiato}.
In each case, $\vec \lambda \in \ZZ^d$ is a vector that is not orthogonal
to any of the $\vec b_{ij}$.
Both substitutions also transform the problem of computing the
constant term in the Laurent expansions about $\vec x = \vec 1$
to that of computing the constant term in the
Laurent expansions about $t = 0$.

Consider now one of the terms in~\eqref{eq:rgf},
$$
g(t) =
    \frac{\sum_{k=1}^{r} e^{a_k t}}
         {\prod_{j=1}^{d}\left(1-e^{c_j t}\right)}
,
$$
with $a_k = \sp{w_{ik}}{\lambda}$ and $c_j = \sp{b_{ij}}{\lambda}$.
We rewrite this equation as
$$
g(t) =
(-1)^d
    \frac{\sum_{k=1}^{r} e^{a_k t}}
         {t^d \prod_{j=1}^d c_j}
    \prod_{j=1}^d \frac{-c_j t}
         {1-e^{c_j t}}
.
$$
The second factor is analytic in a neighborhood of the origin
$x = c_1 = \cdots = c_d = 0$ and therefore has a Taylor series expansion
\begin{equation}
\label{eq:todd}
    \prod_{j=1}^d \frac{-c_j t}
         {1-e^{c_j t}}
=
    \sum_{m=0}^{\infty} \todd_m(-c_1, \ldots, -c_d) t^m
,
\end{equation}
where $\todd_m$ is a homogeneous polynomial of degree $m$ called
the $m$-th \ai{Todd polynomial}~\cite{Barvinok1999}.
Also expanding the numerator in the first factor, we find
$$
g(t)  = \frac{(-1)^d}{t^d \prod_{j=1}^d c_j}
            \left(
                \sum_{n=0}^{\infty}\frac{\sum_{k=1}^{r} a_k^n}{n!}
            \right)
            \left(
                \sum_{m=0}^{\infty} \todd_m(-c_1, \ldots, -c_d) t^m
            \right)
,
$$
with constant term
\begin{multline}
\label{eq:todd:constant}
\frac{(-1)^d}{t^d \prod_{j=1}^d c_j}
    \left(\sum_{i=0}^d \frac{\sum_{k=1}^{r} a_k^i}{i!}
                    \todd_{d-i}(-c_1, \ldots, -c_d)\right)t^d
= \\
\frac{(-1)^d}{\prod_{j=1}^d c_j}
    \sum_{i=0}^d \frac{\sum_{k=1}^{r} a_k^i}{i!} \todd_{d-i}(-c_1, \ldots, -c_d)
.
\end{multline}
To compute the first $d+1$ terms in the Taylor series~\eqref{eq:todd},
we write down the truncated Taylor series
$$
\frac{e^t -1}t \equiv
        \sum_{i=0}^d \frac 1{(i+1)!} t^i \equiv
        \frac 1 {(d+1)!} \sum_{i=0}^d \frac{(d+1)!}{(i+1)!} t^i
\mod t^{d+1}
,
$$
where we have
$$
        \frac 1 {(d+1)!} \sum_{i=0}^d \frac{(d+1)!}{(i+1)!} t^i
                \in \frac 1{(d+1)!} \ZZ[t]
.
$$
Computing the reciprocal as explained in Section~\ref{s:division},
we find
\begin{equation}
\label{eq:t-exp-1}
\frac{t}{e^t-1} = \frac 1{\frac{e^t -1}t}
\equiv (d+1)! \frac 1{\sum_{i=0}^d \frac{(d+1)!}{(i+1)!} t^i}
=: \sum_{i=0}^d b_i t^i
.
\end{equation}
Note that the constant term of the denominator is $1/(d+1)!$.
The denominators of the quotient are therefore $((d+1)!)^{i+1}/(d+1)!$.
Also note that the $b_i$ are independent of the generating function
and can be computed in advance.
An alternative way of computing the $b_i$ is to note that
$$
\frac{t}{e^t-1} = \sum_{i=0}^\infty B_i \frac{t^i}{i!}
,
$$
with $B_i = i! \, b_i$ the \ai{Bernoulli number}s, which can be computed
using the recurrence~\eqref{eq:Bernoulli} (see Section~\ref{s:nested}).

Substituting $t$ by $c_j t$ in~\eqref{eq:t-exp-1}, we have
$$
\frac{-c_j t}{1-e^{c_j t}} = \sum_{i=0}^d b_i c_j^i t^i
.
$$
Multiplication of these truncated Taylor series for each $c_j$
results in the first $d+1$ terms of~\eqref{eq:todd},
$$
    \sum_{m=0}^{d} \todd_m(-c_1, \ldots, -c_d) t^m
=:
    \sum_{m=0}^{d} \frac{\beta_m}{((d+1)!)^m} t^m
,
$$
from which
it is easy to compute the constant term~\eqref{eq:todd:constant}.
Note that this convolution can also be computed without the use
of rational coefficients,
$$
\frac{(-1)^d}{\prod_{j=1}^d c_j}
    \sum_{i=0}^d \frac{\alpha_i}{i!} \frac{\beta_{d-i}}{((d+1)!)^{d-i}}
=
\frac{(-1)^d}{((d+1)!)^d\prod_{j=1}^d c_j}
    \sum_{i=0}^d \left(\frac{((d+1)!)^i}{i!}\alpha_i\right) \beta_{d-i}
,
$$
with $\alpha_i = \sum_{k=1}^{r} a_k^i$.

\begin{example}
Consider the \rgf/
\begin{multline*}
\f T x =
\frac{x_1^2}{(1-x_1^{-1})(1-x_1^{-1}x_2)}
+
\frac{x_2^2}{(1-x_2^{-1})(1-x_1 x_2^{-1})}
+ {} \\
\frac1{(1-x_1)(1-x_2)}
\end{multline*}
from \shortciteN[Example~39]{Verdoolaege2005PhD}.
Since this is a 2-dimensional problem, we first compute the first
3 Todd polynomials (evaluated at $-1$),
$$
\frac{e^t -1}t \equiv
1 + \frac 1 2 t + \frac 1 7 t^2 =
\frac 1 6
\begin{bmatrix}
6 & 3 & 1
\end{bmatrix}
$$
and
$$
\frac t{e^t -1} \equiv
\begin{bmatrix}
\displaystyle\frac{-1}{1} & \displaystyle\frac{3}{6} & \displaystyle\frac{-3}{36}
\end{bmatrix}
,
$$
where we represent each truncated power series by a vector of its
coefficients.
The vector $\vec\lambda = (1, -1)$ is not
orthogonal to any of the rays, so we can use the substitution
$\vec x = e^{(1, -1)t}$
and obtain
$$
\frac{e^{2t}}{(1-e^{-t})(1-e^{-2t})}
+
\frac{e^{-2t}}{(1-e^{t})(1-e^{2t})}
+
\frac1{(1-e^{t})(1-e^{-t})}
.
$$
We have
\begin{align*}
\frac{t}{1-e^{- t}} & = 
\begin{bmatrix}
\displaystyle\frac{-1}{1} & \displaystyle\frac{-3}{6} & \displaystyle\frac{-3}{36}
\end{bmatrix}
\\
\frac{2t}{1-e^{-2 t}} & = 
\begin{bmatrix}
\displaystyle\frac{-1}{1} & \displaystyle\frac{-6}{6} & \displaystyle\frac{-12}{36}
\end{bmatrix}
\\
\frac{-t}{1-e^{t}} & = 
\begin{bmatrix}
\displaystyle\frac{-1}{1} & \displaystyle\frac{3}{6} & \displaystyle\frac{-3}{36}
\end{bmatrix}
\\
\frac{-2t}{1-e^{2t}} & = 
\begin{bmatrix}
\displaystyle\frac{-1}{1} & \displaystyle\frac{6}{6} & \displaystyle\frac{-12}{36}
\end{bmatrix}
.
\end{align*}
The first term in the \rgf/ evaluates to
\begin{align*}
&
\frac 1{-1 \cdot -2}
\begin{bmatrix}
\displaystyle\frac{1}{1} & \displaystyle\frac{2}{1} & \displaystyle\frac{4}{2}
\end{bmatrix}
*
\left(
\begin{bmatrix}
\displaystyle\frac{-1}{1} & \displaystyle\frac{-3}{6} & \displaystyle\frac{-3}{36}
\end{bmatrix}
\begin{bmatrix}
\displaystyle\frac{-1}{1} & \displaystyle\frac{-6}{6} & \displaystyle\frac{-12}{36}
\end{bmatrix}
\right)
\\
= {} &
\frac 1{2}
\begin{bmatrix}
\displaystyle\frac{1}{1} & \displaystyle\frac{2}{1} & \displaystyle\frac{4}{2}
\end{bmatrix}
*
\begin{bmatrix}
\displaystyle\frac{1}{1} & \displaystyle\frac{9}{6} & \displaystyle\frac{33}{36}
\end{bmatrix}
\\
= {} &
\frac 1{72}
\begin{bmatrix}
1 & 2 \cdot 6 & 4 \cdot 18
\end{bmatrix}
*
\begin{bmatrix}
1 & 9 & 33
\end{bmatrix}
= \frac {213}{72} = \frac{71}{24}
.
\end{align*}
Due to symmetry, the second term evaluates to the same value,
while for the third term we find
$$
\frac{1}{-1\cdot 1 \cdot 36}
\begin{bmatrix}
1 & 0 \cdot 6 & 0 \cdot 18
\end{bmatrix}
*
\begin{bmatrix}
1 & 0 & -3
\end{bmatrix}
=
\frac{-3}{-36} = \frac 1{12}
.
$$
The sum is
$$
\frac{71}{24} + \frac{71}{24} + \frac 1{12} = 6
.
$$
\end{example}

Note that the run-time complexities of polynomial and exponential
substitution are basically the same.  The experiments of
\citeN{Koeppe2006primal} are somewhat misleading in this respect
since the polynomial substitution (unlike the exponential
substitution) had not been optimized to take full
advantage of the stopped Barvinok decomposition.
For comparison, \autoref{t:hickerson} shows running times
for the same experiments of that paper, but using
barvinok version \verb+barvinok-0.23-47-gaa9024e+
on an Athlon MP 1500+ with 512MiB internal memory.
This machine appears to be slightly slower than the
machine used in the experiments of \citeN{Koeppe2006primal}
as computing {\tt hickerson-14} using the dual decomposition
with polynomial substitution an maximal index 1
took 2768 seconds on this machine using \LattEmk/.
At this stage, it is not clear yet why the number of
cones in the dual decomposition of {\tt hickerson-13}
differs from that of \LattE/~\shortcite{latte1.1} and
\LattEmk/~\cite{latte-macchiato}.
We conclude from \autoref{t:hickerson} that (our implementation of)
the exponential substitution is always slightly faster than
(our implementation of) the polynomial substitution.
The optimal maximal index for these examples is about 500,
which agrees with the experiments of \citeN{Koeppe2006primal}.

\begin{table}
\begin{center}
\begin{tabular}{rrrrrrr}
\hline
&
\multicolumn{3}{c}{Dual decomposition} &
\multicolumn{3}{c}{Primal decomposition}
\\
&
& \multicolumn{2}{c}{Time (s)} &
& \multicolumn{2}{c}{Time (s)}
\\
\cline{3-4}
\cline{6-7}
Max.\ index & Cones & Poly & Exp & Cones & Poly & Exp \\
\hline
\multicolumn{7}{c}{{\tt hickerson-12}}
\\
\hline
1 & 11625 & 9.24 & 8.90 & 7929 & 4.80 & 4.55
\\
10 & 4251 & 4.32 & 4.19 & 803 & 0.66 & 0.62
\\
100 & 980 & 1.42 & 1.35 & 84 & 0.13 & 0.12
\\
200 & 550 & 1.00 & 0.92 & 76 & 0.12 & 0.12
\\
300 & 474 & 0.93 & 0.86 & 58 & 0.12 & 0.10
\\
500 & 410 & 0.90 & 0.83 & 42 & 0.10 & 0.10
\\
1000 & 130 & 0.42 & 0.38 & 22 & {\bf 0.10} & {\bf 0.07}
\\
2000 & 10 & {\bf 0.10} & {\bf 0.10} & 22 & 0.10 & 0.09
\\
5000 & 7 & 0.12 & 0.11 & 7 & 0.12 & 0.10
\\
\hline
\multicolumn{7}{c}{{\tt hickerson-13}}
\\
\hline
1 & 494836 & 489 & 463 & 483507 & 339 & 315
\\
10 & 296151 & 325 & 309 & 55643 & 51 & 48
\\
100 & 158929 & 203 & 192 & 9158 & 11 & 10
\\
200 & 138296 & 184 & 173 & 6150 & 9 & 8
\\
300 & 110438 & 168 & 157 & 4674 & 8 & 7
\\
500 & 102403 & 163 & 151 & 3381  & {\bf 8} & {\bf 7}
\\
1000 & 83421 & {\bf 163} & {\bf 149} & 2490 & 8 & 7
\\
2000 & 77055 & 170 & 153 & 1857 & 10 & 8
\\
5000 & 57265 & 246 & 211 & 1488 & 13 & 11
\\
10000 & 50963 & 319 & 269 & 1011 & 26 & 21
\\
\hline
\multicolumn{7}{c}{{\tt hickerson-14}}
\\
\hline
1 & 1682743 & 2171 & 2064 & 552065 & 508 & 475
\\
10 & 1027619 & 1453 & 1385 & 49632 & 62 & 59
\\
100 & 455474 & 768 & 730 & 8470 & 14 & 13
\\
200 & 406491 & 699 & 661 & 5554 & 11 & 10
\\
300 & 328340 & 627 & 590 & 4332 & 11 & 9
\\
500 & 303566 & 605 & 565 & 3464 & {\bf 11} & {\bf 9}
\\
1000 & 232626 & {\bf 581} & {\bf 532} & 2384 & 12 & 10
\\
2000 & 195368 & 607 & 545 & 1792 & 14 & 12
\\
5000 & 147496 & 785 & 682 & 1276 & 19 & 16
\\
10000 & 128372 & 966 & 824 & 956 & 29 & 23
\\
\hline
\end{tabular}
\caption{Timing results of dual and primal decomposition with
polynomial or exponential substitution on the Hickerson examples}
\label{t:hickerson}
\end{center}
\end{table}

\subsection{Approximate Enumeration using Nested Sums}
\label{s:nested}

If $P \in \QQ^d$ is a polyhedron and $p(\vec x) \in \QQ[\vec x]$ is a
polynomial and we want to sum $p(\vec x)$ over all integer values
of (a subset of) the variables $\vec x$, then we can do this incrementally
by taking a variable $x_1$ with lower bound $L(\vec{\hat x})$
and upper bound $U(\vec{\hat x})$, with $\vec{\hat x} = (x_2, \ldots, x_d)$,
and computing
\begin{equation}
\label{eq:nested:sum}
Q(\vec{\hat x}) = \sum_{x_1 = L(\vec{\hat x})}^{U(\vec{\hat x})} p(\vec x)
.
\end{equation}
Since $P$ is a polytope, the lower bound is a maximum of affine expressions
in the remaining variables, while the upper bound is a minimum of such expressions.
If the coefficients in these expressions are all integer, then we can
compute $Q(\vec{\hat x})$ exactly as a piecewise polynomial using formulas
for sums of powers, as proposed by, e.g.,
\shortciteN{Tawbi1994,Sakellariou1997sums,VanEngelen2004}.
If some of the coefficients are not integer, we can apply the same formulas
to obtain an approximation, which can is some cases be shown
to be an overapproximation~\shortcite{VanEngelen2004}.
Note that if we take the initial polynomial to be the constant $1$, then
this gives us a method for computing an approximation of the number
of integer points in a (parametric) polytope.

The first step is to compute the chamber decomposition of $P$ when viewed
as a 1-dimensional parametric polytope.  That is, we need to partition
the projection of $P$ onto the remaining variables into polyhedral cells
such that in each cell, both the upper and the lower bound are described
by a single affine expression.  Basically, for each pair of lower and upper
bound, we compute the cell where the chosen lower bound is (strictly)
smaller than all other lower bounds and similarly for the upper bound.

For any given pair of lower and upper bound $(l(\vec {\hat x}), u(\vec{\hat x}))$,
the formula~\eqref{eq:nested:sum} is computed for each monomial of $p(\vec x)$
separately.  For the constant term $\alpha_0$, we have
$$
\sum_{x_1 = l(\vec {\hat x})}^{u(\vec{\hat x})} \alpha_0(\vec{\hat x})
    = \alpha_0(\vec{\hat x}) \left(u(\vec{\hat x}) - l(\vec {\hat x}) + 1\right)
.
$$
For the higher degree monomials, we use the formula
\begin{equation}
\label{eq:summation}
\sum_{k=0}^{m-1} k^n = {1\over{n+1}}\sum_{k=0}^n{n+1\choose{k}} B_k m^{n+1-k}
=: S_n(m)
,
\end{equation}
with $B_i$ the \ai{Bernoulli number}s, which can be computed
using the recurrence
\begin{equation}
\label{eq:Bernoulli}
    \sum_{j=0}^m{m+1\choose{j}}B_j = 0
\qquad B_0 = 1
.
\end{equation}
Note that \eqref{eq:summation} is also valid if $m = 0$,
i.e., $S_n(0) = 0$, a fact
that can be easily shown using Newton series~\shortcite{VanEngelen2004}.

\newcounter{saveenumi}

Since we can only directly apply the summation formula when
the lower bound is zero (or one), we need to consider several
cases.
\begin{enumerate}
\item $l(\vec {\hat x}) \ge 1$

\begin{align*}
\sum_{x_1 = l(\vec {\hat x})}^{u(\vec{\hat x})} \alpha_n(\vec{\hat x}) \, x_1^n
& =
\alpha_n(\vec{\hat x})
\left(
\sum_{x_1 = 1}^{u(\vec{\hat x})} x_1^n
-
\sum_{x_1 = 1}^{l(\vec {\hat x})-1} x_1^n
\right)
\\
& =
\alpha_n(\vec{\hat x})
\left( S_n(u(\vec{\hat x})+1) - S_n(l(\vec {\hat x})) \right)
\end{align*}

\item $u(\vec{\hat x}) \le -1$

\begin{align*}
\sum_{x_1 = l(\vec {\hat x})}^{u(\vec{\hat x})} \alpha_n(\vec{\hat x}) \, x_1^n
& =
\alpha_n(\vec{\hat x}) (-1)^n
\sum_{x_1 = -u(\vec {\hat x})}^{-l(\vec{\hat x})} \alpha_n(\vec{\hat x}) \, x_1^n
\\
& =
\alpha_n(\vec{\hat x}) (-1)^n
\left( S_n(-l(\vec{\hat x})+1) - S_n(-u(\vec {\hat x})) \right)
\end{align*}

\item $l(\vec {\hat x}) \le 0$ and $u(\vec{\hat x}) \ge 0$

\begin{align*}
\sum_{x_1 = l(\vec {\hat x})}^{u(\vec{\hat x})} \alpha_n(\vec{\hat x}) \, x_1^n
& =
\alpha_n(\vec{\hat x})
\left(
\sum_{x_1 = 0}^{u(\vec{\hat x})} x_1^n
+
(-1)^n
\sum_{x_1 = 1}^{-l(\vec {\hat x})} x_1^n
\right)
\\
& =
\alpha_n(\vec{\hat x})
\left(
S_n(u(\vec{\hat x})+1)
+
(-1)^n
S_n(-l(\vec{\hat x})+1)
\right)
\end{align*}

\setcounter{saveenumi}{\value{enumi}}
\end{enumerate}

If the coefficients in the lower and upper bound are all
integer, then the above 3 cases partition (the integer points in)
the projection of $P$ onto the remaining variables.
However, if some of the coefficients are rational, then the lower
and upper bound can lie in the open interval $(0,1)$ for some
values of $\vec{\hat x}$.  We therefore also need to consider
the following two cases.
Note that the results will not be correct in these cases, but
not taking them into account would lead to a greater loss in accuracy.

\begin{enumerate}
\setcounter{enumi}{\value{saveenumi}}
\item $0 < l(\vec {\hat x}) < 1$

$$
\sum_{x_1 = l(\vec {\hat x})}^{u(\vec{\hat x})} \alpha_n(\vec{\hat x}) \, x_1^n
=
\alpha_n(\vec{\hat x})
S_n(u(\vec{\hat x})+1)
$$

\item $0 < -u(\vec {\hat x}) < 1$ and $l(\vec {\hat x}) \le 0$

$$
\sum_{x_1 = l(\vec {\hat x})}^{u(\vec{\hat x})} \alpha_n(\vec{\hat x}) \, x_1^n
=
\alpha_n(\vec{\hat x})
(-1)^n
S_n(-l(\vec{\hat x})+1)
$$

\end{enumerate}

\subsection{Summation using local Euler-Maclaurin formula}
\label{s:euler}

\sindex{local}{Euler-Maclaurin formula}
In this section we provide some implementation details
on using \ai{local Euler-Maclaurin formula} to compute
the sum of a piecewise polynomial evaluated in all integer
points of a two-dimensional parametric polytope.
For the theory behind these formula and a discussion
of the original implementation (for non-parametric simplices),
we refer to \shortciteN{Berline2006local}.

In particular, consider a parametric piecewise polynomial
in $n$ parameters and $m$ variables
$c : \ZZ^n \to \ZZ^m \to \QQ : \vec p \mapsto c(\vec p)$,
with $c(\vec p) : \ZZ^m \to \QQ : \vec x \mapsto c(\vec p)(\vec x)$
and
$$
c_{\vec p}(\vec x) =
\begin{cases}
c_1(\vec p)(\vec x) & \text{if $\vec x \in D_1(\vec p)$}
\\
\vdots
\\
c_r(\vec p)(\vec x) & \text{if $\vec x \in D_r(\vec p)$}
,
\end{cases}
$$
with the $c_i$ polynomials, $c_i \in (\QQ[\vec p])[\vec x]$, and 
the $D_i$ disjoint linearly parametric polytopes.
We want to compute
$$
g(\vec p) = \sum_{\vec x \in \ZZ^m} c(\vec p)(\vec x)
.
$$

\subsubsection{Reduction to the summation of a parametric polynomial
over a parametric polytope with a fixed combinatorial structure}

Since the $D_i$ are disjoint, we can consider each
$(c_i, D_i)$-pair individually and compute
$$
g(\vec p) = \sum_{i=1}^r g_i(\vec p) =
\sum_{i=1}^r \sum_{\vec x \in D_r(\vec p) \cap \ZZ^m} c_r(\vec p)(\vec x)
.
$$
The second step is to compute the \ai{chamber decomposition}
~\shortcite[Section 4.2.3]{Verdoolaege2005PhD} of each parametric
polytope $D_i$.
The result is a subdivision of the parameter space into chambers
$C_{ij}$ such that $D_i$ has a fixed combinatorial structure,
in particular a fixed set of parametric vertices,
on (the interior of) each $C_{ij}$.  Applying \autoref{p:inclusion-exclusion},
this subdivision can be transformed into a partition
$\{\, \tilde C_{ij} \,\}$ by
making some of the facets of the chambers open%
~\shortcite[Section~3.2]{Koeppe2007parametric}.
Since we are only interested in integer parameter values,
any of the resulting open facets $\sp a p + c > 0$,
with $\vec a \in \ZZ^n$ and $c \in \ZZ$,
can then be replaced by  $\sp a p + c-1 \ge 0$.
Again, we have
$$
g_i(\vec p) = \sum_j g_{ij}(\vec p) =
\sum_j \sum_{\vec x \in C_{ij}(\vec p) \cap \ZZ^m} c_r(\vec p)(\vec x)
.
$$

After this reduction, the technique of
\shortciteN{Berline2006local} can be applied practically verbatim
to the parametric polytope with a fixed combinatorial structure.
In principle, we could also handle piecewise quasi-polynomials
using the technique of \shortciteN[Section~4.5.4]{Verdoolaege2005PhD},
except that we only need to create an extra variable for each
distinct floor expression in a monomial, rather than for each
occurrence of a floor expression in a monomial.
However, since we currently only support two-dimensional polytopes,
this reduction has not been implemented yet.

\subsubsection{Summation over a one-dimensional parametric polytope}

The basis for the summation technique is the local
Euler-Maclaurin formula~\cite[Theorem~26]{Berline2006local}
\begin{equation}
\label{eq:EML}
\sum_{\vec x \in P(\vec p) \cap \Lambda} h(\vec p)(\vec x)
= \sum_{F(\vec p) \in {\mathcal F}(P(\vec p))}
            \int_{F(\vec p)} D_{P(\vec p),F(\vec p)} \cdot h(\vec p)
,
\end{equation}
where $P(\vec p)$ is a parametric polytope,
$\Lambda$ is a lattice, ${\mathcal F}(P(\vec p))$
are the faces of $P(\vec p)$, $D_{P(\vec p),F(\vec p)}$ is a
specific differential operator associated to the face of a polytope.
The \ai{Lebesgue measure} used in the integral is such that the
integral of the indicator function of a lattice element of
the lattice $\Lambda \cap (\affhull(F(\vec p)) - F(\vec p))$ is 1,
i.e., the intersection of $\Lambda$ with the linear subspace
parallel to the affine hull of the face $F(\vec p)$.
Note that the original theorem is formulated for a non-parametric
polytope and a non-parametric polynomial.  However, as we will see,
in each of the steps in the computation, the parameters can be
treated as symbolic constants without affecting the validity of the formula,
see also~\shortciteN[Section 6]{Berline2006local}.

The differential operator $D_{P(\vec p),F(\vec p)}$ is obtained
by plugging in the vector $\vec D=(D_1,\ldots,D_m)$ of first
order differential operators, i.e., $D_k$ is the first order
differential operator in the $k$th variable,
in the function $\mu_{P(\vec p),F(\vec p)}$.
This function is determined by the \defindex{transverse cone}
of the polyhedron $P(\vec p)$ along its face $F(\vec p)$,
which is the \ai{supporting cone} of $P(\vec p)$ along $F(\vec p)$
projected into the linear subspace orthogonal to $F(\vec p)$.
The lattice associated to this space is the projection of
$\Lambda$ into this space.

In particular, for a zero-dimensional affine cone in the zero-dimensional
space, we have $\mu = 1$~\cite[Proposition 12]{Berline2006local},
while for a one-dimensional affine
cone $K = (-t + \RR_+) r$ in the one-dimensional space, where
$r$ is a primitive integer vector and $t \in [0,1)$,
we have~\cite[(13)]{Berline2006local}
\begin{equation}
\label{eq:mu:1}
\mu(K)(\xi) = \frac{e^{t y}}{1-e^y} + \frac 1{y}
= -\sum_{n=0}^\infty \frac{b(n+1, t)}{(n+1)!} y^n
,
\end{equation}
with $y = \sps \xi r$ and $b(n,t)$ the \ai{Bernoulli polynomial}s
defined by the generating series
$$
\frac{e^{ty} y}{e^y - 1} = \sum_{n=0}^\infty \frac{b(n,t)}{n!} y^n
.
$$
The constant terms of these Bernoulli polynomials
are the \ai{Bernoulli number}s.

Applying \eqref{eq:EML} to a one-dimensional parametric polytope
$P(\vec p) = [v_1(\vec p), v_2(\vec p)]$, we find
$$
\begin{aligned}
\sum_{x \in P(\vec p) \cap \ZZ} h(\vec p)(x)
= & \int_{P(\vec p)} D_{P(\vec p), P(\vec p)} \cdot h(\vec p)
\\
& + \int_{v_1(\vec p)} D_{P(\vec p), v_1(\vec p)} \cdot h(\vec p)
\\
& + \int_{v_2(\vec p)} D_{P(\vec p), v_2(\vec p)} \cdot h(\vec p)
.
\end{aligned}
$$
The transverse cone of a polytope along the whole polytope is
a zero-dimensional cone in a zero-dimensional space and so
$D_{P(\vec p), P(\vec p)} = \mu_{P(\vec p), P(\vec p)}(D) = 1$.
The transverse cone along $v_1(\vec p)$ is $v_1(\vec p) + \RR_+$
and so $D_{P(\vec p), v_1(\vec p)} = \mu(v_1(\vec p) + \RR_+)(D)$
as in \eqref{eq:mu:1}, with $y = \sps D 1 = D$ and
$t = \ceil{v_1(\vec p)} - v_1(\vec p) =
\fractional{-v_1(\vec p)}$.
Similarly we find
$D_{P(\vec p), v_2(\vec p)} = \mu(v_2(\vec p) - \RR_+)(D)$
as in \eqref{eq:mu:1}, with $y = \sps D {-1} = -D$ and
$t = v_2(\vec p) - \floor{v_2(\vec p)} =
\fractional{v_2(\vec p)}$.
Summarizing, we find
$$
\begin{aligned}
\sum_{x \in P(\vec p) \cap \ZZ} h(\vec p)(x)
= & \int_{v_1(\vec p)}^{v_2(\vec p)} h(\vec p)(t) \, dt
\\
& -\sum_{n=0}^\infty \frac{b(n+1, \fractional{-v_1(\vec p)})}{(n+1)!}
        (D^n h(\vec p))(v_1(\vec p))
\\
& -\sum_{n=0}^\infty (-1)^n \frac{b(n+1, \fractional{v_2(\vec p)})}{(n+1)!}
        (D^n h(\vec p))(v_2(\vec p))
.
\end{aligned}
$$

Note that in order to apply this formula, we need to verify
first that $v_1(\vec p)$ is indeed smaller than (or equal to)
$v_2(\vec p)$.  Since the combinatorial structure of $P(\vec p)$
does not change throughout the interior of the chamber, we only
need to check the order of the two vertices for one value
of the parameters from the interior of the chamber, a point
which we may compute as in \autoref{s:interior}.

\subsubsection{Summation over a two-dimensional parametric polytope}

For two-dimensional polytope, formula~\eqref{eq:EML} has three kinds
of contributions: the integral of the polynomial over the polytope,
contributions along edges and contributions along vertices.
As suggested by~\citeN{Berline2007personal}, the integral can be computed
by applying the Green-Stokes theorem:
$$
    \iint_{P(\vec p)}
    \left(\frac{\partial M}{\partial x} - \frac{\partial L}{\partial y}\right) =
    \int_{\partial P(\vec p)} (L\, dx + M\, dy)
.
$$
In particular, if $M(\vec p)(x,y)$ is such that
$\frac{\partial M}{\partial x}(\vec p)(x,y) = h(\vec p)(x,y)$
then
$$
    \iint_{P(\vec p)} h(\vec p)(x,y) =
    \int_{\partial P(\vec p)} M(\vec p)(x,y) \, dy
.
$$
Care must be taken to integrate over the boundary in the positive
direction.  Assuming the vertices of the polygon are not given
in a predetermined order, we can check the correct orientation
of the vertices of each edge individually.  Let $\vec n = (n_1, n_2)$
be the inner normal of a facet and let $\vec v_1(\vec p)$
and $\vec v_2(\vec p)$ be the two vertices of the facet, then
the vertices are in the correct order if
$$
\begin{vmatrix}
v_{2,1}(\vec p)-v_{1,1}(\vec p) & n_1
\\
v_{2,2}(\vec p)-v_{1,2}(\vec p) & n_2
\end{vmatrix}
\ge 0
.
$$
Since these two vertices belong to the same edge, their order
will not change within a chamber and so we can again perform
this check for a single value of the parameters.
To integrate $M$ over an edge $F$, let $\vec f$ be a primitive
integer vector in the direction of the edge.
Then $\vec v_2(\vec p) = \vec v_1(\vec p) + k(\vec p) \, \vec f$
and any point on the edge can be written as
$\vec v_1(\vec p) + \lambda \vec f$ with
$0 \le \lambda \le k(\vec p)$.
That is,
\begin{equation}
\label{eq:EML:int}
\int_F M(\vec p)(x,y) \, dy
=
\int_0^{k(\vec p)}
    M(\vec p)(v_{1,1}(\vec p) + \lambda f_1,
              v_{1,2}(\vec p) + \lambda f_2)
    f_2 \, d\lambda
.
\end{equation}

For the edges, we can again apply \eqref{eq:mu:1}, but we
must first project the supporting cone at the edge into
the linear subspace orthogonal to the edge.
Let $\vec n = (n_1, n_2)$ be the (primitive integer) inner normal
of this facet $F(\vec p)$, then $\vec f = (-n_2, n_1)$ is parallel
to the facet and we can write one of the vertices $\vec v(\vec p)$
as a linear combination of these two vectors:
\begin{equation}
\label{eq:EML:facet:coordinates}
\vec v(\vec p)
=
\begin{bmatrix}
\vec f & \vec n
\end{bmatrix}
\vec a(\vec p)
=
\begin{bmatrix}
-n_2 & n_1 \\
n_1 & n_2
\end{bmatrix}
\vec a(\vec p)
\end{equation}
or
\begin{equation}
\label{eq:EML:facet:coordinates:2}
\vec a(\vec p)
=
\begin{bmatrix}
-n_2 & n_1 \\
n_1 & n_2
\end{bmatrix}^{-1}
\vec v(\vec p)
=
\begin{bmatrix}
-n_2/d & n_1/d \\
n_1/d & n_2/d
\end{bmatrix}
\vec v(\vec p),
\end{equation}
with $d = n_1^2+n_2^2$.
The lattice associated to the linear subspace orthogonal
to the facet is the projection of $\Lambda$ into this space.
Since $\vec n$ is primitive, a basis for this lattice can be
identified with $\vec n/d$.
The coordinate of the whole facet in this space is therefore
$
d a_2(\vec p) =
\begin{bmatrix}
n_1 & n_2
\end{bmatrix}
\vec v(\vec p)
$, while the transverse cone is $d a_2(\vec p) + \RR_+$.
Similarly, a linear functional $\vec \xi'$ projects onto
a linear functional $\xi = \sp {\xi'} n/d$ in the linear subspace.
Applying \eqref{eq:mu:1}, with $y = \frac{n_1}d D_1 + \frac{n_2}d D_2$
and $t = \fractional{- n_1 v_1(\vec p) - n_2 v_2(\vec p)}$, we therefore
find
\begin{align*}
D_{P(\vec p), F(\vec p)}
& =
-\sum_{n=0}^\infty
    \frac{b(n+1, \fractional{-n_1 v_1(\vec p) - n_2 v_2(\vec p)})}{(n+1)!}
        \left(\frac{n_1}d D_1 + \frac{n_2}d D_2\right)^n
\\
& =
- \sum_{i=0}^\infty \sum_{j=0}^\infty
    \frac{b(i+j+1, \fractional{-n_1 v_1(\vec p) - n_2 v_2(\vec p)})}{(i+j+1)!}
        \frac{n_1^i n_2^j}{d^{i+j}} D_1^i D_2^j
.
\end{align*}
After applying this differential operator to the polynomial
$h(\vec p)(\vec x)$, the resulting polynomial
$
h'(\vec p)(\vec x) = D_{P(\vec p), F(\vec p)} \cdot h(\vec p)(\vec x)
$
needs to be integrated over the facet.
The measure to be used is such that the integral of a lattice tile
in the linear space parallel to the facet is 1, i.e.,
$$
\int_{\vec 0}^{\vec f} 1 = \int_0^1 1 dz = 1,
$$
with $z$ the coordinate along $\vec f$.
Referring to \eqref{eq:EML:facet:coordinates} and
\eqref{eq:EML:facet:coordinates:2}, all points of the facet
have the form $\vec x(\vec p) = z \, \vec f + a_2(\vec p) \, \vec n$,
while the $z$-coordinate of the vertices $\vec v_1(\vec p)$
and $\vec v_2(\vec p)$ are
$(-n_2 v_{1,1} + n_1 v_{1,2})/d$
and
$(-n_2 v_{2,1} + n_1 v_{2,2})/d$, respectively.
That is, the contribution of the facet is equal to
$$
\int_{(-n_2 v_{1,1} + n_1 v_{1,2})/d}^{(-n_2 v_{2,1} + n_1 v_{2,2})/d}
h'(\vec p)\left(z \, \vec f + a_2(\vec p) \, \vec n\right) \, dz
,
$$
where, again, we need to ensure that the lower limit is smaller
than the upper limit using the usual method of plugging in a
particular value of the parameters.

Finally, we consider the contributions of the vertices.
The \ai{transverse cone}s are in this case simply the supporting cones.
Since $\mu$ is a valuation, we may apply \ai{Barvinok's decomposition}
and assume that the cone is unimodular.
For an affine cone
\begin{align*}
K &= \vec v(\vec p) + \RR_+ \vec r_1 + \RR_+ \vec r_2
\\
&= (a_1(\vec p) + \RR_+) \vec r_1 + (a_2(\vec p) + \RR_+) \vec r_2,
\end{align*}
with
$$
\vec a(\vec p) =
\begin{bmatrix}
\vec r_1 & \vec r_2
\end{bmatrix}^{-1}
\vec v(\vec p)
,
$$
we have~\cite[Proposition~31]{Berline2006local},
\begin{equation}
\label{eq:mu:2}
\mu(K)(\vec\xi)
=
\frac{e^{t_1 y_1 + t_2 y_2}}{(1-e^{y_1})(1-e^{y_2})}
+ \frac 1{y_1}B(y_2 - C_1 y_1, t_2)
+ \frac 1{y_2}B(y_1 - C_2 y_2, t_1)
- \frac 1{y_1 y_2},
\end{equation}
with
$$
B(y,t) =
\frac{e^{t y}}{1-e^y} + \frac 1{y}
= -\sum_{n=0}^\infty \frac{b(n+1, t)}{(n+1)!} y^n
,
$$
$y_i = \sps{\vec\xi}{\vec r_i}$,
$C_i = \sps{\vec v_1}{\vec v_2}/\sps{\vec v_i}{\vec v_i}$
and
$t_i = \fractional{-a_i(\vec p)}$.
Expanding \eqref{eq:mu:2}, we find
\begin{align*}
\mu(K)(\vec\xi)
& =
\left(
-\frac{b(0,t1)}{y_1} - \sum_{n=0}^\infty \frac{b(n+1,t_1)}{(n+1)!} y_1^n
\right)
\left(
-\frac{b(0,t2)}{y_2} - \sum_{n=0}^\infty \frac{b(n+1,t_2)}{(n+1)!} y_2^n
\right)
\\
& \phantom{=}
-
\left(
\sum_{n=0}^\infty \frac{b(n+1,t_2)}{(n+1)!} \frac{y_2^n}{y_1}
+
\sum_{n=0}^\infty \frac{b(n+1,t_2)}{(n+1)!} \frac{(y_2-C_1 y_1)^n-y_2^n}{y_1}
\right)
\\
& \phantom{=}
-
\left(
\sum_{n=0}^\infty \frac{b(n+1,t_1)}{(n+1)!} \frac{y_1^n}{y_2}
+
\sum_{n=0}^\infty \frac{b(n+1,t_1)}{(n+1)!} \frac{(y_1-C_2 y_2)^n-y_1^n}{y_2}
\right)
\\
& \phantom{=}
- \frac 1{y_1 y_2}
\\
& =
\sum_{n_1=0}^\infty
\sum_{n_2=0}^\infty
c(C_1, C_2, t_1, t_2; n_1, n_2) \, y_1^n y_2^n
,
\end{align*}
with
\begin{align*}
c(C_1, C_2, t_1, t_2; n_1, n_2)
= &
\frac{b(n_1+1,t_1)}{(n_1+1)!} \frac{b(n_2+1,t_2)}{(n_2+1)!}
\\
&
-
\frac{b(n_1+n_2+1,t_2)}{(n_1+n_2+1)!} {n_1+n_2+1 \choose n_1+1}
\left(-C_1\right)^{n_1+1}
\\
&
-
\frac{b(n_1+n_2+1,t_1)}{(n_1+n_2+1)!} {n_1+n_2+1 \choose n_2+1}
\left(-C_2\right)^{n_2+1}
.
\end{align*}
For $\vec \xi = \vec D$, we have
\begin{align*}
y_1^n y_2^n
& =
\left( r_{1,1} D_1 + r_{1,2} D_2 \right)^{n_1}
\left( r_{2,1} D_1 + r_{2,2} D_2 \right)^{n_2}
\\
& =
\left(
\sum_{k=0}^{n_1} r_{1,1}^k r_{1,2}^{n_1 - k} { n_1 \choose k} D_1^k D_2^{n_1-k}
\right)
\left(
\sum_{l=0}^{n_2} r_{2,1}^l r_{2,2}^{n_2 - l} { n_2 \choose l} D_1^l D_2^{n_2-l}
\right)
\end{align*}
and so
$
D_{P(\vec p), \vec v(\vec p)} = \mu(K)(\vec D)
=
$
$$
\sum_{i=0}^\infty
\sum_{j=0}^\infty
\sum_{\shortstack{$\scriptstyle i+j = n_1+n_2$\\$\scriptstyle n_1 \ge 0$\\$\scriptstyle n_2 \ge 0$}}
\sum_{\shortstack{$\scriptstyle k+l = i$\\$\scriptstyle 0 \le k \le n_1$\\$\scriptstyle 0 \le l \le n_2$}}
c(C_1, C_2, t_1, t_2; n_1, n_2)
r_{1,1}^k r_{1,2}^{n_1 - k}
r_{2,1}^l r_{2,2}^{n_2 - l}
{ n_1 \choose k} { n_2 \choose l} D_1^i D_2^j
.
$$
The contribution of this vertex is then
$$
h'(\vec p)(\vec v(\vec p))
,
$$
with $
h'(\vec p)(\vec x) = D_{P(\vec p), \vec v(\vec p)} \cdot h(\vec p)(\vec x)
$.

\begin{example}
As a simple example, consider the (non-parametric) triangle
in \autoref{f:EML:triangle} and assume we want to compute
$$
\sum_{\vec x \in T \cap \ZZ^2} x_1 x_2
.
$$
Since $T \cap \ZZ^2 = \{\, (2,4), (3,4), (2,5) \, \}$,
the result should be
$$
2 \cdot 4 + 3 \cdot 4 + 2 \cdot 5 = 30
.
$$

\begin{figure}
\intercol=1.2cm
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS@i@={(2,4),(3,4),(2,5),(2,4)},{0*[|(2)]\xypolyline{}}
\POS(2.35,4.25)*{x_1 x_2}
\POS(2,4)*+!U{(2,4)}
\POS(3,4)*+!U{(3,4)}
\POS(2,5)*+!D{(2,5)}
\POS(2,4)*{\cdot}
\POS(3,4)*{\cdot}
\POS(2,5)*{\cdot}
\POS(-1,0)\ar(4,0)
\POS(0,-1)\ar(0,5.5)
\end{xy}
\caption{Sum of polynomial $x_1 x_2$ over the integer points in a triangle $T$}
\label{f:EML:triangle}
\end{figure}

Let us first consider the integral
$$
\iint_T x_1 x_2 = \int_{\partial T} \frac{x_1^2 x_2}2 \, d x_2
.
$$
Integration along each of the edges of the triangle yields
the following.

\marginpar{%
\intercol=1cm
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(0,-1)*\xybox{
\POS@i@={(0,0),(1,0),(0,1),(0,0)},{0*\xypolyline{--}}
\POS(0,0)\ar@[|(2)](1,0)
}
\end{xy}
}%
For the edge in the margin, we have $\vec f = (1,0)$, i.e., $f_2 = 0$.
The contribution of this edge to the integral is therefore zero.

\marginpar{%
\intercol=1cm
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(0,-1)*\xybox{
\POS@i@={(0,0),(1,0),(0,1),(0,0)},{0*\xypolyline{--}}
\POS(1,0)\ar@[|(2)](0,1)
}
\end{xy}
}%
For this edge, we have $\vec f = (-1,1)$.
The contribution of this edge to the integral is therefore
$$
\int_0^1 \frac{(3-\lambda)^2(4+\lambda)}2 d\lambda
= \frac{337}{24}
.
$$

\marginpar{%
\intercol=1cm
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(0,-1)*\xybox{
\POS@i@={(0,0),(1,0),(0,1),(0,0)},{0*\xypolyline{--}}
\POS(0,1)\ar@[|(2)](0,0)
}
\end{xy}
}%
For this edge, we have $\vec f = (0,-1)$.
The contribution of this edge to the integral is therefore
$$
\int_0^1 \frac{2^2(5-\lambda)}2 (-1) d\lambda
= -9
.
$$

The total integral is therefore
$$
\int_{\partial T} \frac{x_1^2 x_2}2 \, d x_2
= 0 + \frac{337}{24} - 9 = \frac{121}{24}
.
$$

Now let us consider the contributions of the edges.
We will need the following \ai{Bernoulli number}s in our
computations.
\begin{align*}
b(1,0) & = - \frac 1 2
\\
b(2,0) & = \frac 1 6
\\
b(3,0) & = 0
\\
b(4,0) & = -\frac 1 {30}
\end{align*}

\marginpar{%
\intercol=1cm
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(0,-1)*\xybox{
\POS@i@={(0,0),(1,0),(0,1),(0,0)},{0*\xypolyline{--}}
\POS(0,0)\ar@[|(2)]@{-}(1,0)
\POS(0.5,0)\ar(0.5,1)
}
\end{xy}
}%
The normal to the facet $F_1$ in the margin is $\vec n = (0,1)$.
The vector $\vec f = (-1,0)$ is parallel to the facet.
We have
$$
\begin{bmatrix}
2 \\ 4
\end{bmatrix}
=
-2
\begin{bmatrix}
-1 \\ 0
\end{bmatrix}
+4
\begin{bmatrix}
0 \\ 1
\end{bmatrix}
\quad\text{and}\quad
\begin{bmatrix}
3 \\ 4
\end{bmatrix}
=
-3
\begin{bmatrix}
-1 \\ 0
\end{bmatrix}
+4
\begin{bmatrix}
0 \\ 1
\end{bmatrix}
.
$$
Therefore $t = \fractional{-4} = 0$, $y = D_2$,
\begin{align*}
D_{T,F_1}
& = 
- \sum_{j=0}^\infty \frac{b(j+1, 0)}{(j+1)!} D_2^j
\\
& =
- \frac{b(1,0)}1 - \frac{b(2,0)}2 D_2 + \cdots
\end{align*}
and
$$
h'(\vec x) =
D_{T,F_1} \cdot x_1 x_2 =
\left(\frac 1 2 - \frac 1{12} D_2\right) \cdot x_1 x_2
=
\frac 1 2 x_1 x_2 - \frac 1{12} x_1
.
$$
With $x_1 = - z$ and $x_2 = 4$, the contribution of this facet
is
$$
\int_{-3}^{-2} - 2 z + \frac 1{12} z \, dz
=
\frac{115}{24}
.
$$

\marginpar{%
\intercol=1cm
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(0,-1)*\xybox{
\POS@i@={(0,0),(1,0),(0,1),(0,0)},{0*\xypolyline{--}}
\POS(0,0)\ar@[|(2)]@{-}(0,1)
\POS(0,0.5)\ar(1,0.5)
}
\end{xy}
}%
The normal to the facet $F_2$ in the margin is $\vec n = (1,0)$.
The vector $\vec f = (0,1)$ is parallel to the facet.
We have
$$
\begin{bmatrix}
2 \\ 4
\end{bmatrix}
=
4
\begin{bmatrix}
0 \\ 1
\end{bmatrix}
+2
\begin{bmatrix}
1 \\ 0
\end{bmatrix}
\quad\text{and}\quad
\begin{bmatrix}
2 \\ 5
\end{bmatrix}
=
5
\begin{bmatrix}
0 \\ 1
\end{bmatrix}
+2
\begin{bmatrix}
1 \\ 0
\end{bmatrix}
.
$$
Therefore $t = \fractional{-2} = 0$, $y = D_1$,
\begin{align*}
D_{T,F_2}
& = 
- \sum_{i=0}^\infty \frac{b(i+1, 0)}{(i+1)!} D_1^i
\\
& =
- \frac{b(1,0)}1 - \frac{b(2,0)}2 D_1 + \cdots
\end{align*}
and
$$
h'(\vec x) =
D_{T,F_2} \cdot x_1 x_2 =
\left(\frac 1 2 - \frac 1{12} D_1\right) \cdot x_1 x_2
=
\frac 1 2 x_1 x_2 - \frac 1{12} x_2
.
$$
With $x_1 = 2$ and $x_2 = z$, the contribution of this facet
is
$$
\int_{4}^{5} z - \frac 1{12} z \, dz
=
\frac{33}{8}
.
$$

\marginpar{%
\intercol=1cm
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(0,-1)*\xybox{
\POS@i@={(0,0),(1,0),(0,1),(0,0)},{0*\xypolyline{--}}
\POS(1,0)\ar@[|(2)]@{-}(0,1)
\POS(0.5,0.5)\ar(-0.5,-0.5)
}
\end{xy}
}%
The normal to the facet $F_3$ in the margin is $\vec n = (-1,-1)$.
The vector $\vec f = (1,-1)$ is parallel to the facet.
We have
$$
\begin{bmatrix}
3 \\ 4
\end{bmatrix}
=
-\frac 1 2
\begin{bmatrix}
1 \\ -1
\end{bmatrix}
-\frac 7 2
\begin{bmatrix}
-1 \\ -1
\end{bmatrix}
\quad\text{and}\quad
\begin{bmatrix}
2 \\ 5
\end{bmatrix}
=
-\frac 3 2
\begin{bmatrix}
1 \\ -1
\end{bmatrix}
-\frac 7 2
\begin{bmatrix}
-1 \\ -1
\end{bmatrix}
.
$$
Therefore $t = \fractional{7} = 0$, $y = -\frac 1 2 D_1 -\frac 1 2 D_2$,
\begin{align*}
D_{T,F_3}
& = 
- \sum_{i=0}^\infty \sum_{j=0}^\infty
    \frac{b(i+j+1, 0)}{(i+j+1)!}
        \frac{(-1)^{i+j}}{2^{i+j}} D_1^i D_2^j
\\
& =
- \frac{b(1,0)}1
+ \frac 1 2 \frac{b(2,0)}2 D_1
+ \frac 1 2 \frac{b(2,0)}2 D_2 + \cdots
\end{align*}
and
$$
h'(\vec x) =
D_{T,F_4} \cdot x_1 x_2 =
\left(\frac 1 2 + \frac 1{24} D_1 + \frac 1{24} D_2\right) \cdot x_1 x_2
=
\frac 1 2 x_1 x_2 + \frac 1{24} x_2 + \frac 1{24} x_1
.
$$
With $x_1 = z + \frac 7 2$ and $x_2 = -z + \frac 7 2$,
the contribution of this facet
is
$$
\int_{-\frac 3 2}^{-\frac 1 2}
\frac 1 2 (z + \frac 7 2)(-z + \frac 7 2)
+ \frac 1{24}(-z + \frac 7 2)
+ \frac 1{24}(z + \frac 7 2) \, dz
=
\frac{47}{8}
.
$$

The total contribution of the edges is therefore
$$
\frac{115}{24}+\frac{33}8+
\frac{47}{8} = \frac{355}{24}
.
$$

Finally, we consider the contributions of the vertices.

\marginpar{%
\intercol=1cm
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(0,-1)*\xybox{
\POS@i@={(0,0),(1,0),(0,1),(0,0)},{0*\xypolyline{--}}
\POS(1,0)\ar@[|(2)](0,1)
\POS(1,0)\ar@[|(2)](0,0)
}
\end{xy}
}%
For the vertex $\vec v = (3,4)$, we have
$\vec r_1 = (-1,0)$ and $\vec r_2 = (-1,1)$.
Since $\vec v$ is integer, we have $t_1 = t_2 = 0$.
Also, $C_1 = 1$, $C_2 = 1/2$, $y_1 = -D_1$ and $y_2 = -D_1 + D_2$.
Since the total degree of the polynomial $x_1 x_2$ is two,
we only need the coefficients of $\mu(K)(\vec \xi)$ up to
$n_1+n_2 = 2$.

\noindent
\begin{tabular}{c|c|c}
$n_1$ & $n_2$
\\
\hline
0 & 0 &
$
\left(
\frac{b(1,0)}{1!}
\frac{b(1,0)}{1!}
-
\frac{b(2,0)}{2!}
{1 \choose 1}(-1)^1
-
\frac{b(2,0)}{2!}
{1 \choose 1}(-\frac 12)^1
\right)
$
\\
1 & 0 &
$
\left(
\frac{b(2,0)}{2!}
\frac{b(1,0)}{1!}
-
\frac{b(3,0)}{3!}
{2 \choose 2}(-1)^2
-
\frac{b(3,0)}{3!}
{2 \choose 1}(-\frac 12)^1
\right)
\left(
-D_1
\right)
$
\\
0 & 1 &
$
\left(
\frac{b(1,0)}{1!}
\frac{b(2,0)}{2!}
-
\frac{b(3,0)}{3!}
{2 \choose 1}(-1)^1
-
\frac{b(3,0)}{3!}
{2 \choose 2}(-\frac 12)^2
\right)
\left(
-D_1 + D_2
\right)
$
\\
2 & 0 &
$
\left(
\frac{b(3,0)}{3!}
\frac{b(1,0)}{1!}
-
\frac{b(4,0)}{4!}
{3 \choose 3}(-1)^3
-
\frac{b(4,0)}{4!}
{3 \choose 1}(-\frac 12)^1
\right)
\left(
-D_1
\right)^2
$
\\
1 & 1 &
$
\left(
\frac{b(2,0)}{2!}
\frac{b(2,0)}{2!}
-
\frac{b(4,0)}{4!}
{3 \choose 2}(-1)^2
-
\frac{b(4,0)}{4!}
{3 \choose 2}(-\frac 12)^2
\right)
\left(
-D_1
\right)
\left(
-D_1 + D_2
\right)
$
\\
0 & 2 &
$
\left(
\frac{b(1,0)}{1!}
\frac{b(3,0)}{3!}
-
\frac{b(4,0)}{4!}
{3 \choose 1}(-1)^1
-
\frac{b(4,0)}{4!}
{3 \choose 3}(-\frac 12)^3
\right)
\left(
-D_1 + D_2
\right)^2
$
\end{tabular}

We find
\begin{align*}
h'(\vec x)
&=
\left(
\frac 3 8 - \frac 1{24} (-D_1) - \frac 1{24} (-D_1 + D_2)
+ \frac 7{576} (-D_1 D_2)
- \frac 5{1152} (-2 D_1 D2)
\right) x_1 x_2
\\
&=
\frac 3 8 x_1 x_2 + \frac 1{24} x_2 - \frac 1{24} (-x_2 + x_1)
+ \frac 7{576} (-1)
- \frac 5{1152} (-2)
.
\end{align*}
The contribution of this vertex is therefore
$$
h'(3,4) = \frac {1355}{288}
.
$$

\marginpar{%
\intercol=1cm
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(0,-1)*\xybox{
\POS@i@={(0,0),(1,0),(0,1),(0,0)},{0*\xypolyline{--}}
\POS(0,1)\ar@[|(2)](1,0)
\POS(0,1)\ar@[|(2)](0,0)
}
\end{xy}
}%
For the vertex $\vec v = (2,5)$, we have
$\vec r_1 = (0,-1)$ and $\vec r_2 = (1,-1)$.
Since $\vec v$ is integer, we have $t_1 = t_2 = 0$.
Also, $C_1 = 1$, $C_2 = 1/2$, $y_1 = -D_2$ and $y_2 = D_1 - D_2$.
We similarly find
$$
h'(\vec x)
=
\frac 3 8 x_1 x_2 + \frac 1{24} x_1 - \frac 1{24} (x_2 - x_1)
+ \frac 7{576} (-1)
- \frac 5{1152} (-2)
.
$$
The contribution of this vertex is therefore
$$
h'(2,5) = \frac {1067}{288}
.
$$

\marginpar{%
\intercol=1cm
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(0,-1)*\xybox{
\POS@i@={(0,0),(1,0),(0,1),(0,0)},{0*\xypolyline{--}}
\POS(0,0)\ar@[|(2)](1,0)
\POS(0,0)\ar@[|(2)](0,1)
}
\end{xy}
}%
For the vertex $\vec v = (2,4)$, we have
$\vec r_1 = (1,0)$ and $\vec r_2 = (0,1)$.
Since $\vec v$ is integer, we have $t_1 = t_2 = 0$.
The computations are easier in this case since
$C_1 = C_2 = 0$, $y_1 = D_1$ and $y_2 = D_2$.
We find
$$
h'(\vec x)
=
\frac 1 4 x_1 x_2 - \frac 1{12} x_2 - \frac 1{12} x_1
+ \frac 1{144} (1)
.
$$
The contribution of this vertex is therefore
$$
h'(2,4) = \frac {253}{144}
.
$$

The total contribution of the vertices is then
$$
\frac {1355}{288} + \frac {1067}{288} + \frac {253}{144}
= \frac {61}6
$$
and the total sum is
$$
\frac{121}{24}+\frac{355}{24}+\frac{61}6 = 30
.
$$

\end{example}

\begin{example}
Consider the parametric polytope
$$
P(n) = \{\, \vec x \mid x_1 \ge 2 \wedge 3 x_1 \le n + 9
\wedge 4 \le x_2 \le 5 \,\}
.
$$
If $n \ge -3$, then the vertices of this polytope are
$(2,4)$, $(2,5)$, $(3+n/3,4)$ and $(3+n/3,5)$.
The contributions of the faces of $P(n)$ to
$$
\sum_{\vec x \in P(n) \cap \ZZ^2} x_1 x_2
$$
for the chamber $n \ge -3$ are shown in \autoref{t:sum:rectangle}.
The final result is
$$
\begin{cases}
\frac{ n^2}{2}
 - 3 n  \fractional{\frac{ n}{3}}
 + \frac{21}{2} n 
 + \frac{9}{2}  \fractional{\frac{ n}{3}}^2
 - \frac{63}{2}  \fractional{\frac{ n}{3}}
 + 45 
& \text{if $  n+3 \ge 0$}.
\end{cases}
$$

\begin{table}
\intercol=1cm
\begin{tabular}{lc}
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-1,-0.5)*\xybox{
\POS@i@={(0,0),(1,0),(1,1),(0,1),(0,0)},{0*[|(2)]\xypolyline{}}
\POS(0,1.1)*{}\POS(0,-0.1)*{}
}
\end{xy}
&
$
\displaystyle
\frac{ n^2}{4}
 + \frac{9}{2} n 
 + \frac{45}{4} 
$
\\
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-1,-0.5)*\xybox{
\POS@i@={(0,0),(1,0),(1,1),(0,1),(0,0)},{0*\xypolyline{--}}
\POS(0,0)\ar@[|(2)]@{-}(0,1)
\POS(0,0.5)*+!L{2}
\POS(0,1.1)*{}\POS(0,-0.1)*{}
}
\end{xy}
&
$
\displaystyle
\frac{33}{8}
$
\\
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-1,-0.5)*\xybox{
\POS@i@={(0,0),(1,0),(1,1),(0,1),(0,0)},{0*\xypolyline{--}}
\POS(1,0)\ar@[|(2)]@{-}(1,1)
\POS(1,0.5)*+!L{3+n/3}
\POS(0,1.1)*{}\POS(0,-0.1)*{}
}
\end{xy}
&
$
\displaystyle
 - \frac{3}{2} n  \fractional{\frac{ n}{3}}
 + \frac{3}{4} n 
 + \frac{9}{4}  \fractional{\frac{ n}{3}}^2
 - \frac{63}{4}  \fractional{\frac{ n}{3}}
 + \frac{57}{8} 
$
\\
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-1,-0.5)*\xybox{
\POS@i@={(0,0),(1,0),(1,1),(0,1),(0,0)},{0*\xypolyline{--}}
\POS(0,0)\ar@[|(2)]@{-}(1,0)
\POS(0.5,0)*+!D{4}
\POS(0,1.1)*{}\POS(0,-0.1)*{}
}
\end{xy}
&
$
\displaystyle
\frac{23}{216} n^2 
 + \frac{23}{12} n 
 + \frac{115}{24} 
$
\\
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-1,-0.5)*\xybox{
\POS@i@={(0,0),(1,0),(1,1),(0,1),(0,0)},{0*\xypolyline{--}}
\POS(0,1)\ar@[|(2)]@{-}(1,1)
\POS(0.5,1)*+!U{5}
\POS(0,1.1)*{}\POS(0,-0.1)*{}
}
\end{xy}
&
$
\displaystyle
\frac{31}{216} n^2 
 + \frac{31}{12} n 
 + \frac{155}{24} 
$
\\
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-1,-0.5)*\xybox{
\POS@i@={(0,0),(1,0),(1,1),(0,1),(0,0)},{0*\xypolyline{--}}
\POS(1,1)\ar@[|(2)](1,0)
\POS(1,1)\ar@[|(2)](0,1)
\POS(1,1)*+!LU{(3+n/3,5)}
\POS(0,1.1)*{}\POS(0,-0.1)*{}
}
\end{xy}
&
$
\displaystyle
 - \frac{31}{36} n  \fractional{\frac{ n}{3}}
 + \frac{31}{72} n 
 + \frac{31}{24}  \fractional{\frac{ n}{3}}^2
 - \frac{217}{24}  \fractional{\frac{ n}{3}}
 + \frac{589}{144} 
$
\\
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-1,-0.5)*\xybox{
\POS@i@={(0,0),(1,0),(1,1),(0,1),(0,0)},{0*\xypolyline{--}}
\POS(0,1)\ar@[|(2)](1,1)
\POS(0,1)\ar@[|(2)](0,0)
\POS(0,1)*+!LU{(2,5)}
\POS(0,1.1)*{}\POS(0,-0.1)*{}
}
\end{xy}
&
$
\displaystyle
\frac{341}{144}
$
\\
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-1,-0.5)*\xybox{
\POS@i@={(0,0),(1,0),(1,1),(0,1),(0,0)},{0*\xypolyline{--}}
\POS(0,0)\ar@[|(2)](1,0)
\POS(0,0)\ar@[|(2)](0,1)
\POS(0,0)*+!LD{(2,4)}
\POS(0,1.1)*{}\POS(0,-0.1)*{}
}
\end{xy}
&
$
\displaystyle
\frac{253}{144}
$
\\
\begin{xy}
<\intercol,0pt>:<0pt,\intercol>::
\POS(-1,-0.5)*\xybox{
\POS@i@={(0,0),(1,0),(1,1),(0,1),(0,0)},{0*\xypolyline{--}}
\POS(1,0)\ar@[|(2)](1,1)
\POS(1,0)\ar@[|(2)](0,0)
\POS(1,0)*+!LD{(3+n/3,4)}
\POS(0,1.1)*{}\POS(0,-0.1)*{}
}
\end{xy}
&
$
\displaystyle
 - \frac{23}{36} n  \fractional{\frac{ n}{3}}
 + \frac{23}{72} n 
 + \frac{23}{24}  \fractional{\frac{ n}{3}}^2
 - \frac{161}{24}  \fractional{\frac{ n}{3}}
 + \frac{437}{144} 
$
\end{tabular}
\caption{Contributions of the faces of $P(n)$ to the sum of $x_1 x_2$ over
the integer points of $P(n)$}
\label{t:sum:rectangle}
\end{table}

\end{example}

\subsection{Using TOPCOM to compute Chamber Decompositions}

In this section, we describe how to use the correspondence
between the \ai{regular triangulation}s of a point set
and the chambers of the \ai{Gale transform}
of the point set~\shortcite{Gelfand1994}
to compute the chamber decomposition of a parametric polytope.
This correspondence was also used by \shortciteN{Pfeifle2003}
\shortciteN{Eisenschmidt2007integrally}.

Let us first assume that the parametric polytope can be written as
\begin{equation}
\label{eq:TOPCOM:polytope}
\begin{cases}
    \begin{aligned}
\vec x &\ge 0
\\
A \, \vec x &\le \vec b(\vec p)
,
    \end{aligned}
\end{cases}
\end{equation}
where the right hand side $\vec b(\vec p)$ is arbitrary and
may depend on the parameters.
The first step is to add slack variables $\vec s$ to obtain
the \ai{vector partition} problem
$$
\begin{cases}
    \begin{aligned}
A \, \vec x + I \, \vec s & = \vec b(\vec p)
\\
\vec x, \vec s &\ge 0
,
    \end{aligned}
\end{cases}
$$
with $I$ the identity matrix.
Then we compute the (right) kernel $K$ of the matrix
$\begin{bmatrix}
A & I
\end{bmatrix}$, i.e.,
$$
\begin{bmatrix}
A & I
\end{bmatrix}
K
=
0
$$
and use \ai[\tt]{TOPCOM}'s \ai[\tt]{points2triangs} to
compute the \ai{regular triangulation}s of the points specified
by the rows of $K$.
Each of the resulting triangulations corresponds to a chamber
in the chamber complex of the above vector partition problem.
Each simplex in a triangulation corresponds to a parametric
vertex active on the corresponding chamber and
each point in the simplex (i.e., a row of $K$) corresponds
to a variable ($x_j$ or $s_j$) that is set to zero to obtain
this parametric vertex.
In the original formulation of the problem~\eqref{eq:TOPCOM:polytope}
each such variable set to zero reflects the saturation of the
corresponding constraint ($x_j = 0$ for $x_j = 0$ and
$\sps {\vec a_j}{\vec x} = b_j(\vec p)$ for $s_j = 0$).
A description of the chamber can then be obtained by plugging
in the parametric vertices in the remaining constraints.

\begin{example}
Consider the parametric polytope
$$
P(p,q,r) = \{\,
(i,j) \mid 0 \le i \le p \wedge
0 \le j \le 2 i + q \wedge
0 \le k \le i - p + r \wedge
p \ge 0 \wedge
q \ge 0 \wedge
r \ge 0
\,\}
.
$$
The constraints involving the variables are
$$
\begin{cases}
    \begin{aligned}
\begin{bmatrix}
1
\\
& 1
\\
& & 1
\end{bmatrix}
\begin{bmatrix}
i \\ j \\ k
\end{bmatrix}
&
\begin{matrix}
\ge
\\
\ge
\\
\ge
\end{matrix}
\begin{array}{l}
0
\\
0
\\
0
\end{array}
\\
\begin{bmatrix}
1 & 0 & 0
\\
-1 & 0 & 1
\\
-2 & 1 & 0
\end{bmatrix}
\begin{bmatrix}
i \\ j \\ k
\end{bmatrix}
&
\begin{matrix}
\le
\\
\le
\\
\le
\end{matrix}
\begin{array}{l}
p
\\
q
\\
-p + r
\end{array}
    \end{aligned}
\end{cases}
$$
We have
$$
\begin{bmatrix}
   1 &   0  &  0 &   1  &  0  &  0  \\
  -1 &   0  &  1 &   0  &  1  &  0  \\
  -2 &   1  &  0 &   0  &  0  &  1  \\
\end{bmatrix}
\begin{bmatrix}
  -1  &  0 &   0  \\
  -2  &  0 &  -1  \\
  -1  & -1 &   0  \\
   1  &  0 &   0  \\
   0  &  1 &   0  \\
   0  &  0 &   1  \\
\end{bmatrix}
= 0
$$

Computing the \ai{regular triangulation}s of the rows of $K$
using \ai[\tt]{TOPCOM}, we obtain
\begin{verbatim}
> cat e2.topcom
[
[  -1    0    0 ]
[  -2    0   -1 ]
[  -1   -1    0 ]
[   1    0    0 ]
[   0    1    0 ]
[   0    0    1 ]
]
> points2triangs --regular < e2.topcom 
T[1]:={{0,1,2},{1,2,3},{0,1,4},{1,3,4},{0,2,5},{2,3,5},{0,4,5},{3,4,5}};
T[2]:={{1,2,3},{1,3,4},{2,3,5},{3,4,5},{1,2,5},{1,4,5}};
T[3]:={{1,2,3},{1,3,4},{2,3,5},{3,4,5},{1,2,4},{2,4,5}};
\end{verbatim}

We see that we have three chambers in the decomposition,
one with 8 vertices and two with 6 vertices.
Take the second vertex (``\verb+{1,2,3}+'') of the first chamber.
This vertex corresponds
to the saturation of the constraints $j \ge 0$, $k \ge 0$
and $i \le p$, i.e., $(i,j,k) = (p,0,0)$.  Plugging in this
vertex in the remaining constraints, we see that it is only valid
in case $p \ge 0$, $r \ge 0$ and $2p + q \ge 0$.
For the remaining vertices of the first chamber, we similarly find
\\
\begin{tabular}{ccc}
% e0
\verb+{0,1,2}+ & $(0,0,0)$ & $p \ge 0$, $-q + r \ge 0$ and $q \ge 0$
\\
% 70
\verb+{1,2,3}+ & $(p,0,0)$ & $p \ge 0$, $r \ge 0$ and $2p + q \ge 0$
\\
% c8
\verb+{0,1,4}+ & $(0,0,-p+r)$ & $-q + r \ge 0$, $p \ge 0$ and $q \ge 0$
\\
% 58
\verb+{1,3,4}+ & $(p,0,r)$ & $p \ge 0$, $r \ge 0$ and $2p + q \ge 0$
\\
% a4
\verb+{0,2,5}+ & $(0,q,0)$ & $q \ge 0$, $p \ge 0$ and $-q + r \ge 0$
\\
% 34
\verb+{2,3,5}+ & $(p, 2p+q, 0)$ & $p \ge 0$, $2p + q \ge 0$ and $r \ge 0$
\\
% 8c
\verb+{0,4,5}+ & $(0, q, -p+r)$ & $q \ge 0$, $-q + r \ge 0$ and $p \ge 0$
\\
% 1c
\verb+{3,4,5}+ & $(p, 2p+q, r)$ & $p \ge 0$, $2p + q \ge 0$ and $r \ge 0$
\end{tabular}
\\
Combining these constraints with the initial constraints of the problem
on the parameters
$p \ge 0$, $q \ge 0$ and $r \ge 0$, we find the chamber
$
\{\,
(p,q,r) \mid p \ge 0 \wedge -p + r \ge 0 \wedge q \ge 0
\,\}
$.
For the second chamber, we have
\\
\begin{tabular}{ccc}
% 70
\verb+{1,2,3}+ & $(p,0,0)$ & $p \ge 0$, $r \ge 0$ and $2p + q \ge 0$
\\
% 58
\verb+{1,3,4}+ & $(p,0,r)$ & $p \ge 0$, $r \ge 0$ and $2p + q \ge 0$
\\
% 34
\verb+{2,3,5}+ & $(p, 2p+q, 0)$ & $p \ge 0$, $2p + q \ge 0$ and $r \ge 0$
\\
% 1c
\verb+{3,4,5}+ & $(p, 2p+q, r)$ & $p \ge 0$, $2p + q \ge 0$ and $r \ge 0$
\\
% 64
\verb+{1,2,5}+ & $(-\frac q 2,0,0)$ &
        $-q \ge 0$, $2p + q \ge 0$ and $-2p -q+2r \ge 0$
\\
% 4c
\verb+{1,4,5}+ & $(-\frac q 2,0,-p-\frac q 2+r)$ &
        $-q \ge 0$, $-2p -q+2r \ge 0$ and $2p + q \ge 0$
\end{tabular}
\\
The chamber is therefore
$
\{\,
(p,q,r) \mid q = 0 \wedge p \ge 0 \wedge -p +r \ge 0
\,\}
$.
Note that by intersecting with the initial constraints this chamber
is no longer full-dimensional and can therefore be discarded.
Finally, for the third chamber, we have
\\
\begin{tabular}{ccc}
% 70
\verb+{1,2,3}+ & $(p,0,0)$ & $p \ge 0$, $r \ge 0$ and $2p + q \ge 0$
\\
% 58
\verb+{1,3,4}+ & $(p,0,r)$ & $p \ge 0$, $r \ge 0$ and $2p + q \ge 0$
\\
% 34
\verb+{2,3,5}+ & $(p, 2p+q, 0)$ & $p \ge 0$, $2p + q \ge 0$ and $r \ge 0$
\\
% 1c
\verb+{3,4,5}+ & $(p, 2p+q, r)$ & $p \ge 0$, $2p + q \ge 0$ and $r \ge 0$
\\
% 68
\verb+{1,2,4}+ & $(p-r,0,0)$ &
        $p -r \ge 0$, $r \ge 0$ and $2p +q -2r \ge 0$
\\
% 2c
\verb+{2,4,5}+ & $(p-r,2p+q-2r, 0)$ &
        $p -r \ge 0$, $2p +q -2r \ge 0$ and $r \ge 0$
\end{tabular}
\\
The chamber is therefore
$
\{\,
(p,q,r) \mid p - r \ge 0 \wedge q \ge 0 \wedge r \ge 0
\,\}
$.
\end{example}

Now let us consider general parametric polytopes.
First note that we can follow the same procedure as above
if we replace $\vec x$ by $\vec x' - \vec c(\vec p)$
in \eqref{eq:TOPCOM:polytope}, i.e.,
if our problem has the form
\begin{equation}
\label{eq:TOPCOM:polytope:2}
\begin{cases}
    \begin{aligned}
\vec x' &\ge \vec c(\vec p)
\\
A \, \vec x' &\le \vec b(\vec p) + A \vec c(\vec p)
,
    \end{aligned}
\end{cases}
\end{equation}
as saturating a constraint $x_i = 0$ is equivalent
to saturating the constraint $x_i' = c_i(\vec p)$
and, similarly, $\sps {\vec a_j}{\vec x} = b_j(\vec p)$
is equivalent to
$\sps {\vec a_j}{\vec x'} = b_j(\vec p) + \sps {\vec a_j}{\vec c(\vec p)}$.

In the general case, the problem has the form
$$
A \vec x \ge \vec b(\vec p)
.
$$
Let $A'$ be a non-singular square submatrix of $A$ with the same number
of columns and compute the (left) \indac{HNF} $H = A' U$ with $U$ unimodular
and $H$ lower-triangular with non-positive elements below the diagonal.
Replacing $\vec x$ by $U \vec x'$, we obtain
$$
\begin{cases}
    \begin{aligned}
H \vec x' &\ge \vec b'(\vec p)
\\
-A''U \, \vec x' &\le -\vec b''(\vec p)
,
    \end{aligned}
\end{cases}
$$
with $A''$ the remaining rows of $A$ and $\vec b(\vec p)$ split
in the same way.
If $H$ happens to be the identity matrix, then our problem is
of the form \eqref{eq:TOPCOM:polytope:2} and we already know how
to solve this problem.
Note that, again, saturating any of the transformed constraints
in $\vec x'$ is equivalent to saturating the corresponding constraint
in $\vec x$.  We therefore only need to compute $-A'' U$ for the
computation of the kernel $K$.  To construct the parametric vertices
in the original coordinate system, we can simply use the original
constraints.
The same reasoning holds if $H$ is any diagonal matrix, since
we can virtually replace $H \vec x$ by $\vec x'$ without affecting
the non-negativity of the variables.

If $H$ is not diagonal, then we can introduce new constraints
$x'_j \ge d(\vec p)$, where $d(\vec p)$ is some symbolic constant.
These constraints do not remove any solutions
since each row in $H$ expresses that the corresponding variable is
greater than or equal to a non-negative combination of the
previous variables plus some constant.
We can then proceed as before.  However, to reduce unnecessary computations
we may remove from $K$ the rows that correspond to these new rows.
Any solution saturating the new constraint, would also saturate
the corresponding constraint $\vec h_j^\T$ and all
the constraints corresponding to the non-zero
entries in $\vec h_j^\T$.
If a chamber contains a vertex obtained by saturating such a new
constraint, it would appear multiple times in the same chamber,
each time combined with different constraints from the above set.
Furthermore, there would also be another (as it turns out, identical)
chamber where the vertex is only defined by the other constraints.

\begin{example}
Consider the parametric polytope
$$
P(n) = \{\,
(i,j) \mid
1 \le i \wedge 2 i \le 3 j \wedge j \le n
\,\}
.
$$
The constraints are
$$
\begin{bmatrix}
1 & 0 \\
-2 & 3 \\
0 & -1
\end{bmatrix}
\begin{bmatrix}
i \\ j
\end{bmatrix}
\ge
\begin{bmatrix}
1 \\
0 \\
-n
\end{bmatrix}
.
$$
The top $2 \times 2$ submatrix is already in \indac{HNF}.
We have $3 j \ge 2i \ge 2$, so we can add a constraint
of the form $j \ge c(n)$ and obtain
$$
\begin{bmatrix}
A & I
\end{bmatrix}
=
\begin{bmatrix}
0 & 1 & 1 & 0
\\
2 & -3 & 0 & 1
\end{bmatrix}
,
$$
while $K$ with $\begin{bmatrix}A & I\end{bmatrix} K = 0$ is given
by
$$
\begin{bmatrix}
0 & 1 & 1 & 0
\\
2 & -3 & 0 & 1
\end{bmatrix}
\begin{bmatrix}
1 & 0 \\
0 & 1 \\
0 & -1 \\
-2 & 3
\end{bmatrix}
.
$$
The second row of $K$ corresponds to the second variable,
which in turn corresponds to the newly added constraint.
Passing all rows of $K$ to \ai[\tt]{TOPCOM} we would get
\begin{verbatim}
> points2triangs --regular <<EOF
> [[1 0],[0,1],[0,-1],[-2,3]]
> EOF
T[1]:={{0,1},{0,2},{1,3},{2,3}};
T[2]:={{0,2},{2,3},{0,3}};
T[3]:={};
\end{verbatim}
The first vertex in the first chamber saturates the second row
(row 1) and therefore saturates both the first (0) and fourth (3)
and it appears a second time as \verb+{1,3}+.  Combining
these ``two'' vertices into one as \verb+{0,3}+ results in the
second (identical) chamber.
Removing the row corresponding to the new constraint from $K$
we remove the duplicates
\begin{verbatim}
> points2triangs --regular <<EOF
> [[1 0],[0,-1],[-2,3]]
> EOF
T[1]:={{0,1},{1,2},{0,2}};
T[2]:={};
\end{verbatim}
Note that in this example, we also could have interchanged
the second and the third constraint and then have replaced $j$ by $-j'$.
\end{example}

In practice, this method of computing a \ai{chamber decomposition}
does not seem to perform very well, mostly because
\ai[\tt]{TOPCOM} can not exploit all available information
about the parametric polytopes and will therefore compute
many redundant triangulations/chambers.
In particular, any chamber that does not intersect with
the parameter domain of the parametric polytope, or only
intersects in a face of this parameter domain, is completely redundant.
Furthermore, if the parametric polytope is not simple, then many
different combinations of the constraints will lead to the same parametric
vertex.  Many triangulations will therefore correspond to one and the
same chamber in the chamber complex of the parametric polytope.
For example, for a dilated octahedron, \ai[\tt]{TOPCOM} will
compute 150 triangulations/chambers, 104 of which are empty,
while the remaining 46 refer to the same single chamber.