evalue.c: add evalue_add_constant

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

\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 decription 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.