doc/implementation.tex

   1 \section{Implementation details}
   2
   3 \subsection{The integer points in the fundamental parallelepiped of a simple cone}
   4
   5 \label{s:fundamental}
   6
   7 This section is based on \shortciteN[Lemma 5.1]{Barvinok1992volume} and
   8 \shortciteN{Koeppe2006experiments}.
   9
  10 \sindex{simple}{cone}
  11 \sindex{open}{facet}
  12 \sindex{open}{ray}
  13 \sindex{explicit}{representation}
  14 In this section we will deal exclusively with \ai{simple cone}s,
  15 i.e. $d$-dimensional cones with $d$ extremal rays and $d$ facets.
  16 \index{open facet}%
  17 Some of the facets of these cones may be open.
  18 Since we will mostly be dealing with cones in their
  19 \ai{explicit representation}, we will have occasion to speak of
  20 ``\ai{open ray}s'', by which we will mean that the facet not
  21 containing the ray is open.  (There is only one such facet because the cone
  22 is simple.)
  23
  24 \sindex{fundamental}{parallelepiped}
  25 \begin{definition}[Fundamental parallelepiped]
  26 Let $K = \vec v + \poshull \lb\, \vec u_i \,\rb$ be
  27 a closed (shifted) cone, then the \defindex{fundamental parallelepiped} $\Pi$
  28 of $K$ is
  29 $$
  30 \Pi = \vec v +
  31 \lb\, \sum_i \alpha_i \vec u_i \mid 0 \leq \alpha_i < 1 \,\rb
  32 .
  33 $$
  34 If some of the rays $\vec u_i$ of $K$ are open, then the constraints on
  35 the corresponding coefficient $\alpha_i$ are such that $0 < \alpha_i \le 1$.
  36 \end{definition}
  37
  38 \begin{lemma}[Integer points in the fundamental parallelepiped of a simple cone]
  39 Let $K = \vec v + \poshull \lb\, \vec u_i \,\rb$ be a closed simple cone
  40 and let $A$ be the matrix with the generators $\vec u_i$ of $K$
  41 as rows.
  42 Furthermore let $V A W^{-1} = S = \diag \vec s$ be the \indac{SNF} of $A$.
  43 Then the integer points in the fundamental parallelepiped of $K$ are given
  44 by
  45 \begin{eqnarray}
  46 \label{eq:parallelepiped}
  47 \vec w^\T & = & \vec v^\T + \fractional{(\vec k^\T W - \vec v^\T) A^{-1}} A
  48 \\
  49 \nonumber
  50 & = &
  51 \vec v^T +
  52 \sum_{i=1}^d
  53     \fractional{\sps{\sum_{j=1}^d k_j \vec w^T_j - \vec v^\T}{\vec u^*_i}} \vec u_i,
  54 \end{eqnarray}
  55 where $\vec u^*_i$ are the columns of $A^{-1}$ and $k_j \in \ZZ$ ranges
  56 over $0 \le k_j < s_j$.
  57 \end{lemma}
  58
  59 \begin{proof}
  60 Since $0 \le \fractional{x} < 1$, it is clear that each such $\vec w$
  61 lies inside the fundamental parallelepiped.
  62 Furthermore,
  63 \begin{eqnarray*}
  64 \vec w^\T & = & \vec v^\T + \fractional{(\vec k^\T W - \vec v^\T) A^{-1}} A
  65 \\
  66 & = &
  67 \vec v^T +
  68 \left(
  69 (\vec k^\T W - \vec v^\T) A^{-1} - \floor{(\vec k^\T W - \vec v^\T) A^{-1}}
  70 \right) A
  71 \\
  72 & = &
  73 \underbrace{\vec k^\T W\mathstrut}_{\in \ZZ^{1\times d}}
  74 +
  75 \underbrace{\floor{(\vec k^\T W - \vec v^\T) A^{-1}}}_{\in \ZZ^{1\times d}}
  76 \underbrace{A\mathstrut}_{\in \ZZ^{d\times d}} \in \ZZ^{1\times d}.
  77 \end{eqnarray*}
  78 Finally, if two such $\vec w$ are equal, i.e., $\vec w_1 = \vec w_2$,
  79 then
  80 \begin{eqnarray*}
  81 \vec 0^\T = \vec w_1^\T - \vec w_2^\T
  82 & = & \vec k_1^\T W - \vec k_2^\T W + \vec p^\T A
  83 \\
  84 & = & \left(\vec k_1^\T - \vec k_2^\T \right) W + \vec p^\T V^{-1} S W,
  85 \end{eqnarray*}
  86 with $\vec p \in \ZZ^d$,
  87 or $\vec k_1 \equiv \vec k_2 \mod \vec s$, i.e., $\vec k_1 = \vec k_2$.
  88 Since $\det S = \det A$, we obtain all points in the fundamental parallelepiped
  89 by taking all $\vec k \in \ZZ^d$ satisfying $0 \le k_j < s_j$.
  90 \end{proof}
  91
  92 If the cone $K$ is not closed then the coefficients of the open rays
  93 should be in $(0,1]$ rather than in $[0,1)$.
  94 In (\ref{eq:parallelepiped}),
  95 we therefore need to replace the fractional part $\fractional{x} = x - \floor{x}$
  96 by $\cractional{x} = x - \ceil{x-1}$ for the open rays.
  97
  98 \begin{figure}
  99 \intercol=1.2cm
 100 \begin{xy}
 101 <\intercol,0pt>:<0pt,\intercol>::
 102 \POS@i@={(0,-3),(0,0),(4,2),(4,-3)},{0*[grey]\xypolyline{*}}
 103 \POS@i@={(0,-3),(0,0),(4,2)},{0*[|(2)]\xypolyline{}}
 104 \POS(-1,0)\ar(4.5,0)
 105 \POS(0,-3)\ar(0,3)
 106 \POS(0,0)\ar@[|(3)](0,-1)
 107 \POS(0,0)\ar@[|(3)](2,1)
 108 \POS(0,-1)\ar@{--}@[|(2)](2,0)
 109 \POS(2,1)\ar@{--}@[|(2)](2,0)
 110 \POS(0,0)*{\bullet}
 111 \POS(1,0)*{\bullet}
 112 \end{xy}
 113 \caption{The integer points in the fundamental parallelepiped of $K$}
 114 \label{f:parallelepiped}
 115 \end{figure}
 116
 117 \begin{example}
 118 Let $K$ be the cone
 119 $$
 120 K = \sm{0 \\ 0} + \poshull \lb\, \sm{2 \\ 1}, \sm{0 \\ -1} \,\rb
 121 ,
 122 $$
 123 shown in Figure~\ref{f:parallelepiped}.
 124 Then
 125 $$
 126 A = \sm{2 & 1\\0 & -1} \qquad A^{-1} = \sm{1/2 & 1/2 \\ 0 & -1 }
 127 $$
 128 and
 129 $$
 130 \sm{1 & 0 \\ 1 & 1 } \sm{2 & 1\\0 & -1} = \sm{1 & 0 \\ 0 & 2} \sm{2 & 1 \\ 1 & 0}.
 131 $$
 132 We have $\det A = \det S = 2$ and
 133 $\vec k_1^\T = \sm{0 & 0}$ and $\vec k_2^\T = \sm{0 & 1}$.
 134 Therefore,
 135 $$
 136 \vec w_1^\T = \fractional{\sm{0 & 0} \sm{2 & 1 \\ 1 & 0} \sm{1/2 & 1/2 \\ 0 & -1 }}
 137 \sm{2 & 1\\0 & -1} = \sm{0 & 0}
 138 $$
 139 and
 140 \begin{eqnarray*}
 141 \vec w_2^\T & = &
 142 \fractional{\sm{0 & 1} \sm{2 & 1 \\ 1 & 0} \sm{1/2 & 1/2 \\ 0 & -1 }}
 143 \sm{2 & 1\\0 & -1}
 144 \\
 145 & = &
 146 \sm{1/2 & 1/2} \sm{2 & 1\\0 & -1} = \sm{1 & 0}.
 147 \end{eqnarray*}
 148 \end{example}
 149
 150
 151
 152
 153 \subsection{Barvinok's decomposition of simple cones in primal space}
 154
 155 As described by \shortciteN{DeLoera2003effective}, the first
 156 implementation of Barvinok's counting algorithm applied
 157 \ai{Barvinok's decomposition} \shortcite{Barvinok1994} in the \ai{dual space}.
 158 \ai{Brion's polarization trick} \shortcite{Brion88} then ensures that you
 159 do not need to worry about lower-dimensional faces in the decomposition.
 160 Another way of avoiding the lower-dimensional faces, in the \ai{primal space},
 161 is to perturb the vertex of the cone such that none of the lower-dimensional
 162 face encountered contain any integer points \shortcite{Koeppe2006primal}.
 163 In this section, we describe another technique that is based on allowing
 164 some of the facets of the cone to be open.
 165
 166 The basic step in Barvinok's decomposition is to replace a
 167 $d$-dimensional simple cone
 168 $K = \poshull \lb\, \vec u_i \,\rb_{i=1}^d \subset \QQ^d$
 169 by a signed sum of (at most) $d$ cones $K_j$
 170 with a smaller determinant (in absolute value).
 171 The cones are obtained by successively replacing each generator
 172 of $K$ by an appropriately chosen
 173 $\vec w = \sum_{i=1}^d \alpha_i \vec u_i$, i.e.,
 174 \begin{equation}
 175 \label{eq:K_j}
 176 K_j =
 177 \poshull \left(\lb\, \vec u_i \,\rb_{i=1}^d
 178 \setminus \{\, \vec u_j \,\} \cup \{\, \vec w \,\}\right)
 179 .
 180 \end{equation}
 181 To see that we can use these $K_j$ to perform a decomposition,
 182 rearrange the $\vec u_i$ such that for all $1 \le i \le k$ we have
 183 $\alpha_i < 0$ and for all $k+1 \le i \le d'$ we have $\alpha_i > 0$,
 184 with $d - d'$ the number of zero $\alpha_i$.
 185 We may assume $k < d'$; otherwise replace $\vec w \in B$ by
 186 $-\vec w \in B$.  We have
 187 $$
 188 \vec w + \sum_{i=1}^k (-\alpha_i) \vec u_i =
 189 \sum_{i=k+1}^{d'} \alpha_i \vec u_i
 190 $$
 191 or
 192 \begin{equation}
 193 \label{eq:sub}
 194 \sum_{i=0}^k \beta_i \vec u_i =
 195 \sum_{i=k+1}^{d'} \alpha_i \vec u_i
 196 ,
 197 \end{equation}
 198 with $\vec u_0 = \vec w$, $\beta_0 = 1$ and $\beta_i = -\alpha_i > 0$
 199 for $1 \le i \le k$.  Any two $\vec u_j$ and $\vec u_l$ on the same side
 200 of the equality are on opposite sides of the linear hull $H$ of
 201 the other $\vec u_i$s since there exists a convex combination
 202 of $\vec u_j$ and $\vec u_l$ on this hyperplane.
 203 In particular, since $\alpha_j$ and $\alpha_l$ have the same sign,
 204 we have
 205 \begin{equation}
 206 \label{eq:opposite}
 207 \frac {\alpha_j}{\alpha_j+\alpha_l} \vec u_j
 208 +
 209 \frac {\alpha_l}{\alpha_j+\alpha_l} \vec u_l
 210 \in H
 211 \qquad\text{for $\alpha_i \alpha_l > 0$}
 212 .
 213 \end{equation}
 214 The corresponding cones $K_j$ and $K_l$ (with $K_0 = K$)
 215 therefore intersect in a common face $F \subset H$.
 216 Let
 217 $$
 218 K' :=
 219 \poshull \left(\lb\, \vec u_i \,\rb_{i=1}^d \cup \{\, \vec w \,\}\right)
 220 ,
 221 $$
 222 then any $\vec x \in K'$ lies both in some cone $K_i$ with
 223 $0 \le i \le k$ and in some cone $K_i$ with $k+1 \le i \le d'$.
 224 (Just subtract an appropriate multiple of Equation~(\ref{eq:sub}).)
 225 The cones
 226 $\{\, K_i \,\}_{i=0}^k$
 227 and
 228 $\{\, K_i \,\}_{i=k+1}^{d'}$
 229 therefore both form a triangulation of $K'$ and hence
 230 \begin{equation}
 231 \label{eq:triangulations}
 232 \indf{K'}
 233 =
 234 \indf{K} + \sum_{i=1}^k \indf{K_i} - \sum_{j\in J_1} \indf{F_j}
 235 =
 236 \sum_{i=k+1}^{d'} \indf{K_i} - \sum_{j\in J_2} \indf{F_j}
 237 \end{equation}
 238 or
 239 \begin{equation}
 240 \label{eq:decomposition}
 241 \indf{K} = \sum_{i=1}^{d'} \varepsilon_i \indf{K_i} + \sum_j \delta_j \indf{F_j}
 242 ,
 243 \end{equation}
 244 with $\varepsilon_i = -1$ for $1 \le i \le k$,
 245 $\varepsilon_i = 1$ for $k+1 \le i \le d'$,
 246 $\delta_j \in \{ -1, 1 \}$ and $F_j$ some lower-dimensional faces.
 247 Figure~\ref{fig:w} shows the possible configurations
 248 in the case of a $3$-dimensional cone.
 249
 250 \begin{figure}
 251 \intercol=0.48cm
 252 \begin{center}
 253 \begin{minipage}{0cm}
 254 \begin{xy}
 255 <\intercol,0pt>:<0pt,\intercol>::
 256 *
 257 \xybox{
 258 \POS(-2,-1)="a"*+!U{+}
 259 \POS(2,0)="b"*+!L{+}
 260 \POS(0,2)="c"*+!D{+}
 261 \POS(0,0)="w"*+!DR{\vec w}
 262 \POS"a"\ar@{-}"b"
 263 \POS"b"\ar@{-}"c"
 264 \POS"c"\ar@{-}"a"
 265 \POS"a"\ar@{--}"w"
 266 \POS"b"\ar@{--}"w"
 267 \POS"c"\ar@{--}"w"
 268 }="a"
 269 +R+(2,0)*!L
 270 \xybox{
 271 \POS(-2,-1)="a"*+!U{+}
 272 \POS(2,0)="b"*+!L{-}
 273 \POS(0,2)="c"*+!D{+}
 274 \POS(-3,1)="w"*+!DR{\vec w}
 275 \POS"a"\ar@{-}"b"
 276 \POS"b"\ar@{-}"c"
 277 \POS"c"\ar@{-}"a"
 278 \POS"a"\ar@{--}"w"
 279 \POS"b"\ar@{--}"w"
 280 \POS"c"\ar@{--}"w"
 281 }="b"
 282 +R+(2,0)*!L
 283 \xybox{
 284 \POS(-2,-1)="a"*+!U{-}
 285 \POS(2,0)="b"*+!U{+}
 286 \POS(0,2)="c"*+!D{-}
 287 \POS(5,-1)="w"*+!L{\vec w}
 288 \POS"a"\ar@{-}"b"
 289 \POS"b"\ar@{-}"c"
 290 \POS"c"\ar@{-}"a"
 291 \POS"a"\ar@{--}"w"
 292 \POS"b"\ar@{--}"w"
 293 \POS"c"\ar@{--}"w"
 294 }
 295 \POS"a"
 296 +D-(0,1)*!U
 297 \xybox{
 298 \POS(-2,-1)="a"*+!U{0}
 299 \POS(2,0)="b"*+!L{+}
 300 \POS(0,2)="c"*+!D{+}
 301 \POS(1,1)="w"*+!DL{\vec w}
 302 \POS"a"\ar@{-}"b"
 303 \POS"b"\ar@{-}"c"
 304 \POS"c"\ar@{-}"a"
 305 \POS"a"\ar@{--}"w"
 306 }
 307 \POS"b"
 308 +DL-(0,1)*!UL
 309 \xybox{
 310 \POS(-2,-1)="a"*+!U{0}
 311 \POS(2,0)="b"*+!U{+}
 312 \POS(0,2)="c"*+!D{-}
 313 \POS(4,-2)="w"*+!L{\vec w}
 314 \POS"a"\ar@{-}"b"
 315 \POS"b"\ar@{-}"c"
 316 \POS"c"\ar@{-}"a"
 317 \POS"a"\ar@{--}"w"
 318 \POS"b"\ar@{--}"w"
 319 }
 320 \end{xy}
 321 \end{minipage}
 322 \end{center}
 323 \caption[Possible locations of the vector $\vec w$ with respect to the rays
 324 of a $3$-dimensional cone.]
 325 {Possible locations of $\vec w$ with respect to the rays
 326 of a $3$-dimensional cone.  The figure shows a section of the cones.}
 327 \label{fig:w}
 328 \end{figure}
 329
 330 As explained above there are several ways of avoiding the lower-dimensional
 331 faces in (\ref{eq:decomposition}).  Here we will apply the following proposition.
 332 \begin{proposition}[\shortciteN{Koeppe2007parametric}]
 333   Let
 334   \begin{equation}
 335     \label{eq:full-source-identity}
 336     \sum_{i\in {I_1}} \epsilon_i [P_i] + \sum_{i\in {I_2}} \delta_k [P_i] = 0
 337   \end{equation}
 338   be a (finite) linear identity of indicator functions of closed
 339   polyhedra~$P_i\subseteq\QQ^d$, where the
 340   polyhedra~$P_i$ with $i \in I_1$ are full-dimensional and those with $i \in I_2$
 341   lower-dimensional.  Let each closed polyhedron be given as
 342 $$
 343     P_i = \left\{\, \vec x \mid \sp{b^*_{i,j}}{x} \ge \beta_{i,j} \text{
 344       for $j\in J_i$}\,\right\}
 345   .
 346 $$
 347   Let $\vec y\in\QQ^d$ be a vector such that $\langle \vec b^*_{i,j}, \vec
 348   y\rangle \neq 0$ for all $i\in I_1\cup I_2$, $j\in J_i$.
 349   For each $i\in I_1$, we define the half-open polyhedron
 350   \begin{equation}
 351     \label{eq:half-open-by-y}
 352     \begin{aligned}
 353       \tilde P_i = \Bigl\{\, \vec x\in\QQ^d \mid {}&
 354             \sp{b^*_{i,j}}{x} \ge \beta_{i,j}
 355       \text{ for $j\in J_i$ with $\sp{b^*_{i,j}}{y} > 0$,} \\
 356       & \sp{b^*_{i,j}}{x} > \beta_{i,j}
 357       \text{ for $j\in J_i$ with $\sp{b^*_{i,j}}{y} < 0$} \,\Bigr\}.
 358     \end{aligned}
 359   \end{equation}
 360   Then
 361   \begin{equation}
 362     \label{eq:target-identity}
 363     \sum_{i\in I_1} \epsilon_i [\tilde P_i] = 0.
 364   \end{equation}
 365 \end{proposition}
 366 When applying this proposition to (\ref{eq:decomposition}), we obtain
 367 \begin{equation}
 368 \label{eq:decomposition:2}
 369 \indf{\tilde K} = \sum_{i=1}^{d'} \varepsilon_i \indf{\tilde K_i}
 370 ,
 371 \end{equation}
 372 where we start out
 373 from a given $\tilde K$, which may be $K$ itself, i.e., a fully closed cone,
 374 or the result of a previous application of the proposition.
 375 In either case, a suitable $\vec y$ is available, either as an interior
 376 point of the cone or as the vector used in the previous application
 377 (which may require a slight perturbation if it happens to lie on one of
 378 the new facets of the cones $K_i$).
 379 We are, however, free to construct a new $\vec y$ on each application
 380 of the proposition.
 381 In fact, we will not even construct such a vector explicitly, but
 382 rather apply a set of rules that is equivalent to a valid choice of $\vec y$.
 383
 384 The vector $\vec y$ has to satisfy $\sp{b^*_j}y > 0$ for normals $\vec b^*_j$
 385 of closed facets and $\sp{b^*_j}y < 0$ for normals $\vec b^*_j$ of open facets of
 386 $\tilde K$.
 387 These constraints delineate a non-empty open cone $R$ from which
 388 $\vec y$ should be selected.  For some of the new facets of the cones
 389 $\tilde K_j$, the cone $R$ will not be cut by the affine hull of the facet.
 390 The closedness of these facets is therefore predetermined by $\tilde K$.
 391 For the other facets, a choice will have to be made.
 392 To be able to make the choice based on local information and without
 393 computing an explicit vector $\vec y$, we use the following convention.
 394 We first assign an arbitrary total order to the rays.
 395 If (the affine hull of) a facet separates the two rays not on the facet $\vec u_i$
 396 and $\vec u_j$, i.e., $\alpha_i \alpha_j > 0$ (\ref{eq:opposite}), then
 397 we choose $\vec y$ to lie on the side of the smallest ray, according
 398 to the chosen order.
 399 That is, $\sp{{\tilde n}_{ij}}y > 0$, for
 400 $\vec {\tilde n}_{ij}$ the normal of the facet pointing towards this smallest ray.
 401 Otherwise, i.e., if $\alpha_i \alpha_j < 0$,
 402 the interior of $K$ will lie on one side
 403 of the facet and then we choose $\vec y$ to lie on the other side.
 404 That is, $\sp{{\tilde n}_{ij}}y > 0$, for
 405 $\vec {\tilde n}_{ij}$ the normal of the facet pointing away from the cone $K$.
 406 Figure~\ref{fig:primal:examples} shows some example decompositions with
 407 an explicitly marked $\vec y$.
 408
 409 \begin{figure}
 410 \begin{align*}
 411 \intercol=0.48cm
 412 \begin{xy}
 413 <\intercol,0pt>:<0pt,\intercol>::
 414 \POS(-2,-1)="a"*+!U{+}
 415 \POS(2,0)="b"*+!L{+}
 416 \POS(0,2)="c"*+!D{+}
 417 \POS"a"\ar@{-}@[|(3)]"b"
 418 \POS"b"\ar@{-}@[|(3)]"c"
 419 \POS"c"\ar@{-}@[|(3)]"a"
 420 \POS(0.3,0.6)*{\bullet},*+!L{\vec y}
 421 \end{xy}
 422 & =
 423 \intercol=0.48cm
 424 \begin{xy}
 425 <\intercol,0pt>:<0pt,\intercol>::
 426 \POS(2,0)="b"*+!L{+}
 427 \POS(0,2)="c"*+!D{+}
 428 \POS(0,0)="w"*+!DR{\vec w}
 429 \POS"b"\ar@{-}@[|(3)]"c"
 430 \POS"b"\ar@{-}@[|(3)]"w"
 431 \POS"c"\ar@{-}@[|(3)]"w"
 432 \POS(0.3,0.6)*{\bullet},*+!L{\vec y}
 433 \end{xy}
 434 +
 435 \begin{xy}
 436 <\intercol,0pt>:<0pt,\intercol>::
 437 \POS(-2,-1)="a"*+!U{+}
 438 \POS(0,2)="c"*+!D{+}
 439 \POS(0,0)="w"*+!DR{\vec w}
 440 \POS"c"\ar@{-}@[|(3)]"a"
 441 \POS"a"\ar@{-}@[|(3)]"w"
 442 \POS"c"\ar@{--}@[|(3)]"w"
 443 \POS(0.3,0.6)*{\bullet},*+!L{\vec y}
 444 \end{xy}
 445 +
 446 \begin{xy}
 447 <\intercol,0pt>:<0pt,\intercol>::
 448 \POS(-2,-1)="a"*+!U{+}
 449 \POS(2,0)="b"*+!L{+}
 450 \POS(0,0)="w"*+!DR{\vec w}
 451 \POS"a"\ar@{-}@[|(3)]"b"
 452 \POS"a"\ar@{--}@[|(3)]"w"
 453 \POS"b"\ar@{--}@[|(3)]"w"
 454 \POS(0.3,0.6)*{\bullet},*+!L{\vec y}
 455 \end{xy}
 456 \\
 457 \intercol=0.48cm
 458 \begin{xy}
 459 <\intercol,0pt>:<0pt,\intercol>::
 460 \POS(-2,-1)="a"*+!U{+}
 461 \POS(2,0)="b"*+!L{+}
 462 \POS(0,2)="c"*+!D{+}
 463 \POS"a"\ar@{--}@[|(3)]"b"
 464 \POS"b"\ar@{-}@[|(3)]"c"
 465 \POS"c"\ar@{--}@[|(3)]"a"
 466 \POS(-2.5,-1.5)*{\bullet},*+!U{\vec y}
 467 \end{xy}
 468 & =
 469 \intercol=0.48cm
 470 \begin{xy}
 471 <\intercol,0pt>:<0pt,\intercol>::
 472 \POS(2,0)="b"*+!L{+}
 473 \POS(0,2)="c"*+!D{+}
 474 \POS(0,0)="w"*+!DR{\vec w}
 475 \POS"b"\ar@{-}@[|(3)]"c"
 476 \POS"b"\ar@{--}@[|(3)]"w"
 477 \POS"c"\ar@{--}@[|(3)]"w"
 478 \POS(-2.5,-1.5)*{\bullet},*+!U{\vec y}
 479 \end{xy}
 480 +
 481 \begin{xy}
 482 <\intercol,0pt>:<0pt,\intercol>::
 483 \POS(-2,-1)="a"*+!U{+}
 484 \POS(0,2)="c"*+!D{+}
 485 \POS(0,0)="w"*+!DR{\vec w}
 486 \POS"c"\ar@{--}@[|(3)]"a"
 487 \POS"a"\ar@{--}@[|(3)]"w"
 488 \POS"c"\ar@{-}@[|(3)]"w"
 489 \POS(-2.5,-1.5)*{\bullet},*+!U{\vec y}
 490 \end{xy}
 491 +
 492 \begin{xy}
 493 <\intercol,0pt>:<0pt,\intercol>::
 494 \POS(-2,-1)="a"*+!U{+}
 495 \POS(2,0)="b"*+!L{+}
 496 \POS(0,0)="w"*+!DR{\vec w}
 497 \POS"a"\ar@{--}@[|(3)]"b"
 498 \POS"a"\ar@{-}@[|(3)]"w"
 499 \POS"b"\ar@{-}@[|(3)]"w"
 500 \POS(-2.5,-1.5)*{\bullet},*+!U{\vec y}
 501 \end{xy}
 502 \\
 503 \intercol=0.48cm
 504 \begin{xy}
 505 <\intercol,0pt>:<0pt,\intercol>::
 506 \POS(-2,-1)="a"*+!U{+}
 507 \POS(2,0)="b"*+!L{+}
 508 \POS(0,2)="c"*+!D{+}
 509 \POS"a"\ar@{--}@[|(3)]"b"
 510 \POS"b"\ar@{-}@[|(3)]"c"
 511 \POS"c"\ar@{-}@[|(3)]"a"
 512 \POS(1,-1.5)*{\bullet},*+!L{\vec y}
 513 \end{xy}
 514 & =
 515 \intercol=0.48cm
 516 \begin{xy}
 517 <\intercol,0pt>:<0pt,\intercol>::
 518 \POS(2,0)="b"*+!L{+}
 519 \POS(0,2)="c"*+!D{+}
 520 \POS(0,0)="w"*+!DR{\vec w}
 521 \POS"b"\ar@{-}@[|(3)]"c"
 522 \POS"b"\ar@{--}@[|(3)]"w"
 523 \POS"c"\ar@{-}@[|(3)]"w"
 524 \POS(1,-1.5)*{\bullet},*+!L{\vec y}
 525 \end{xy}
 526 +
 527 \begin{xy}
 528 <\intercol,0pt>:<0pt,\intercol>::
 529 \POS(-2,-1)="a"*+!U{+}
 530 \POS(0,2)="c"*+!D{+}
 531 \POS(0,0)="w"*+!DR{\vec w}
 532 \POS"c"\ar@{-}@[|(3)]"a"
 533 \POS"a"\ar@{--}@[|(3)]"w"
 534 \POS"c"\ar@{--}@[|(3)]"w"
 535 \POS(1,-1.5)*{\bullet},*+!L{\vec y}
 536 \end{xy}
 537 +
 538 \begin{xy}
 539 <\intercol,0pt>:<0pt,\intercol>::
 540 \POS(-2,-1)="a"*+!U{+}
 541 \POS(2,0)="b"*+!L{+}
 542 \POS(0,0)="w"*+!DR{\vec w}
 543 \POS"a"\ar@{--}@[|(3)]"b"
 544 \POS"a"\ar@{-}@[|(3)]"w"
 545 \POS"b"\ar@{-}@[|(3)]"w"
 546 \POS(1,-1.5)*{\bullet},*+!L{\vec y}
 547 \end{xy}
 548 \\
 549 \intercol=0.48cm
 550 \begin{xy}
 551 <\intercol,0pt>:<0pt,\intercol>::
 552 \POS(-2,-1)="a"*+!U{0}
 553 \POS(2,0)="b"*+!L{+}
 554 \POS(0,2)="c"*+!D{+}
 555 \POS"a"\ar@{-}@[|(3)]"b"
 556 \POS"b"\ar@{-}@[|(3)]"c"
 557 \POS"c"\ar@{-}@[|(3)]"a"
 558 \POS(1,0.2)*{\bullet},*+!R{\vec y}
 559 \end{xy}
 560 & =
 561 \intercol=0.48cm
 562 \begin{xy}
 563 <\intercol,0pt>:<0pt,\intercol>::
 564 \POS(-2,-1)="a"*+!U{0}
 565 \POS(2,0)="b"*+!L{+}
 566 \POS(1,1)="w"*+!DL{\vec w}
 567 \POS"a"\ar@{-}@[|(3)]"b"
 568 \POS"a"\ar@{-}@[|(3)]"w"
 569 \POS"b"\ar@{-}@[|(3)]"w"
 570 \POS(1,0.2)*{\bullet},*+!R{\vec y}
 571 \end{xy}
 572 +
 573 \begin{xy}
 574 <\intercol,0pt>:<0pt,\intercol>::
 575 \POS(-2,-1)="a"*+!U{0}
 576 \POS(0,2)="c"*+!D{+}
 577 \POS(1,1)="w"*+!DL{\vec w}
 578 \POS"c"\ar@{-}@[|(3)]"a"
 579 \POS"a"\ar@{--}@[|(3)]"w"
 580 \POS"c"\ar@{-}@[|(3)]"w"
 581 \POS(1,0.2)*{\bullet},*+!R{\vec y}
 582 \end{xy}
 583 \\
 584 \intercol=0.48cm
 585 \begin{xy}
 586 <\intercol,0pt>:<0pt,\intercol>::
 587 \POS(-2,-1)="a"*+!U{0}
 588 \POS(2,0)="b"*+!U{+}
 589 \POS(0,2)="c"*+!D{-}
 590 \POS"a"\ar@{-}@[|(3)]"b"
 591 \POS"b"\ar@{--}@[|(3)]"c"
 592 \POS"c"\ar@{-}@[|(3)]"a"
 593 \POS(1.5,1.5)*{\bullet},*+!D{\vec y}
 594 \end{xy}
 595 & =
 596 \intercol=0.48cm
 597 \begin{xy}
 598 <\intercol,0pt>:<0pt,\intercol>::
 599 \POS(-2,-1)="a"*+!U{0}
 600 \POS(0,2)="c"*+!D{-}
 601 \POS(4,-2)="w"*+!L{\vec w}
 602 \POS"c"\ar@{-}@[|(3)]"a"
 603 \POS"a"\ar@{-}@[|(3)]"w"
 604 \POS"c"\ar@{--}@[|(3)]"w"
 605 \POS(1.5,1.5)*{\bullet},*+!D{\vec y}
 606 \end{xy}
 607 -
 608 \begin{xy}
 609 <\intercol,0pt>:<0pt,\intercol>::
 610 \POS(-2,-1)="a"*+!U{0}
 611 \POS(2,0)="b"*+!U{+}
 612 \POS(4,-2)="w"*+!L{\vec w}
 613 \POS"a"\ar@{--}@[|(3)]"b"
 614 \POS"a"\ar@{-}@[|(3)]"w"
 615 \POS"b"\ar@{--}@[|(3)]"w"
 616 \POS(1.5,1.5)*{\bullet},*+!D{\vec y}
 617 \end{xy}
 618 \\
 619 \intercol=0.48cm
 620 \begin{xy}
 621 <\intercol,0pt>:<0pt,\intercol>::
 622 \POS(-2,-1)="a"*+!U{0}
 623 \POS(2,0)="b"*+!U{+}
 624 \POS(0,2)="c"*+!D{-}
 625 \POS"a"\ar@{--}@[|(3)]"b"
 626 \POS"b"\ar@{--}@[|(3)]"c"
 627 \POS"c"\ar@{-}@[|(3)]"a"
 628 \POS(4.7,-2.5)*{\bullet},*+!R{\vec y}
 629 \end{xy}
 630 & =
 631 \intercol=0.48cm
 632 \begin{xy}
 633 <\intercol,0pt>:<0pt,\intercol>::
 634 \POS(-2,-1)="a"*+!U{0}
 635 \POS(0,2)="c"*+!D{-}
 636 \POS(4,-2)="w"*+!L{\vec w}
 637 \POS"c"\ar@{-}@[|(3)]"a"
 638 \POS"a"\ar@{--}@[|(3)]"w"
 639 \POS"c"\ar@{--}@[|(3)]"w"
 640 \POS(4.7,-2.5)*{\bullet},*+!R{\vec y}
 641 \end{xy}
 642 -
 643 \begin{xy}
 644 <\intercol,0pt>:<0pt,\intercol>::
 645 \POS(-2,-1)="a"*+!U{0}
 646 \POS(2,0)="b"*+!U{+}
 647 \POS(4,-2)="w"*+!L{\vec w}
 648 \POS"a"\ar@{-}@[|(3)]"b"
 649 \POS"a"\ar@{--}@[|(3)]"w"
 650 \POS"b"\ar@{--}@[|(3)]"w"
 651 \POS(4.7,-2.5)*{\bullet},*+!R{\vec y}
 652 \end{xy}
 653 \end{align*}
 654 \caption{Examples of decompositions in primal space.}
 655 \label{fig:primal:examples}
 656 \end{figure}
 657
 658 To see that there is a $\vec y$ satisfying the above constraints,
 659 we need to show that $R \cap S$ is non-empty, with
 660 $S = \{ \vec y \mid \sp{{\tilde n}_{i_kj_k}}y > 0 \text{ for all $k$}\}$.
 661 It will be easier to show this set is non-empty when the $\vec u_i$ form
 662 an ortogonal basis.  Applying a non-singular linear transformation $T$
 663 does not change the decomposition of $\vec w$ in terms of the $\vec u_i$ (i.e., the
 664 $\alpha_i$ remain unchanged), nor does this change
 665 any of the scalar products in the constraints that define $R \cap S$
 666 (the normals are transformed by $\left(T^{-1}\right)^\T$).
 667 Finding a vector $\vec y \in T(R \cap S)$ ensures that
 668 $T^{-1}(\vec y) \in R \cap S$.
 669 Without loss of generality, we can therefore assume for the purpose of
 670 showing that $R \cap S$ is non-empty that
 671 the $\vec u_i$ indeed form an orthogonal basis.
 672
 673 In the orthogonal basis, we have $\vec b_i^* = \vec u_i$
 674 and the corresponding inward normal $\vec N_i$ is either
 675 $\vec u_i$ or $-\vec u_i$.
 676 Furthermore, each normal of a facet of $S$ of the first type is of the
 677 form $\vec {\tilde n}_{i_kj_k} = a_k \vec u_{i_k} - b_k \vec u_{j_k}$, with
 678 $a_k, b_k > 0$ and ${i_k} < {j_k}$,
 679 while for the second type each normal is of the form
 680 $\vec {\tilde n}_{i_kj_k} = -a_k \vec u_{i_k} - b_k \vec u_{j_k}$, with
 681 $a_k, b_k > 0$.
 682 If $\vec {\tilde n}_{i_kj_k} = a_k \vec u_{i_k} - b_k \vec u_{j_k}$
 683 is the normal of a facet of $S$
 684 then either
 685 $(\vec N_{i_k}, \vec N_{j_k}) = (\vec u_{i_k}, \vec u_{j_k})$
 686 or
 687 $(\vec N_{i_k}, \vec N_{j_k}) = (-\vec u_{i_k}, -\vec u_{j_k})$.
 688 Otherwise, the facet would not cut $R$.
 689 Similarly,
 690 if $\vec {\tilde n}_{i_kj_k} = -a_k \vec u_{i_k} - b_k \vec u_{j_k}$
 691 is the normal of a facet of $S$
 692 then either
 693 $(\vec N_{i_k}, \vec N_{j_k}) = (\vec u_{i_k}, -\vec u_{j_k})$
 694 or
 695 $(\vec N_{i_k}, \vec N_{j_k}) = (-\vec u_{i_k}, \vec u_{j_k})$.
 696 Assume now that $R \cap S$ is empty, then there exist
 697 $\lambda_k, \mu_i \ge 0$ not all zero
 698 such that
 699 $\sum_k \lambda_k \vec {\tilde n}_{i_kj_k} + \sum_l \mu_i \vec N_i = \vec 0$.
 700 Assume $\lambda_k > 0$ for some facet of the first type.
 701 If $\vec N_{j_k} = -\vec u_{j_k}$, then $-b_k$ can only be cancelled
 702 by another facet $k'$ of the first type with $j_k = i_{k'}$, but then
 703 also $\vec N_{j_{k'}} = -\vec u_{j_{k'}}$.  Since the $j_k$ are strictly
 704 increasing, this sequence has to stop with a strictly positive coefficient
 705 for the largest $\vec u_{j_k}$ in this sequence.
 706 If, on the other hand, $\vec N_{i_k} = \vec u_{i_k}$, then $a_k$ can only
 707 be cancelled by the normal of a facet $k'$ of the second kind
 708 with $i_k = j_{k'}$, but then
 709 $\vec N_{i_{k'}} = -\vec u_{i_{k'}}$ and we return to the first case.
 710 Finally, if $\lambda_k > 0$ only for normals of facets of the second type,
 711 then either $\vec N_{i_k} = -\vec u_{i_k}$ or $\vec N_{j_k} = -\vec u_{j_k}$
 712 and so the coefficient of one of these basis vectors will be strictly
 713 negative.
 714 That is, the sum of the normals will never be zero and
 715 the set $R \cap S$ is non-empty.
 716
 717 For each ray $\vec u_j$ of cone $K_i$, i.e., the cone with $\vec u_i$ replaced
 718 by $\vec w$, we now need to determine whether the facet not containing this
 719 ray is closed or not.  We denote the (inward) normal of this cone by
 720 $\vec n_{ij}$.  Note that cone $K_j$ (if it appears in (\ref{eq:triangulations}),
 721 i.e., $\alpha_j \ne 0$) has the same facet opposite $\vec u_i$
 722 and its normal $\vec n_{ji}$ will be equal to either $\vec n_{ij}$ or
 723 $-\vec n_{ij}$, depending on whether we are dealing with an ``external'' facet,
 724 i.e., a facet of $K'$, or an ``internal'' facet.
 725 If, on the other hand, $\alpha_j = 0$, then $\vec n_{ij} = \vec n_{0j}$.
 726 If $\sp{n_{ij}}y > 0$, then the facet is closed.
 727 Otherwise it is open.
 728 It follows that the two (or more) occurrences of external facets are either all open
 729 or all closed, while for internal facets, exactly one is closed.
 730
 731 First consider the facet not containing $\vec u_0 = \vec w$.
 732 If $\alpha_i > 0$, then $\vec u_i$ and $\vec w$ are on the same side of the facet
 733 and so $\vec n_{i0} = \vec n_{0i}$.  Otherwise, $\vec n_{i0} = -\vec n_{i0}$.
 734 Second, if $\alpha_j = 0$, then replacing $\vec u_i$ by $\vec w$ does not
 735 change the affine hull of the facet and so $\vec n_{ij} = \vec n_{0j}$.
 736 Now consider the case that $\alpha_i \alpha_j < 0$, i.e., $\vec u_i$
 737 and $\vec u_j$ are on the same side of the hyperplane through the other rays.
 738 If we project $\vec u_i$, $\vec u_j$ and $\vec w$ onto a plane orthogonal
 739 to the ridge through the other rays, then the possible locations of $\vec w$
 740 with respect to $\vec u_i$ and $\vec u_j$ are shown in Figure~\ref{fig:w:same}.
 741 If both $\vec n_{0i}$ and $\vec n_{0j}$ are closed then $\vec y$ lies in region~1
 742 and therefore $\vec n_{ij}$ (as well as $\vec n_{ji}$) is closed too.
 743 Similarly, if both $\vec n_{0i}$ and $\vec n_{0j}$ are open then so is
 744 $\vec n_{ij}$.  If only one of the facets is closed, then, as explained above,
 745 we choose $\vec n_{ij}$ to be open, i.e., we take $\vec y$ to lie in region~3
 746 or~5.
 747 Figure~\ref{fig:w:opposite} shows the possible configurations
 748 for the case that $\alpha_i \alpha_j > 0$.
 749 If exactly one of $\vec n_{0i}$ and $\vec n_{0j}$ is closed, then
 750 $\vec y$ lies in region~3 or region~5 and therefore $\vec n_{ij}$ is closed iff
 751 $\vec n_{0j}$ is closed.
 752 Otherwise, as explained above, we choose $\vec n_{ij}$ to be closed if $i < j$.
 753
 754 \begin{figure}
 755 \intercol=0.7cm
 756 \begin{minipage}{0cm}
 757 \begin{xy}
 758 <\intercol,0pt>:<0pt,\intercol>::
 759 \POS(-4,0)="a"\ar@[|(2)]@{-}(4,0)
 760 \POS?(0)/\intercol/="b"\POS(2,0)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0i}}(0,0.75)}
 761 \POS(2,0)*{\bullet},*+!U{\vec u_j}
 762 \POS(-2,-3)="a"\ar@[|(2)]@{-}(2,3)
 763 \POS?(0)/\intercol/="b"\POS(1,1.5)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0j}}(0,-0.75)}
 764 \POS(1,1.5)*{\bullet},*+!R{\vec u_i}
 765 \POS(-2,3)="a"\ar@{-}(2,-3)
 766 \POS?(0)/\intercol/="b"\POS(1.5,-2.25)
 767 *\xybox{"b"-"a":(0,0)\ar_{\vec n_{ji}}^{\vec n_{ij}}(0,+0.75)}
 768 \POS(1.5,-2.25)*{\bullet},*+!R{\vec u_0 = \vec w}
 769 \POS(0,0)*{\bullet}
 770 \POS(3,1.5)*+[o][F]{\scriptstyle 1}
 771 \POS(0,2.5)*+[o][F]{\scriptstyle 2}
 772 \POS(-3,1.5)*+[o][F]{\scriptstyle 3}
 773 \POS(-3,-1.5)*+[o][F]{\scriptstyle 4}
 774 \POS(0,-3)*+[o][F]{\scriptstyle 5}
 775 \POS(3,-1.5)*+[o][F]{\scriptstyle 6}
 776 \end{xy}
 777 \end{minipage}
 778 \begin{minipage}{0cm}
 779 \begin{xy}
 780 <\intercol,0pt>:<0pt,\intercol>::
 781 \POS(-4,0)="a"\ar@[|(2)]@{-}(4,0)
 782 \POS?(0)/\intercol/="b"\POS(2,0)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0i}}(0,0.75)}
 783 \POS(2,0)*{\bullet},*+!U{\vec u_j}
 784 \POS(-2,-3)="a"\ar@[|(2)]@{-}(2,3)
 785 \POS?(0)/\intercol/="b"\POS(1,1.5)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0j}}(0,-0.75)}
 786 \POS(1,1.5)*{\bullet},*+!R{\vec u_i}
 787 \POS(-2,3)="a"\ar@{-}(2,-3)
 788 \POS?(0)/\intercol/="b"\POS(-1.5,2.25)
 789 *\xybox{"b"-"a":(0,0)\ar_{\vec n_{ji}}^{\vec n_{ij}}(0,+0.75)}
 790 \POS(-1.5,2.25)*{\bullet},*+!R{\vec u_0 = \vec w}
 791 \POS(0,0)*{\bullet}
 792 \POS(3,1.5)*+[o][F]{\scriptstyle 1}
 793 \POS(0,2.5)*+[o][F]{\scriptstyle 2}
 794 \POS(-3,1.5)*+[o][F]{\scriptstyle 3}
 795 \POS(-3,-1.5)*+[o][F]{\scriptstyle 4}
 796 \POS(0,-3)*+[o][F]{\scriptstyle 5}
 797 \POS(3,-1.5)*+[o][F]{\scriptstyle 6}
 798 \end{xy}
 799 \end{minipage}
 800 \caption{Possible locations of $\vec w$ with respect to $\vec u_i$ and
 801 $\vec u_j$, projected onto a plane orthogonal to the other rays, when
 802 $\alpha_i \alpha_j < 0$.}
 803 \label{fig:w:same}
 804 \end{figure}
 805
 806 \begin{figure}
 807 \intercol=0.7cm
 808 \begin{minipage}{0cm}
 809 \begin{xy}
 810 <\intercol,0pt>:<0pt,\intercol>::
 811 \POS(-4,0)="a"\ar@[|(2)]@{-}(4,0)
 812 \POS?(0)/\intercol/="b"\POS(2,0)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0i}}(0,0.75)}
 813 \POS(2,0)*{\bullet},*+!U{\vec u_j}
 814 \POS(-2,3)="a"\ar@[|(2)]@{-}(2,-3)
 815 \POS?(0)/\intercol/="b"\POS(-1,1.5)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0j}}(0,+0.75)}
 816 \POS(-1,1.5)*{\bullet},*+!R{\vec u_i}
 817 \POS(-2,-3)="a"\ar@{-}(2,3)
 818 \POS?(0)/\intercol/="b"\POS(1.5,2.25)
 819 *\xybox{"b"-"a":(0,0)\ar^{\vec n_{ji}}(0,+0.75)
 820 \POS(0,0)\ar_{\vec n_{ij}}(0,-0.75)}
 821 \POS(1.5,2.25)*{\bullet},*+!L{\vec u_0 = \vec w}
 822 \POS(0,0)*{\bullet}
 823 \POS(3,1.5)*+[o][F]{\scriptstyle 1}
 824 \POS(0,2.5)*+[o][F]{\scriptstyle 2}
 825 \POS(-3,1.5)*+[o][F]{\scriptstyle 3}
 826 \POS(-3,-1.5)*+[o][F]{\scriptstyle 4}
 827 \POS(0,-3)*+[o][F]{\scriptstyle 5}
 828 \POS(3,-1.5)*+[o][F]{\scriptstyle 6}
 829 \end{xy}
 830 \end{minipage}
 831 \begin{minipage}{0cm}
 832 \begin{xy}
 833 <\intercol,0pt>:<0pt,\intercol>::
 834 \POS(-4,0)="a"\ar@[|(2)]@{-}(4,0)
 835 \POS?(0)/\intercol/="b"\POS(2,0)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0i}}(0,0.75)}
 836 \POS(2,0)*{\bullet},*+!U{\vec u_j}
 837 \POS(-2,3)="a"\ar@[|(2)]@{-}(2,-3)
 838 \POS?(0)/\intercol/="b"\POS(-1,1.5)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0j}}(0,+0.75)}
 839 \POS(-1,1.5)*{\bullet},*+!R{\vec u_i}
 840 \POS(-2,-3)="a"\ar@{-}(2,3)
 841 \POS?(0)/\intercol/="b"\POS(-1.5,-2.25)
 842 *\xybox{"b"-"a":(0,0)\ar^{\vec n_{ji}}(0,+0.75)
 843 \POS(0,0)\ar_{\vec n_{ij}}(0,-0.75)}
 844 \POS(-1.5,-2.25)*{\bullet},*+!L{\vec u_0 = \vec w}
 845 \POS(0,0)*{\bullet}
 846 \POS(3,1.5)*+[o][F]{\scriptstyle 1}
 847 \POS(0,2.5)*+[o][F]{\scriptstyle 2}
 848 \POS(-3,1.5)*+[o][F]{\scriptstyle 3}
 849 \POS(-3,-1.5)*+[o][F]{\scriptstyle 4}
 850 \POS(0,-3)*+[o][F]{\scriptstyle 5}
 851 \POS(3,-1.5)*+[o][F]{\scriptstyle 6}
 852 \end{xy}
 853 \end{minipage}
 854 \caption{Possible locations of $\vec w$ with respect to $\vec u_i$ and
 855 $\vec u_j$, projected onto a plane orthogonal to the other rays, when
 856 $\alpha_i \alpha_j > 0$.}
 857 \label{fig:w:opposite}
 858 \end{figure}
 859
 860 The algorithm is summarized in Algorithm~\ref{alg:closed}, where
 861 we use the convention that in cone $K_i$, $\vec u_i$ refers to
 862 $\vec u_0 = \vec w$.
 863 Note that we do not need any of the rays or normals in this code.
 864 The only information we need is the closedness of the facets in the
 865 original cone and the signs of the $\alpha_i$.
 866
 867 \begin{algorithm}
 868 \begin{tabbing}
 869 next \= next \= next \= \kill
 870 if $\alpha_j = 0$ \\
 871 \> closed[$K_i$][$\vec u_j$] := closed[$\tilde K$][$\vec u_j$] \\
 872 else if $i = j$ \\
 873 \> if $\alpha_j > 0$ \\
 874 \> \> closed[$K_i$][$\vec u_j$] := closed[$\tilde K$][$\vec u_j$] \\
 875 \> else \\
 876 \> \> closed[$K_i$][$\vec u_j$] := $\lnot$closed[$\tilde K$][$\vec u_j$] \\
 877 else if $\alpha_i \alpha_j > 0$ \\
 878 \> if closed[$\tilde K$][$\vec u_i$] = closed[$\tilde K$][$\vec u_j$] \\
 879 \> \> closed[$K_i$][$\vec u_j$] := $i < j$ \\
 880 \> else \\
 881 \> \> closed[$K_i$][$\vec u_j$] := closed[$\tilde K$][$\vec u_j$] \\
 882 else \\
 883 \> closed[$K_i$][$\vec u_j$] := closed[$\tilde K$][$\vec u_i$] and
 884 closed[$\tilde K$][$\vec u_j$]
 885 \end{tabbing}
 886 \caption{Determine whether the facet opposite $\vec u_j$ is closed in $K_i$.}
 887 \label{alg:closed}
 888 \end{algorithm}