barvinok_maximize: move --(bernstein-)recurse option to main library

[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,
\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}