doc: add a note on the accuracy of our approximation
[platform/upstream/isl.git] / doc / implementation.tex
1 \section{Sets and Relations}
2
3 \begin{definition}[Polyhedral Set]
4 A {\em polyhedral set}\index{polyhedral set} $S$ is a finite union of basic sets
5 $S = \bigcup_i S_i$, each of which can be represented using affine
6 constraints
7 $$
8 S_i : \Z^n \to 2^{\Z^d} : \vec s \mapsto
9 S_i(\vec s) =
10 \{\, \vec x \in \Z^d \mid \exists \vec z \in \Z^e :
11 A \vec x + B \vec s + D \vec z + \vec c \geq \vec 0 \,\}
12 ,
13 $$
14 with $A \in \Z^{m \times d}$,
15 $B \in \Z^{m \times n}$,
16 $D \in \Z^{m \times e}$
17 and $\vec c \in \Z^m$.
18 \end{definition}
19
20 \begin{definition}[Parameter Domain of a Set]
21 Let $S \in \Z^n \to 2^{\Z^d}$ be a set.
22 The {\em parameter domain} of $S$ is the set
23 $$\pdom S \coloneqq \{\, \vec s \in \Z^n \mid S(\vec s) \ne \emptyset \,\}.$$
24 \end{definition}
25
26 \begin{definition}[Polyhedral Relation]
27 A {\em polyhedral relation}\index{polyhedral relation}
28 $R$ is a finite union of basic relations
29 $R = \bigcup_i R_i$ of type
30 $\Z^n \to 2^{\Z^{d_1+d_2}}$,
31 each of which can be represented using affine
32 constraints
33 $$
34 R_i = \vec s \mapsto
35 R_i(\vec s) =
36 \{\, \vec x_1 \to \vec x_2 \in \Z^{d_1} \times \Z^{d_2}
37 \mid \exists \vec z \in \Z^e :
38 A_1 \vec x_1 + A_2 \vec x_2 + B \vec s + D \vec z + \vec c \geq \vec 0 \,\}
39 ,
40 $$
41 with $A_i \in \Z^{m \times d_i}$,
42 $B \in \Z^{m \times n}$,
43 $D \in \Z^{m \times e}$
44 and $\vec c \in \Z^m$.
45 \end{definition}
46
47 \begin{definition}[Parameter Domain of a Relation]
48 Let $R \in \Z^n \to 2^{\Z^{d+d}}$ be a relation.
49 The {\em parameter domain} of $R$ is the set
50 $$\pdom R \coloneqq \{\, \vec s \in \Z^n \mid R(\vec s) \ne \emptyset \,\}.$$
51 \end{definition}
52
53 \begin{definition}[Domain of a Relation]
54 Let $R \in \Z^n \to 2^{\Z^{d+d}}$ be a relation.
55 The {\em domain} of $R$ is the polyhedral set
56 $$\domain R \coloneqq \vec s \mapsto
57 \{\, \vec x_1 \in \Z^{d_1} \mid \exists \vec x_2 \in \Z^{d_2} :
58 (\vec x_1, \vec x_2) \in R(\vec s) \,\}
59 .
60 $$
61 \end{definition}
62
63 \begin{definition}[Range of a Relation]
64 Let $R \in \Z^n \to 2^{\Z^{d+d}}$ be a relation.
65 The {\em range} of $R$ is the polyhedral set
66 $$
67 \range R \coloneqq \vec s \mapsto
68 \{\, \vec x_2 \in \Z^{d_2} \mid \exists \vec x_1 \in \Z^{d_1} :
69 (\vec x_1, \vec x_2) \in R(\vec s) \,\}
70 .
71 $$
72 \end{definition}
73
74 \begin{definition}[Composition of Relations]
75 Let $R \in \Z^n \to 2^{\Z^{d_1+d_2}}$ and
76 $S \in \Z^n \to 2^{\Z^{d_2+d_3}}$ be two relations,
77 then the composition of
78 $R$ and $S$ is defined as
79 $$
80 S \circ R \coloneqq
81 \vec s \mapsto
82 \{\, \vec x_1 \to \vec x_3 \in \Z^{d_1} \times \Z^{d_3}
83 \mid \exists \vec x_2 \in \Z^{d_2} :
84 \vec x_1 \to \vec x_2 \in R(\vec s) \wedge
85 \vec x_2 \to \vec x_3 \in S(\vec s)
86 \,\}
87 .
88 $$
89 \end{definition}
90
91 \begin{definition}[Difference Set of a Relation]
92 Let $R \in \Z^n \to 2^{\Z^{d+d}}$ be a relation.
93 The difference set ($\Delta \, R$) of $R$ is the set
94 of differences between image elements and the corresponding
95 domain elements,
96 $$
97 \diff R \coloneqq
98 \vec s \mapsto
99 \{\, \vec \delta \in \Z^{d} \mid \exists \vec x \to \vec y \in R :
100 \vec \delta = \vec y - \vec x
101 \,\}
102 $$
103 \end{definition}
104
105 \section{Simple Hull}\label{s:simple hull}
106
107 It is sometimes useful to have a single
108 basic set or basic relation that contains a given set or relation.
109 For rational sets, the obvious choice would be to compute the
110 (rational) convex hull.  For integer sets, the obvious choice
111 would be the integer hull.
112 However, {\tt isl} currently does not support an integer hull operation
113 and even if it did, it would be fairly expensive to compute.
114 The convex hull operation is supported, but it is also fairly
115 expensive to compute given only an implicit representation.
116
117 Usually, it is not required to compute the exact integer hull,
118 and an overapproximation of this hull is sufficient.
119 The ``simple hull'' of a set is such an overapproximation
120 and it is defined as the (inclusion-wise) smallest basic set
121 that is described by constraints that are translates of
122 the constraints in the input set.
123 This means that the simple hull is relatively cheap to compute
124 and that the number of constraints in the simple hull is no
125 larger than the number of constraints in the input.
126 \begin{definition}[Simple Hull of a Set]
127 The {\em simple hull} of a set
128 $S = \bigcup_{1 \le i \le v} S_i$, with
129 $$
130 S : \Z^n \to 2^{\Z^d} : \vec s \mapsto
131 S(\vec s) =
132 \left\{\, \vec x \in \Z^d \mid \exists \vec z \in \Z^e :
133 \bigvee_{1 \le i \le v}
134 A_i \vec x + B_i \vec s + D_i \vec z + \vec c_i \geq \vec 0 \,\right\}
135 $$
136 is the set
137 $$
138 H : \Z^n \to 2^{\Z^d} : \vec s \mapsto
139 S(\vec s) =
140 \left\{\, \vec x \in \Z^d \mid \exists \vec z \in \Z^e :
141 \bigwedge_{1 \le i \le v}
142 A_i \vec x + B_i \vec s + D_i \vec z + \vec c_i + \vec K_i \geq \vec 0
143 \,\right\}
144 ,
145 $$
146 with $\vec K_i$ the (component-wise) smallest non-negative integer vectors
147 such that $S \subseteq H$.
148 \end{definition}
149 The $\vec K_i$ can be obtained by solving a number of
150 LP problems, one for each element of each $\vec K_i$.
151 If any LP problem is unbounded, then the corresponding constraint
152 is dropped.
153
154 \section{Parametric Integer Programming}
155
156 \subsection{Introduction}\label{s:intro}
157
158 Parametric integer programming \shortcite{Feautrier88parametric}
159 is used to solve many problems within the context of the polyhedral model.
160 Here, we are mainly interested in dependence analysis \shortcite{Fea91}
161 and in computing a unique representation for existentially quantified
162 variables.  The latter operation has been used for counting elements
163 in sets involving such variables
164 \shortcite{BouletRe98,Verdoolaege2005experiences} and lies at the core
165 of the internal representation of {\tt isl}.
166
167 Parametric integer programming was first implemented in \texttt{PipLib}.
168 An alternative method for parametric integer programming
169 was later implemented in {\tt barvinok} \cite{barvinok-0.22}.
170 This method is not based on Feautrier's algorithm, but on rational
171 generating functions \cite{Woods2003short} and was inspired by the
172 ``digging'' technique of \shortciteN{DeLoera2004Three} for solving
173 non-parametric integer programming problems.
174
175 In the following sections, we briefly recall the dual simplex
176 method combined with Gomory cuts and describe some extensions
177 and optimizations.  The main algorithm is applied to a matrix
178 data structure known as a tableau.  In case of parametric problems,
179 there are two tableaus, one for the main problem and one for
180 the constraints on the parameters, known as the context tableau.
181 The handling of the context tableau is described in \autoref{s:context}.
182
183 \subsection{The Dual Simplex Method}
184
185 Tableaus can be represented in several slightly different ways.
186 In {\tt isl}, the dual simplex method uses the same representation
187 as that used by its incremental LP solver based on the \emph{primal}
188 simplex method.  The implementation of this LP solver is based
189 on that of {\tt Simplify} \shortcite{Detlefs2005simplify}, which, in turn,
190 was derived from the work of \shortciteN{Nelson1980phd}.
191 In the original \shortcite{Nelson1980phd}, the tableau was implemented
192 as a sparse matrix, but neither {\tt Simplify} nor the current
193 implementation of {\tt isl} does so.
194
195 Given some affine constraints on the variables,
196 $A \vec x + \vec b \ge \vec 0$, the tableau represents the relationship
197 between the variables $\vec x$ and non-negative variables
198 $\vec y = A \vec x + \vec b$ corresponding to the constraints.
199 The initial tableau contains $\begin{pmatrix}
200 \vec b & A
201 \end{pmatrix}$ and expresses the constraints $\vec y$ in the rows in terms
202 of the variables $\vec x$ in the columns.  The main operation defined
203 on a tableau exchanges a column and a row variable and is called a pivot.
204 During this process, some coefficients may become rational.
205 As in the \texttt{PipLib} implementation,
206 {\tt isl} maintains a shared denominator per row.
207 The sample value of a tableau is one where each column variable is assigned
208 zero and each row variable is assigned the constant term of the row.
209 This sample value represents a valid solution if each constraint variable
210 is assigned a non-negative value, i.e., if the constant terms of
211 rows corresponding to constraints are all non-negative.
212
213 The dual simplex method starts from an initial sample value that
214 may be invalid, but that is known to be (lexicographically) no
215 greater than any solution, and gradually increments this sample value
216 through pivoting until a valid solution is obtained.
217 In particular, each pivot exchanges a row variable
218 $r = -n + \sum_i a_i \, c_i$ with negative
219 sample value $-n$ with a column variable $c_j$
220 such that $a_j > 0$.  Since $c_j = (n + r - \sum_{i\ne j} a_i \, c_i)/a_j$,
221 the new row variable will have a positive sample value $n$.
222 If no such column can be found, then the problem is infeasible.
223 By always choosing the column that leads to the (lexicographically)
224 smallest increment in the variables $\vec x$,
225 the first solution found is guaranteed to be the (lexicographically)
226 minimal solution \cite{Feautrier88parametric}.
227 In order to be able to determine the smallest increment, the tableau
228 is (implicitly) extended with extra rows defining the original
229 variables in terms of the column variables.
230 If we assume that all variables are non-negative, then we know
231 that the zero vector is no greater than the minimal solution and
232 then the initial extended tableau looks as follows.
233 $$
234 \begin{tikzpicture}
235 \matrix (m) [matrix of math nodes]
236 {
237 & {} & 1 & \vec c \\
238 \vec x && |(top)| \vec 0 & I \\
239 \vec r && \vec b & |(bottom)|A \\
240 };
241 \begin{pgfonlayer}{background}
242 \node (core) [inner sep=0pt,fill=black!20,right delimiter=),left delimiter=(,fit=(top)(bottom)] {};
243 \end{pgfonlayer}
244 \end{tikzpicture}
245 $$
246 Each column in this extended tableau is lexicographically positive
247 and will remain so because of the column choice explained above.
248 It is then clear that the value of $\vec x$ will increase in each step.
249 Note that there is no need to store the extra rows explicitly.
250 If a given $x_i$ is a column variable, then the corresponding row
251 is the unit vector $e_i$.  If, on the other hand, it is a row variable,
252 then the row already appears somewhere else in the tableau.
253
254 In case of parametric problems, the sign of the constant term
255 may depend on the parameters.  Each time the constant term of a constraint row
256 changes, we therefore need to check whether the new term can attain
257 negative and/or positive values over the current set of possible
258 parameter values, i.e., the context.
259 If all these terms can only attain non-negative values, the current
260 state of the tableau represents a solution.  If one of the terms
261 can only attain non-positive values and is not identically zero,
262 the corresponding row can be pivoted.
263 Otherwise, we pick one of the terms that can attain both positive
264 and negative values and split the context into a part where
265 it only attains non-negative values and a part where it only attains
266 negative values.
267
268 \subsection{Gomory Cuts}
269
270 The solution found by the dual simplex method may have
271 non-integral coordinates.  If so, some rational solutions
272 (including the current sample value), can be cut off by
273 applying a (parametric) Gomory cut.
274 Let $r = b(\vec p) + \sp {\vec a} {\vec c}$ be the row
275 corresponding to the first non-integral coordinate of $\vec x$,
276 with $b(\vec p)$ the constant term, an affine expression in the
277 parameters $\vec p$, i.e., $b(\vec p) = \sp {\vec f} {\vec p} + g$.
278 Note that only row variables can attain
279 non-integral values as the sample value of the column variables is zero.
280 Consider the expression
281 $b(\vec p) - \ceil{b(\vec p)} + \sp {\fract{\vec a}} {\vec c}$,
282 with $\ceil\cdot$ the ceiling function and $\fract\cdot$ the
283 fractional part.  This expression is negative at the sample value
284 since $\vec c = \vec 0$ and $r = b(\vec p)$ is fractional, i.e.,
285 $\ceil{b(\vec p)} > b(\vec p)$.  On the other hand, for each integral
286 value of $r$ and $\vec c \ge 0$, the expression is non-negative
287 because $b(\vec p) - \ceil{b(\vec p)} > -1$.
288 Imposing this expression to be non-negative therefore does not
289 invalidate any integral solutions, while it does cut away the current
290 fractional sample value.  To be able to formulate this constraint,
291 a new variable $q = \floor{-b(\vec p)} = - \ceil{b(\vec p)}$ is added
292 to the context.  This integral variable is uniquely defined by the constraints
293 $0 \le -d \, b(\vec p) - d \, q \le d - 1$, with $d$ the common
294 denominator of $\vec f$ and $g$.  In practice, the variable
295 $q' = \floor{\sp {\fract{-f}} {\vec p} + \fract{-g}}$ is used instead
296 and the coefficients of the new constraint are adjusted accordingly.
297 The sign of the constant term of this new constraint need not be determined
298 as it is non-positive by construction.
299 When several of these extra context variables are added, it is important
300 to avoid adding duplicates.
301 Recent versions of {\tt PipLib} also check for such duplicates.
302
303 \subsection{Negative Unknowns and Maximization}
304
305 There are two places in the above algorithm where the unknowns $\vec x$
306 are assumed to be non-negative: the initial tableau starts from
307 sample value $\vec x = \vec 0$ and $\vec c$ is assumed to be non-negative
308 during the construction of Gomory cuts.
309 To deal with negative unknowns, \shortciteN[Appendix A.2]{Fea91}
310 proposed to use a ``big parameter'', say $M$, that is taken to be
311 an arbitrarily large positive number.  Instead of looking for the
312 lexicographically minimal value of $\vec x$, we search instead
313 for the lexicographically minimal value of $\vec x' = \vec M + \vec x$.
314 The sample value $\vec x' = \vec 0$ of the initial tableau then
315 corresponds to $\vec x = -\vec M$, which is clearly not greater than
316 any potential solution.  The sign of the constant term of a row
317 is determined lexicographically, with the coefficient of $M$ considered
318 first.  That is, if the coefficient of $M$ is not zero, then its sign
319 is the sign of the entire term.  Otherwise, the sign is determined
320 by the remaining affine expression in the parameters.
321 If the original problem has a bounded optimum, then the final sample
322 value will be of the form $\vec M + \vec v$ and the optimal value
323 of the original problem is then $\vec v$.
324 Maximization problems can be handled in a similar way by computing
325 the minimum of $\vec M - \vec x$.
326
327 When the optimum is unbounded, the optimal value computed for
328 the original problem will involve the big parameter.
329 In the original implementation of {\tt PipLib}, the big parameter could
330 even appear in some of the extra variables $\vec q$ created during
331 the application of a Gomory cut.  The final result could then contain
332 implicit conditions on the big parameter through conditions on such
333 $\vec q$ variables.  This problem was resolved in later versions
334 of {\tt PipLib} by taking $M$ to be divisible by any positive number.
335 The big parameter can then never appear in any $\vec q$ because
336 $\fract {\alpha M } = 0$.  It should be noted, though, that an unbounded
337 problem usually (but not always)
338 indicates an incorrect formulation of the problem.
339
340 The original version of {\tt PipLib} required the user to ``manually''
341 add a big parameter, perform the reformulation and interpret the result
342 \shortcite{Feautrier02}.  Recent versions allow the user to simply
343 specify that the unknowns may be negative or that the maximum should
344 be computed and then these transformations are performed internally.
345 Although there are some application, e.g.,
346 that of \shortciteN{Feautrier92multi},
347 where it is useful to have explicit control over the big parameter,
348 negative unknowns and maximization are by far the most common applications
349 of the big parameter and we believe that the user should not be bothered
350 with such implementation issues.
351 The current version of {\tt isl} therefore does not
352 provide any interface for specifying big parameters.  Instead, the user
353 can specify whether a maximum needs to be computed and no assumptions
354 are made on the sign of the unknowns.  Instead, the sign of the unknowns
355 is checked internally and a big parameter is automatically introduced when
356 needed.  For compatibility with {\tt PipLib}, the {\tt isl\_pip} tool
357 does explicitly add non-negativity constraints on the unknowns unless
358 the \verb+Urs_unknowns+ option is specified.
359 Currently, there is also no way in {\tt isl} of expressing a big
360 parameter in the output.  Even though
361 {\tt isl} makes the same divisibility assumption on the big parameter
362 as recent versions of {\tt PipLib}, it will therefore eventually
363 produce an error if the problem turns out to be unbounded.
364
365 \subsection{Preprocessing}
366
367 In this section, we describe some transformations that are
368 or can be applied in advance to reduce the running time
369 of the actual dual simplex method with Gomory cuts.
370
371 \subsubsection{Feasibility Check and Detection of Equalities}
372
373 Experience with the original {\tt PipLib} has shown that Gomory cuts
374 do not perform very well on problems that are (non-obviously) empty,
375 i.e., problems with rational solutions, but no integer solutions.
376 In {\tt isl}, we therefore first perform a feasibility check on
377 the original problem considered as a non-parametric problem
378 over the combined space of unknowns and parameters.
379 In fact, we do not simply check the feasibility, but we also
380 check for implicit equalities among the integer points by computing
381 the integer affine hull.  The algorithm used is the same as that
382 described in \autoref{s:GBR} below.
383 Computing the affine hull is fairly expensive, but it can
384 bring huge benefits if any equalities can be found or if the problem
385 turns out to be empty.
386
387 \subsubsection{Constraint Simplification}
388
389 If the coefficients of the unknown and parameters in a constraint
390 have a common factor, then this factor should be removed, possibly
391 rounding down the constant term.  For example, the constraint
392 $2 x - 5 \ge 0$ should be simplified to $x - 3 \ge 0$.
393 {\tt isl} performs such simplifications on all sets and relations.
394 Recent versions of {\tt PipLib} also perform this simplification
395 on the input.
396
397 \subsubsection{Exploiting Equalities}\label{s:equalities}
398
399 If there are any (explicit) equalities in the input description,
400 {\tt PipLib} converts each into a pair of inequalities.
401 It is also possible to write $r$ equalities as $r+1$ inequalities
402 \shortcite{Feautrier02}, but it is even better to \emph{exploit} the
403 equalities to reduce the dimensionality of the problem.
404 Given an equality involving at least one unknown, we pivot
405 the row corresponding to the equality with the column corresponding
406 to the last unknown with non-zero coefficient.  The new column variable
407 can then be removed completely because it is identically zero,
408 thereby reducing the dimensionality of the problem by one.
409 The last unknown is chosen to ensure that the columns of the initial
410 tableau remain lexicographically positive.  In particular, if
411 the equality is of the form $b + \sum_{i \le j} a_i \, x_i = 0$ with
412 $a_j \ne 0$, then the (implicit) top rows of the initial tableau
413 are changed as follows
414 $$
415 \begin{tikzpicture}
416 \matrix [matrix of math nodes]
417 {
418  & {} & |(top)| 0 & I_1 & |(j)| &  \\
419 j && 0 & & 1 & \\
420   && 0 & & & |(bottom)|I_2 \\
421 };
422 \node[overlay,above=2mm of j,anchor=south]{j};
423 \begin{pgfonlayer}{background}
424 \node (m) [inner sep=0pt,fill=black!20,right delimiter=),left delimiter=(,fit=(top)(bottom)] {};
425 \end{pgfonlayer}
426 \begin{scope}[xshift=4cm]
427 \matrix [matrix of math nodes]
428 {
429  & {} & |(top)| 0 & I_1 &  \\
430 j && |(left)| -b/a_j & -a_i/a_j & \\
431   && 0 & & |(bottom)|I_2 \\
432 };
433 \begin{pgfonlayer}{background}
434 \node (m2) [inner sep=0pt,fill=black!20,right delimiter=),left delimiter=(,fit=(top)(bottom)(left)] {};
435 \end{pgfonlayer}
436 \end{scope}
437  \draw [shorten >=7mm,-to,thick,decorate,
438         decoration={snake,amplitude=.4mm,segment length=2mm,
439                     pre=moveto,pre length=5mm,post length=8mm}]
440    (m) -- (m2);
441 \end{tikzpicture}
442 $$
443 Currently, {\tt isl} also eliminates equalities involving only parameters
444 in a similar way, provided at least one of the coefficients is equal to one.
445 The application of parameter compression (see below)
446 would obviate the need for removing parametric equalities.
447
448 \subsubsection{Offline Symmetry Detection}\label{s:offline}
449
450 Some problems, notably those of \shortciteN{Bygde2010licentiate},
451 have a collection of constraints, say
452 $b_i(\vec p) + \sp {\vec a} {\vec x} \ge 0$,
453 that only differ in their (parametric) constant terms.
454 These constant terms will be non-negative on different parts
455 of the context and this context may have to be split for each
456 of the constraints.  In the worst case, the basic algorithm may
457 have to consider all possible orderings of the constant terms.
458 Instead, {\tt isl} introduces a new parameter, say $u$, and
459 replaces the collection of constraints by the single
460 constraint $u + \sp {\vec a} {\vec x} \ge 0$ along with
461 context constraints $u \le b_i(\vec p)$.
462 Any solution to the new system is also a solution
463 to the original system since
464 $\sp {\vec a} {\vec x} \ge -u \ge -b_i(\vec p)$.
465 Conversely, $m = \min_i b_i(\vec p)$ satisfies the constraints
466 on $u$ and therefore extends a solution to the new system.
467 It can also be plugged into a new solution.
468 See \autoref{s:post} for how this substitution is currently performed
469 in {\tt isl}.
470 The method described in this section can only detect symmetries
471 that are explicitly available in the input.
472 See \autoref{s:online} for the detection
473 and exploitation of symmetries that appear during the course of
474 the dual simplex method.
475
476 \subsubsection{Parameter Compression}\label{s:compression}
477
478 It may in some cases be apparent from the equalities in the problem
479 description that there can only be a solution for a sublattice
480 of the parameters.  In such cases ``parameter compression''
481 \shortcite{Meister2004PhD,Meister2008} can be used to replace
482 the parameters by alternative ``dense'' parameters.
483 For example, if there is a constraint $2x = n$, then the system
484 will only have solutions for even values of $n$ and $n$ can be replaced
485 by $2n'$.  Similarly, the parameters $n$ and $m$ in a system with
486 the constraint $2n = 3m$ can be replaced by a single parameter $n'$
487 with $n=3n'$ and $m=2n'$.
488 It is also possible to perform a similar compression on the unknowns,
489 but it would be more complicated as the compression would have to
490 preserve the lexicographical order.  Moreover, due to our handling
491 of equalities described above there should be
492 no need for such variable compression.
493 Although parameter compression has been implemented in {\tt isl},
494 it is currently not yet used during parametric integer programming.
495
496 \subsection{Postprocessing}\label{s:post}
497
498 The output of {\tt PipLib} is a quast (quasi-affine selection tree).
499 Each internal node in this tree corresponds to a split of the context
500 based on a parametric constant term in the main tableau with indeterminate
501 sign.  Each of these nodes may introduce extra variables in the context
502 corresponding to integer divisions.  Each leaf of the tree prescribes
503 the solution in that part of the context that satisfies all the conditions
504 on the path leading to the leaf.
505 Such a quast is a very economical way of representing the solution, but
506 it would not be suitable as the (only) internal representation of
507 sets and relations in {\tt isl}.  Instead, {\tt isl} represents
508 the constraints of a set or relation in disjunctive normal form.
509 The result of a parametric integer programming problem is then also
510 converted to this internal representation.  Unfortunately, the conversion
511 to disjunctive normal form can lead to an explosion of the size
512 of the representation.
513 In some cases, this overhead would have to be paid anyway in subsequent
514 operations, but in other cases, especially for outside users that just
515 want to solve parametric integer programming problems, we would like
516 to avoid this overhead in future.  That is, we are planning on introducing
517 quasts or a related representation as one of several possible internal
518 representations and on allowing the output of {\tt isl\_pip} to optionally
519 be printed as a quast.
520
521 Currently, {\tt isl} also does not have an internal representation
522 for expressions such as $\min_i b_i(\vec p)$ from the offline
523 symmetry detection of \autoref{s:offline}.
524 Assume that one of these expressions has $n$ bounds $b_i(\vec p)$.
525 If the expression
526 does not appear in the affine expression describing the solution,
527 but only in the constraints, and if moreover, the expression
528 only appears with a positive coefficient, i.e.,
529 $\min_i b_i(\vec p) \ge f_j(\vec p)$, then each of these constraints
530 can simply be reduplicated $n$ times, once for each of the bounds.
531 Otherwise, a conversion to disjunctive normal form
532 leads to $n$ cases, each described as $u = b_i(\vec p)$ with constraints
533 $b_i(\vec p) \le b_j(\vec p)$ for $j > i$
534 and
535 $b_i(\vec p)  < b_j(\vec p)$ for $j < i$.
536 Note that even though this conversion leads to a size increase
537 by a factor of $n$, not detecting the symmetry could lead to
538 an increase by a factor of $n!$ if all possible orderings end up being
539 considered.
540
541 \subsection{Context Tableau}\label{s:context}
542
543 The main operation that a context tableau needs to provide is a test
544 on the sign of an affine expression over the elements of the context.
545 This sign can be determined by solving two integer linear feasibility
546 problems, one with a constraint added to the context that enforces
547 the expression to be non-negative and one where the expression is
548 negative.  As already mentioned by \shortciteN{Feautrier88parametric},
549 any integer linear feasibility solver could be used, but the {\tt PipLib}
550 implementation uses a recursive call to the dual simplex with Gomory
551 cuts algorithm to determine the feasibility of a context.
552 In {\tt isl}, two ways of handling the context have been implemented,
553 one that performs the recursive call and one, used by default, that
554 uses generalized basis reduction.
555 We start with some optimizations that are shared between the two
556 implementations and then discuss additional details of each of them.
557
558 \subsubsection{Maintaining Witnesses}\label{s:witness}
559
560 A common feature of both integer linear feasibility solvers is that
561 they will not only say whether a set is empty or not, but if the set
562 is non-empty, they will also provide a \emph{witness} for this result,
563 i.e., a point that belongs to the set.  By maintaining a list of such
564 witnesses, we can avoid many feasibility tests during the determination
565 of the signs of affine expressions.  In particular, if the expression
566 evaluates to a positive number on some of these points and to a negative
567 number on some others, then no feasibility test needs to be performed.
568 If all the evaluations are non-negative, we only need to check for the
569 possibility of a negative value and similarly in case of all
570 non-positive evaluations.  Finally, in the rare case that all points
571 evaluate to zero or at the start, when no points have been collected yet,
572 one or two feasibility tests need to be performed depending on the result
573 of the first test.
574
575 When a new constraint is added to the context, the points that
576 violate the constraint are temporarily removed.  They are reconsidered
577 when we backtrack over the addition of the constraint, as they will
578 satisfy the negation of the constraint.  It is only when we backtrack
579 over the addition of the points that they are finally removed completely.
580 When an extra integer division is added to the context,
581 the new coordinates of the
582 witnesses can easily be computed by evaluating the integer division.
583 The idea of keeping track of witnesses was first used in {\tt barvinok}.
584
585 \subsubsection{Choice of Constant Term on which to Split}
586
587 Recall that if there are no rows with a non-positive constant term,
588 but there are rows with an indeterminate sign, then the context
589 needs to be split along the constant term of one of these rows.
590 If there is more than one such row, then we need to choose which row
591 to split on first.  {\tt PipLib} uses a heuristic based on the (absolute)
592 sizes of the coefficients.  In particular, it takes the largest coefficient
593 of each row and then selects the row where this largest coefficient is smaller
594 than those of the other rows.
595
596 In {\tt isl}, we take that row for which non-negativity of its constant
597 term implies non-negativity of as many of the constant terms of the other
598 rows as possible.  The intuition behind this heuristic is that on the
599 positive side, we will have fewer negative and indeterminate signs,
600 while on the negative side, we need to perform a pivot, which may
601 affect any number of rows meaning that the effect on the signs
602 is difficult to predict.  This heuristic is of course much more
603 expensive to evaluate than the heuristic used by {\tt PipLib}.
604 More extensive tests are needed to evaluate whether the heuristic is worthwhile.
605
606 \subsubsection{Dual Simplex + Gomory Cuts}
607
608 When a new constraint is added to the context, the first steps
609 of the dual simplex method applied to this new context will be the same
610 or at least very similar to those taken on the original context, i.e.,
611 before the constraint was added.  In {\tt isl}, we therefore apply
612 the dual simplex method incrementally on the context and backtrack
613 to a previous state when a constraint is removed again.
614 An initial implementation that was never made public would also
615 keep the Gomory cuts, but the current implementation backtracks
616 to before the point where Gomory cuts are added before adding
617 an extra constraint to the context.
618 Keeping the Gomory cuts has the advantage that the sample value
619 is always an integer point and that this point may also satisfy
620 the new constraint.  However, due to the technique of maintaining
621 witnesses explained above,
622 we would not perform a feasibility test in such cases and then
623 the previously added cuts may be redundant, possibly resulting
624 in an accumulation of a large number of cuts.
625
626 If the parameters may be negative, then the same big parameter trick
627 used in the main tableau is applied to the context.  This big parameter
628 is of course unrelated to the big parameter from the main tableau.
629 Note that it is not a requirement for this parameter to be ``big'',
630 but it does allow for some code reuse in {\tt isl}.
631 In {\tt PipLib}, the extra parameter is not ``big'', but this may be because
632 the big parameter of the main tableau also appears
633 in the context tableau.
634
635 Finally, it was reported by \shortciteN{Galea2009personal}, who
636 worked on a parametric integer programming implementation
637 in {\tt PPL} \shortcite{PPL},
638 that it is beneficial to add cuts for \emph{all} rational coordinates
639 in the context tableau.  Based on this report,
640 the initial {\tt isl} implementation was adapted accordingly.
641
642 \subsubsection{Generalized Basis Reduction}\label{s:GBR}
643
644 The default algorithm used in {\tt isl} for feasibility checking
645 is generalized basis reduction \shortcite{Cook1991implementation}.
646 This algorithm is also used in the {\tt barvinok} implementation.
647 The algorithm is fairly robust, but it has some overhead.
648 We therefore try to avoid calling the algorithm in easy cases.
649 In particular, we incrementally keep track of points for which
650 the entire unit hypercube positioned at that point lies in the context.
651 This set is described by translates of the constraints of the context
652 and if (rationally) non-empty, any rational point
653 in the set can be rounded up to yield an integer point in the context.
654
655 A restriction of the algorithm is that it only works on bounded sets.
656 The affine hull of the recession cone therefore needs to be projected
657 out first.  As soon as the algorithm is invoked, we then also
658 incrementally keep track of this recession cone.  The reduced basis
659 found by one call of the algorithm is also reused as initial basis
660 for the next call.
661
662 Some problems lead to the
663 introduction of many integer divisions.  Within a given context,
664 some of these integer divisions may be equal to each other, even
665 if the expressions are not identical, or they may be equal to some
666 affine combination of other variables.
667 To detect such cases, we compute the affine hull of the context
668 each time a new integer division is added.  The algorithm used
669 for computing this affine hull is that of \shortciteN{Karr1976affine},
670 while the points used in this algorithm are obtained by performing
671 integer feasibility checks on that part of the context outside
672 the current approximation of the affine hull.
673 The list of witnesses is used to construct an initial approximation
674 of the hull, while any extra points found during the construction
675 of the hull is added to this list.
676 Any equality found in this way that expresses an integer division
677 as an \emph{integer} affine combination of other variables is
678 propagated to the main tableau, where it is used to eliminate that
679 integer division.
680
681 \subsection{Experiments}
682
683 \autoref{t:comparison} compares the execution times of {\tt isl}
684 (with both types of context tableau)
685 on some more difficult instances to those of other tools,
686 run on an Intel Xeon W3520 @ 2.66GHz.
687 Easier problems such as the
688 test cases distributed with {\tt Pip\-Lib} can be solved so quickly
689 that we would only be measuring overhead such as input/output and conversions
690 and not the running time of the actual algorithm.
691 We compare the following versions:
692 {\tt piplib-1.4.0-5-g0132fd9},
693 {\tt barvinok-0.32.1-73-gc5d7751},
694 {\tt isl-0.05.1-82-g3a37260}
695 and {\tt PPL} version 0.11.2.
696
697 The first test case is the following dependence analysis problem
698 originating from the Phideo project \shortcite{Verhaegh1995PhD}
699 that was communicated to us by Bart Kienhuis:
700 \begin{lstlisting}[flexiblecolumns=true,breaklines=true]{}
701 lexmax { [j1,j2] -> [i1,i2,i3,i4,i5,i6,i7,i8,i9,i10] : 1 <= i1,j1 <= 8 and 1 <= i2,i3,i4,i5,i6,i7,i8,i9,i10 <= 2 and 1 <= j2 <= 128 and i1-1 = j1-1 and i2-1+2*i3-2+4*i4-4+8*i5-8+16*i6-16+32*i7-32+64*i8-64+128*i9-128+256*i10-256=3*j2-3+66 };
702 \end{lstlisting}
703 This problem was the main inspiration
704 for some of the optimizations in \autoref{s:GBR}.
705 The second group of test cases are projections used during counting.
706 The first nine of these come from \shortciteN{Seghir2006minimizing}.
707 The remaining two come from \shortciteN{Verdoolaege2005experiences} and
708 were used to drive the first, Gomory cuts based, implementation
709 in {\tt isl}.
710 The third and final group of test cases are borrowed from
711 \shortciteN{Bygde2010licentiate} and inspired the offline symmetry detection
712 of \autoref{s:offline}.  Without symmetry detection, the running times
713 are 11s and 5.9s.
714 All running times of {\tt barvinok} and {\tt isl} include a conversion
715 to disjunctive normal form.  Without this conversion, the final two
716 cases can be solved in 0.07s and 0.21s.
717 The {\tt PipLib} implementation has some fixed limits and will
718 sometimes report the problem to be too complex (TC), while on some other
719 problems it will run out of memory (OOM).
720 The {\tt barvinok} implementation does not support problems
721 with a non-trivial lineality space (line) nor maximization problems (max).
722 The Gomory cuts based {\tt isl} implementation was terminated after 1000
723 minutes on the first problem.  The gbr version introduces some
724 overhead on some of the easier problems, but is overall the clear winner.
725
726 \begin{table}
727 \begin{center}
728 \begin{tabular}{lrrrrr}
729     & {\tt PipLib} & {\tt barvinok} & {\tt isl} cut & {\tt isl} gbr & {\tt PPL} \\
730 \hline
731 \hline
732 % bart.pip
733 Phideo & TC    & 793m   & $>$999m &   2.7s  & 372m \\
734 \hline
735 e1 & 0.33s & 3.5s & 0.08s & 0.11s & 0.18s \\
736 e3 & 0.14s & 0.13s & 0.10s & 0.10s & 0.17s \\
737 e4 & 0.24s & 9.1s & 0.09s & 0.11s & 0.70s \\
738 e5 & 0.12s & 6.0s & 0.06s & 0.14s & 0.17s \\
739 e6 & 0.10s & 6.8s & 0.17s & 0.08s & 0.21s \\
740 e7 & 0.03s & 0.27s & 0.04s & 0.04s & 0.03s \\
741 e8 & 0.03s & 0.18s & 0.03s & 0.04s & 0.01s \\
742 e9 & OOM & 70m & 2.6s & 0.94s & 22s \\
743 vd & 0.04s & 0.10s & 0.03s & 0.03s & 0.03s \\
744 bouleti & 0.25s & line & 0.06s & 0.06s & 0.15s \\
745 difficult & OOM & 1.3s & 1.7s & 0.33s & 1.4s \\
746 \hline
747 cnt/sum & TC & max & 2.2s & 2.2s & OOM \\
748 jcomplex & TC & max & 3.7s & 3.9s & OOM \\
749 \end{tabular}
750 \caption{Comparison of Execution Times}
751 \label{t:comparison}
752 \end{center}
753 \end{table}
754
755 \subsection{Online Symmetry Detection}\label{s:online}
756
757 Manual experiments on small instances of the problems of
758 \shortciteN{Bygde2010licentiate} and an analysis of the results
759 by the approximate MPA method developed by \shortciteN{Bygde2010licentiate}
760 have revealed that these problems contain many more symmetries
761 than can be detected using the offline method of \autoref{s:offline}.
762 In this section, we present an online detection mechanism that has
763 not been implemented yet, but that has shown promising results
764 in manual applications.
765
766 Let us first consider what happens when we do not perform offline
767 symmetry detection.  At some point, one of the
768 $b_i(\vec p) + \sp {\vec a} {\vec x} \ge 0$ constraints,
769 say the $j$th constraint, appears as a column
770 variable, say $c_1$, while the other constraints are represented
771 as rows of the form $b_i(\vec p) - b_j(\vec p) + c$.
772 The context is then split according to the relative order of
773 $b_j(\vec p)$ and one of the remaining $b_i(\vec p)$.
774 The offline method avoids this split by replacing all $b_i(\vec p)$
775 by a single newly introduced parameter that represents the minimum
776 of these $b_i(\vec p)$.
777 In the online method the split is similarly avoided by the introduction
778 of a new parameter.  In particular, a new parameter is introduced
779 that represents
780 $\left| b_j(\vec p) - b_i(\vec p) \right|_+ =
781 \max(b_j(\vec p) - b_i(\vec p), 0)$.
782
783 In general, let $r = b(\vec p) + \sp {\vec a} {\vec c}$ be a row
784 of the tableau such that the sign of $b(\vec p)$ is indeterminate
785 and such that exactly one of the elements of $\vec a$ is a $1$,
786 while all remaining elements are non-positive.
787 That is, $r = b(\vec p) + c_j - f$ with $f = -\sum_{i\ne j} a_i c_i \ge 0$.
788 We introduce a new parameter $t$ with
789 context constraints $t \ge -b(\vec p)$ and $t \ge 0$ and replace
790 the column variable $c_j$ by $c' + t$.  The row $r$ is now equal
791 to $b(\vec p) + t + c' - f$.  The constant term of this row is always
792 non-negative because any negative value of $b(\vec p)$ is compensated
793 by $t \ge -b(\vec p)$ while and non-negative value remains non-negative
794 because $t \ge 0$.
795
796 We need to show that this transformation does not eliminate any valid
797 solutions and that it does not introduce any spurious solutions.
798 Given a valid solution for the original problem, we need to find
799 a non-negative value of $c'$ satisfying the constraints.
800 If $b(\vec p) \ge 0$, we can take $t = 0$ so that
801 $c' = c_j - t = c_j \ge 0$.
802 If $b(\vec p) < 0$, we can take $t = -b(\vec p)$.
803 Since $r = b(\vec p) + c_j - f \ge 0$ and $f \ge 0$, we have 
804 $c' = c_j + b(\vec p) \ge 0$.
805 Note that these choices amount to plugging in
806 $t = \left|-b(\vec p)\right|_+ = \max(-b(\vec p), 0)$.
807 Conversely, given a solution to the new problem, we need to find
808 a non-negative value of $c_j$, but this is easy since $c_j = c' + t$
809 and both of these are non-negative.
810
811 Plugging in $t = \max(-b(\vec p), 0)$ can be performed as in
812 \autoref{s:post}, but, as in the case of offline symmetry detection,
813 it may be better to provide a direct representation for such
814 expressions in the internal representation of sets and relations
815 or at least in a quast-like output format.
816
817 \section{Coalescing}\label{s:coalescing}
818
819 See \shortciteN{Verdoolaege2009isl}, for now.
820 More details will be added later.
821
822 \section{Transitive Closure}
823
824 \subsection{Introduction}
825
826 \begin{definition}[Power of a Relation]
827 Let $R \in \Z^n \to 2^{\Z^{d+d}}$ be a relation and
828 $k \in \Z_{\ge 1}$
829 a positive number, then power $k$ of relation $R$ is defined as
830 \begin{equation}
831 \label{eq:transitive:power}
832 R^k \coloneqq
833 \begin{cases}
834 R & \text{if $k = 1$}
835 \\
836 R \circ R^{k-1} & \text{if $k \ge 2$}
837 .
838 \end{cases}
839 \end{equation}
840 \end{definition}
841
842 \begin{definition}[Transitive Closure of a Relation]
843 Let $R \in \Z^n \to 2^{\Z^{d+d}}$ be a relation,
844 then the transitive closure $R^+$ of $R$ is the union
845 of all positive powers of $R$,
846 $$
847 R^+ \coloneqq \bigcup_{k \ge 1} R^k
848 .
849 $$
850 \end{definition}
851 Alternatively, the transitive closure may be defined
852 inductively as
853 \begin{equation}
854 \label{eq:transitive:inductive}
855 R^+ \coloneqq R \cup \left(R \circ R^+\right)
856 .
857 \end{equation}
858
859 Since the transitive closure of a polyhedral relation
860 may no longer be a polyhedral relation \shortcite{Kelly1996closure},
861 we can, in the general case, only compute an approximation
862 of the transitive closure.
863 Whereas \shortciteN{Kelly1996closure} compute underapproximations,
864 we, like \shortciteN{Beletska2009}, compute overapproximations.
865 That is, given a relation $R$, we will compute a relation $T$
866 such that $R^+ \subseteq T$.  Of course, we want this approximation
867 to be as close as possible to the actual transitive closure
868 $R^+$ and we want to detect the cases where the approximation is
869 exact, i.e., where $T = R^+$.
870
871 For computing an approximation of the transitive closure of $R$,
872 we follow the same general strategy as \shortciteN{Beletska2009}
873 and first compute an approximation of $R^k$ for $k \ge 1$ and then project
874 out the parameter $k$ from the resulting relation.
875
876 \begin{example}
877 As a trivial example, consider the relation
878 $R = \{\, x \to x + 1 \,\}$.  The $k$th power of this map
879 for arbitrary $k$ is
880 $$
881 R^k = k \mapsto \{\, x \to x + k \mid k \ge 1 \,\}
882 .
883 $$
884 The transitive closure is then
885 $$
886 \begin{aligned}
887 R^+ & = \{\, x \to y \mid \exists k \in \Z_{\ge 1} : y = x + k \,\}
888 \\
889 & = \{\, x \to y \mid y \ge x + 1 \,\}
890 .
891 \end{aligned}
892 $$
893 \end{example}
894
895 \subsection{Computing an Approximation of $R^k$}
896 \label{s:power}
897
898 There are some special cases where the computation of $R^k$ is very easy.
899 One such case is that where $R$ does not compose with itself,
900 i.e., $R \circ R = \emptyset$ or $\domain R \cap \range R = \emptyset$.
901 In this case, $R^k$ is only non-empty for $k=1$ where it is equal
902 to $R$ itself.
903
904 In general, it is impossible to construct a closed form
905 of $R^k$ as a polyhedral relation.
906 We will therefore need to make some approximations.
907 As a first approximations, we will consider each of the basic
908 relations in $R$ as simply adding one or more offsets to a domain element
909 to arrive at an image element and ignore the fact that some of these
910 offsets may only be applied to some of the domain elements.
911 That is, we will only consider the difference set $\Delta\,R$ of the relation.
912 In particular, we will first construct a collection $P$ of paths
913 that move through
914 a total of $k$ offsets and then intersect domain and range of this
915 collection with those of $R$.
916 That is, 
917 \begin{equation}
918 \label{eq:transitive:approx}
919 K = P \cap \left(\domain R \to \range R\right)
920 ,
921 \end{equation}
922 with
923 \begin{equation}
924 \label{eq:transitive:path}
925 P = \vec s \mapsto \{\, \vec x \to \vec y \mid
926 \exists k_i \in \Z_{\ge 0}, \vec\delta_i \in k_i \, \Delta_i(\vec s) :
927 \vec y = \vec x + \sum_i \vec\delta_i
928 \wedge
929 \sum_i k_i = k > 0
930 \,\}
931 \end{equation}
932 and with $\Delta_i$ the basic sets that compose
933 the difference set $\Delta\,R$.
934 Note that the number of basic sets $\Delta_i$ need not be
935 the same as the number of basic relations in $R$.
936 Also note that since addition is commutative, it does not
937 matter in which order we add the offsets and so we are allowed
938 to group them as we did in \eqref{eq:transitive:path}.
939
940 If all the $\Delta_i$s are singleton sets
941 $\Delta_i = \{\, \vec \delta_i \,\}$ with $\vec \delta_i \in \Z^d$,
942 then \eqref{eq:transitive:path} simplifies to
943 \begin{equation}
944 \label{eq:transitive:singleton}
945 P = \{\, \vec x \to \vec y \mid
946 \exists k_i \in \Z_{\ge 0} :
947 \vec y = \vec x + \sum_i k_i \, \vec \delta_i
948 \wedge
949 \sum_i k_i = k > 0
950 \,\}
951 \end{equation}
952 and then the approximation computed in \eqref{eq:transitive:approx}
953 is essentially the same as that of \shortciteN{Beletska2009}.
954 If some of the $\Delta_i$s are not singleton sets or if
955 some of $\vec \delta_i$s are parametric, then we need
956 to resort to further approximations.
957
958 To ease both the exposition and the implementation, we will for
959 the remainder of this section work with extended offsets
960 $\Delta_i' = \Delta_i \times \{\, 1 \,\}$.
961 That is, each offset is extended with an extra coordinate that is
962 set equal to one.  The paths constructed by summing such extended
963 offsets have the length encoded as the difference of their
964 final coordinates.  The path $P'$ can then be decomposed into
965 paths $P_i'$, one for each $\Delta_i$,
966 \begin{equation}
967 \label{eq:transitive:decompose}
968 P' = \left(
969 (P_m' \cup \identity) \circ \cdots \circ
970 (P_2' \cup \identity) \circ
971 (P_1' \cup \identity)
972 \right) \cap
973 \{\,
974 \vec x' \to \vec y' \mid y_{d+1} - x_{d+1} = k > 0
975 \,\}
976 ,
977 \end{equation}
978 with
979 $$
980 P_i' = \vec s \mapsto \{\, \vec x' \to \vec y' \mid
981 \exists k \in \Z_{\ge 1}, \vec \delta \in k \, \Delta_i'(\vec s) :
982 \vec y' = \vec x' + \vec \delta
983 \,\}
984 .
985 $$
986 Note that each $P_i'$ contains paths of length at least one.
987 We therefore need to take the union with the identity relation
988 when composing the $P_i'$s to allow for paths that do not contain
989 any offsets from one or more $\Delta_i'$.
990 The path that consists of only identity relations is removed
991 by imposing the constraint $y_{d+1} - x_{d+1} > 0$.
992 Taking the union with the identity relation means that
993 that the relations we compose in \eqref{eq:transitive:decompose}
994 each consist of two basic relations.  If there are $m$
995 disjuncts in the input relation, then a direct application
996 of the composition operation may therefore result in a relation
997 with $2^m$ disjuncts, which is prohibitively expensive.
998 It is therefore crucial to apply coalescing (\autoref{s:coalescing})
999 after each composition.
1000
1001 Let us now consider how to compute an overapproximation of $P_i'$.
1002 Those that correspond to singleton $\Delta_i$s are grouped together
1003 and handled as in \eqref{eq:transitive:singleton}.
1004 Note that this is just an optimization.  The procedure described
1005 below would produce results that are at least as accurate.
1006 For simplicity, we first assume that no constraint in $\Delta_i'$
1007 involves any existentially quantified variables.
1008 We will return to existentially quantified variables at the end
1009 of this section.
1010 Without existentially quantified variables, we can classify
1011 the constraints of $\Delta_i'$ as follows
1012 \begin{enumerate}
1013 \item non-parametric constraints
1014 \begin{equation}
1015 \label{eq:transitive:non-parametric}
1016 A_1 \vec x + \vec c_1 \geq \vec 0
1017 \end{equation}
1018 \item purely parametric constraints
1019 \begin{equation}
1020 \label{eq:transitive:parametric}
1021 B_2 \vec s + \vec c_2 \geq \vec 0
1022 \end{equation}
1023 \item negative mixed constraints
1024 \begin{equation}
1025 \label{eq:transitive:mixed}
1026 A_3 \vec x + B_3 \vec s + \vec c_3 \geq \vec 0
1027 \end{equation}
1028 such that for each row $j$ and for all $\vec s$,
1029 $$
1030 \Delta_i'(\vec s) \cap
1031 \{\, \vec \delta' \mid B_{3,j} \vec s + c_{3,j} > 0 \,\}
1032 = \emptyset
1033 $$
1034 \item positive mixed constraints
1035 $$
1036 A_4 \vec x + B_4 \vec s + \vec c_4 \geq \vec 0
1037 $$
1038 such that for each row $j$, there is at least one $\vec s$ such that
1039 $$
1040 \Delta_i'(\vec s) \cap
1041 \{\, \vec \delta' \mid B_{4,j} \vec s + c_{4,j} > 0 \,\}
1042 \ne \emptyset
1043 $$
1044 \end{enumerate}
1045 We will use the following approximation $Q_i$ for $P_i'$:
1046 \begin{equation}
1047 \label{eq:transitive:Q}
1048 \begin{aligned}
1049 Q_i = \vec s \mapsto
1050 \{\,
1051 \vec x' \to \vec y'
1052 \mid {} & \exists k \in \Z_{\ge 1}, \vec f \in \Z^d :
1053 \vec y' = \vec x' + (\vec f, k)
1054 \wedge {}
1055 \\
1056 &
1057 A_1 \vec f + k \vec c_1 \geq \vec 0
1058 \wedge
1059 B_2 \vec s + \vec c_2 \geq \vec 0
1060 \wedge
1061 A_3 \vec f + B_3 \vec s + \vec c_3 \geq \vec 0
1062 \,\}
1063 .
1064 \end{aligned}
1065 \end{equation}
1066 To prove that $Q_i$ is indeed an overapproximation of $P_i'$,
1067 we need to show that for every $\vec s \in \Z^n$, for every
1068 $k \in \Z_{\ge 1}$ and for every $\vec f \in k \, \Delta_i(\vec s)$
1069 we have that
1070 $(\vec f, k)$ satisfies the constraints in \eqref{eq:transitive:Q}.
1071 If $\Delta_i(\vec s)$ is non-empty, then $\vec s$ must satisfy
1072 the constraints in \eqref{eq:transitive:parametric}.
1073 Each element $(\vec f, k) \in k \, \Delta_i'(\vec s)$ is a sum
1074 of $k$ elements $(\vec f_j, 1)$ in $\Delta_i'(\vec s)$.
1075 Each of these elements satisfies the constraints in
1076 \eqref{eq:transitive:non-parametric}, i.e.,
1077 $$
1078 \left[
1079 \begin{matrix}
1080 A_1 & \vec c_1
1081 \end{matrix}
1082 \right]
1083 \left[
1084 \begin{matrix}
1085 \vec f_j \\ 1
1086 \end{matrix}
1087 \right]
1088 \ge \vec 0
1089 .
1090 $$
1091 The sum of these elements therefore satisfies the same set of inequalities,
1092 i.e., $A_1 \vec f + k \vec c_1 \geq \vec 0$.
1093 Finally, the constraints in \eqref{eq:transitive:mixed} are such
1094 that for any $\vec s$ in the parameter domain of $\Delta$,
1095 we have $-\vec r(\vec s) \coloneqq B_3 \vec s + \vec c_3 \le \vec 0$,
1096 i.e., $A_3 \vec f_j \ge \vec r(\vec s) \ge \vec 0$
1097 and therefore also $A_3 \vec f \ge \vec r(\vec s)$.
1098 Note that if there are no mixed constraints and if the
1099 rational relaxation of $\Delta_i(\vec s)$, i.e.,
1100 $\{\, \vec x \in \Q^d \mid A_1 \vec x + \vec c_1 \ge \vec 0\,\}$,
1101 has integer vertices, then the approximation is exact, i.e.,
1102 $Q_i = P_i'$.  In this case, the vertices of $\Delta'_i(\vec s)$
1103 generate the rational cone
1104 $\{\, \vec x' \in \Q^{d+1} \mid \left[
1105 \begin{matrix}
1106 A_1 & \vec c_1
1107 \end{matrix}
1108 \right] \vec x' \,\}$ and therefore $\Delta'_i(\vec s)$ is
1109 a Hilbert basis of this cone \shortcite[Theorem~16.4]{Schrijver1986}.
1110
1111 Note however that, as pointed out by \shortciteN{DeSmet2010personal},
1112 if there \emph{are} any mixed constraints, then the above procedure may
1113 not compute the most accurate affine approximation of
1114 $k \, \Delta_i(\vec s)$ with $k \ge 1$.
1115 In particular, we only consider the negative mixed constraints that
1116 happen to appear in the description of $\Delta_i(\vec s)$, while we
1117 should instead consider \emph{all} valid such constraints.
1118 It is also sufficient to consider those constraints because any
1119 constraint that is valid for $k \, \Delta_i(\vec s)$ is also
1120 valid for $1 \, \Delta_i(\vec s) = \Delta_i(\vec s)$.
1121 Take therefore any constraint
1122 $\spv a x + \spv b s + c \ge 0$ valid for $\Delta_i(\vec s)$.
1123 This constraint is also valid for $k \, \Delta_i(\vec s)$ iff
1124 $k \, \spv a x + \spv b s + c \ge 0$.
1125 If $\spv b s + c$ can attain any positive value, then $\spv a x$
1126 may be negative for some elements of $\Delta_i(\vec s)$.
1127 We then have $k \, \spv a x < \spv a x$ for $k > 1$ and so the constraint
1128 is not valid for $k \, \Delta_i(\vec s)$.
1129 We therefore need to impose $\spv b s + c \le 0$ for all values
1130 of $\vec s$ such that $\Delta_i(\vec s)$ is non-empty, i.e.,
1131 $\vec b$ and $c$ need to be such that $- \spv b s - c \ge 0$ is a valid
1132 constraint of $\Delta_i(\vec s)$.  That is, $(\vec b, c)$ are the opposites
1133 of the coefficients of a valid constraint of $\Delta_i(\vec s)$.
1134 The approximation of $k \, \Delta_i(\vec s)$ can therefore be obtained
1135 using three applications of Farkas' lemma.  The first obtains the coefficients
1136 of constraints valid for $\Delta_i(\vec s)$.  The second obtains
1137 the coefficients of constraints valid for the projection of $\Delta_i(\vec s)$
1138 onto the parameters.  The opposite of the second set is then computed
1139 and intersected with the first set.  The result is the set of coefficients
1140 of constraints valid for $k \, \Delta_i(\vec s)$.  A final application
1141 of Farkas' lemma is needed to obtain the approximation of
1142 $k \, \Delta_i(\vec s)$ itself.
1143
1144 \begin{example}
1145 Consider the relation
1146 $$
1147 n \to \{\, (x, y) \to (1 + x, 1 - n + y) \mid n \ge 2 \,\}
1148 .
1149 $$
1150 The difference set of this relation is
1151 $$
1152 \Delta = n \to \{\, (1, 1 - n) \mid n \ge 2 \,\}
1153 .
1154 $$
1155 Using our approach, we would only consider the mixed constraint
1156 $y - 1 + n \ge 0$, leading to the following approximation of the
1157 transitive closure:
1158 $$
1159 n \to \{\, (x, y) \to (o_0, o_1) \mid n \ge 2 \wedge o_1 \le 1 - n + y \wedge o_0 \ge 1 + x \,\}
1160 .
1161 $$
1162 If, instead, we apply Farkas's lemma to $\Delta$, i.e.,
1163 \begin{verbatim}
1164 D := [n] -> { [1, 1 - n] : n >= 2 };
1165 CD := coefficients D;
1166 CD;
1167 \end{verbatim}
1168 we obtain
1169 \begin{verbatim}
1170 { rat: coefficients[[c_cst, c_n] -> [i2, i3]] : i3 <= c_n and
1171   i3 <= c_cst + 2c_n + i2 }
1172 \end{verbatim}
1173 The pure-parametric constraints valid for $\Delta$,
1174 \begin{verbatim}
1175 P := { [a,b] -> [] }(D);
1176 CP := coefficients P;
1177 CP;
1178 \end{verbatim}
1179 are
1180 \begin{verbatim}
1181 { rat: coefficients[[c_cst, c_n] -> []] : c_n >= 0 and 2c_n >= -c_cst }
1182 \end{verbatim}
1183 Negating these coefficients and intersecting with \verb+CD+,
1184 \begin{verbatim}
1185 NCP := { rat: coefficients[[a,b] -> []]
1186               -> coefficients[[-a,-b] -> []] }(CP);
1187 CK := wrap((unwrap CD) * (dom (unwrap NCP)));
1188 CK;
1189 \end{verbatim}
1190 we obtain
1191 \begin{verbatim}
1192 { rat: [[c_cst, c_n] -> [i2, i3]] : i3 <= c_n and
1193   i3 <= c_cst + 2c_n + i2 and c_n <= 0 and 2c_n <= -c_cst }
1194 \end{verbatim}
1195 The approximation for $k\,\Delta$,
1196 \begin{verbatim}
1197 K := solutions CK;
1198 K;
1199 \end{verbatim}
1200 is then
1201 \begin{verbatim}
1202 [n] -> { rat: [i0, i1] : i1 <= -i0 and i0 >= 1 and i1 <= 2 - n - i0 }
1203 \end{verbatim}
1204 Finally, the computed approximation for $R^+$,
1205 \begin{verbatim}
1206 T := unwrap({ [dx,dy] -> [[x,y] -> [x+dx,y+dy]] }(K));
1207 R := [n] -> { [x,y] -> [x+1,y+1-n] : n >= 2 };
1208 T := T * ((dom R) -> (ran R));
1209 T;
1210 \end{verbatim}
1211 is
1212 \begin{verbatim}
1213 [n] -> { [x, y] -> [o0, o1] : o1 <= x + y - o0 and
1214          o0 >= 1 + x and o1 <= 2 - n + x + y - o0 and n >= 2 }
1215 \end{verbatim}
1216 \end{example}
1217
1218 Existentially quantified variables can be handled by
1219 classifying them into variables that are uniquely
1220 determined by the parameters, variables that are independent
1221 of the parameters and others.  The first set can be treated
1222 as parameters and the second as variables.  Constraints involving
1223 the other existentially quantified variables are removed.
1224
1225 \begin{example}
1226 Consider the relation
1227 $$
1228 R =
1229 n \to \{\, x \to y \mid \exists \, \alpha_0, \alpha_1: 7\alpha_0 = -2 + n \wedge 5\alpha_1 = -1 - x + y \wedge y \ge 6 + x \,\}
1230 .
1231 $$
1232 The difference set of this relation is
1233 $$
1234 \Delta = \Delta \, R =
1235 n \to \{\, x \mid \exists \, \alpha_0, \alpha_1: 7\alpha_0 = -2 + n \wedge 5\alpha_1 = -1 + x \wedge x \ge 6 \,\}
1236 .
1237 $$
1238 The existentially quantified variables can be defined in terms
1239 of the parameters and variables as
1240 $$
1241 \alpha_0 = \floor{\frac{-2 + n}7}
1242 \qquad
1243 \text{and}
1244 \qquad
1245 \alpha_1 = \floor{\frac{-1 + x}5}
1246 .
1247 $$
1248 $\alpha_0$ can therefore be treated as a parameter,
1249 while $\alpha_1$ can be treated as a variable.
1250 This in turn means that $7\alpha_0 = -2 + n$ can be treated as
1251 a purely parametric constraint, while the other two constraints are
1252 non-parametric.
1253 The corresponding $Q$~\eqref{eq:transitive:Q} is therefore
1254 $$
1255 \begin{aligned}
1256 n \to \{\, (x,z) \to (y,w) \mid
1257 \exists\, \alpha_0, \alpha_1, k, f : {} &
1258 k \ge 1 \wedge
1259 y = x + f \wedge
1260 w = z + k \wedge {} \\
1261 &
1262 7\alpha_0 = -2 + n \wedge
1263 5\alpha_1 = -k + x \wedge
1264 x \ge 6 k
1265 \,\}
1266 .
1267 \end{aligned}
1268 $$
1269 Projecting out the final coordinates encoding the length of the paths,
1270 results in the exact transitive closure
1271 $$
1272 R^+ =
1273 n \to \{\, x \to y \mid \exists \, \alpha_0, \alpha_1: 7\alpha_1 = -2 + n \wedge 6\alpha_0 \ge -x + y \wedge 5\alpha_0 \le -1 - x + y \,\}
1274 .
1275 $$
1276 \end{example}
1277
1278 The fact that we ignore some impure constraints clearly leads
1279 to a loss of accuracy.  In some cases, some of this loss can be recovered
1280 by not considering the parameters in a special way.
1281 That is, instead of considering the set
1282 $$
1283 \Delta = \diff R =
1284 \vec s \mapsto
1285 \{\, \vec \delta \in \Z^{d} \mid \exists \vec x \to \vec y \in R :
1286 \vec \delta = \vec y - \vec x
1287 \,\}
1288 $$
1289 we consider the set
1290 $$
1291 \Delta' = \diff R' =
1292 \{\, \vec \delta \in \Z^{n+d} \mid \exists
1293 (\vec s, \vec x) \to (\vec s, \vec y) \in R' :
1294 \vec \delta = (\vec s - \vec s, \vec y - \vec x)
1295 \,\}
1296 .
1297 $$
1298 The first $n$ coordinates of every element in $\Delta'$ are zero.
1299 Projecting out these zero coordinates from $\Delta'$ is equivalent
1300 to projecting out the parameters in $\Delta$.
1301 The result is obviously a superset of $\Delta$, but all its constraints
1302 are of type \eqref{eq:transitive:non-parametric} and they can therefore
1303 all be used in the construction of $Q_i$.
1304
1305 \begin{example}
1306 Consider the relation
1307 $$
1308 % [n] -> { [x, y] -> [1 + x, 1 - n + y] | n >= 2 }
1309 R = n \to \{\, (x, y) \to (1 + x, 1 - n + y) \mid n \ge 2 \,\}
1310 .
1311 $$
1312 We have
1313 $$
1314 \diff R = n \to \{\, (1, 1 - n) \mid n \ge 2 \,\}
1315 $$
1316 and so, by treating the parameters in a special way, we obtain
1317 the following approximation for $R^+$:
1318 $$
1319 n \to \{\, (x, y) \to (x', y') \mid n \ge 2 \wedge y' \le 1 - n + y \wedge x' \ge 1 + x \,\}
1320 .
1321 $$
1322 If we consider instead
1323 $$
1324 R' = \{\, (n, x, y) \to (n, 1 + x, 1 - n + y) \mid n \ge 2 \,\}
1325 $$
1326 then
1327 $$
1328 \diff R' = \{\, (0, 1, y) \mid y \le -1 \,\}
1329 $$
1330 and we obtain the approximation
1331 $$
1332 n \to \{\, (x, y) \to (x', y') \mid n \ge 2 \wedge x' \ge 1 + x \wedge y' \le x + y - x' \,\}
1333 .
1334 $$
1335 If we consider both $\diff R$ and $\diff R'$, then we obtain
1336 $$
1337 n \to \{\, (x, y) \to (x', y') \mid n \ge 2 \wedge y' \le 1 - n + y \wedge x' \ge 1 + x \wedge y' \le x + y - x' \,\}
1338 .
1339 $$
1340 Note, however, that this is not the most accurate affine approximation that
1341 can be obtained.  That would be
1342 $$
1343 n \to \{\, (x, y) \to (x', y') \mid y' \le 2 - n + x + y - x' \wedge n \ge 2 \wedge x' \ge 1 + x \,\}
1344 .
1345 $$
1346 \end{example}
1347
1348 \subsection{Checking Exactness}
1349
1350 The approximation $T$ for the transitive closure $R^+$ can be obtained
1351 by projecting out the parameter $k$ from the approximation $K$
1352 \eqref{eq:transitive:approx} of the power $R^k$.
1353 Since $K$ is an overapproximation of $R^k$, $T$ will also be an
1354 overapproximation of $R^+$.
1355 To check whether the results are exact, we need to consider two
1356 cases depending on whether $R$ is {\em cyclic}, where $R$ is defined
1357 to be cyclic if $R^+$ maps any element to itself, i.e.,
1358 $R^+ \cap \identity \ne \emptyset$.
1359 If $R$ is acyclic, then the inductive definition of
1360 \eqref{eq:transitive:inductive} is equivalent to its completion,
1361 i.e.,
1362 $$
1363 R^+ = R \cup \left(R \circ R^+\right)
1364 $$
1365 is a defining property.
1366 Since $T$ is known to be an overapproximation, we only need to check
1367 whether
1368 $$
1369 T \subseteq R \cup \left(R \circ T\right)
1370 .
1371 $$
1372 This is essentially Theorem~5 of \shortciteN{Kelly1996closure}.
1373 The only difference is that they only consider lexicographically
1374 forward relations, a special case of acyclic relations.
1375
1376 If, on the other hand, $R$ is cyclic, then we have to resort
1377 to checking whether the approximation $K$ of the power is exact.
1378 Note that $T$ may be exact even if $K$ is not exact, so the check
1379 is sound, but incomplete.
1380 To check exactness of the power, we simply need to check
1381 \eqref{eq:transitive:power}.  Since again $K$ is known
1382 to be an overapproximation, we only need to check whether
1383 $$
1384 \begin{aligned}
1385 K'|_{y_{d+1} - x_{d+1} = 1} & \subseteq R'
1386 \\
1387 K'|_{y_{d+1} - x_{d+1} \ge 2} & \subseteq R' \circ K'|_{y_{d+1} - x_{d+1} \ge 1}
1388 ,
1389 \end{aligned}
1390 $$
1391 where $R' = \{\, \vec x' \to \vec y' \mid \vec x \to \vec y \in R
1392 \wedge y_{d+1} - x_{d+1} = 1\,\}$, i.e., $R$ extended with path
1393 lengths equal to 1.
1394
1395 All that remains is to explain how to check the cyclicity of $R$.
1396 Note that the exactness on the power is always sound, even
1397 in the acyclic case, so we only need to be careful that we find
1398 all cyclic cases.  Now, if $R$ is cyclic, i.e.,
1399 $R^+ \cap \identity \ne \emptyset$, then, since $T$ is
1400 an overapproximation of $R^+$, also
1401 $T \cap \identity \ne \emptyset$.  This in turn means
1402 that $\Delta \, K'$ contains a point whose first $d$ coordinates
1403 are zero and whose final coordinate is positive.
1404 In the implementation we currently perform this test on $P'$ instead of $K'$.
1405 Note that if $R^+$ is acyclic and $T$ is not, then the approximation
1406 is clearly not exact and the approximation of the power $K$
1407 will not be exact either.
1408
1409 \subsection{Decomposing $R$ into strongly connected components}
1410
1411 If the input relation $R$ is a union of several basic relations
1412 that can be partially ordered
1413 then the accuracy of the approximation may be improved by computing
1414 an approximation of each strongly connected components separately.
1415 For example, if $R = R_1 \cup R_2$ and $R_1 \circ R_2 = \emptyset$,
1416 then we know that any path that passes through $R_2$ cannot later
1417 pass through $R_1$, i.e.,
1418 \begin{equation}
1419 \label{eq:transitive:components}
1420 R^+ = R_1^+ \cup R_2^+ \cup \left(R_2^+ \circ R_1^+\right)
1421 .
1422 \end{equation}
1423 We can therefore compute (approximations of) transitive closures
1424 of $R_1$ and $R_2$ separately.
1425 Note, however, that the condition $R_1 \circ R_2 = \emptyset$
1426 is actually too strong.
1427 If $R_1 \circ R_2$ is a subset of $R_2 \circ R_1$
1428 then we can reorder the segments
1429 in any path that moves through both $R_1$ and $R_2$ to
1430 first move through $R_1$ and then through $R_2$.
1431
1432 This idea can be generalized to relations that are unions
1433 of more than two basic relations by constructing the
1434 strongly connected components in the graph with as vertices
1435 the basic relations and an edge between two basic relations
1436 $R_i$ and $R_j$ if $R_i$ needs to follow $R_j$ in some paths.
1437 That is, there is an edge from $R_i$ to $R_j$ iff
1438 \begin{equation}
1439 \label{eq:transitive:edge}
1440 R_i \circ R_j
1441 \not\subseteq
1442 R_j \circ R_i
1443 .
1444 \end{equation}
1445 The components can be obtained from the graph by applying
1446 Tarjan's algorithm \shortcite{Tarjan1972}.
1447
1448 In practice, we compute the (extended) powers $K_i'$ of each component
1449 separately and then compose them as in \eqref{eq:transitive:decompose}.
1450 Note, however, that in this case the order in which we apply them is
1451 important and should correspond to a topological ordering of the
1452 strongly connected components.  Simply applying Tarjan's
1453 algorithm will produce topologically sorted strongly connected components.
1454 The graph on which Tarjan's algorithm is applied is constructed on-the-fly.
1455 That is, whenever the algorithm checks if there is an edge between
1456 two vertices, we evaluate \eqref{eq:transitive:edge}.
1457 The exactness check is performed on each component separately.
1458 If the approximation turns out to be inexact for any of the components,
1459 then the entire result is marked inexact and the exactness check
1460 is skipped on the components that still need to be handled.
1461
1462 It should be noted that \eqref{eq:transitive:components}
1463 is only valid for exact transitive closures.
1464 If overapproximations are computed in the right hand side, then the result will
1465 still be an overapproximation of the left hand side, but this result
1466 may not be transitively closed.  If we only separate components based
1467 on the condition $R_i \circ R_j = \emptyset$, then there is no problem,
1468 as this condition will still hold on the computed approximations
1469 of the transitive closures.  If, however, we have exploited
1470 \eqref{eq:transitive:edge} during the decomposition and if the
1471 result turns out not to be exact, then we check whether
1472 the result is transitively closed.  If not, we recompute
1473 the transitive closure, skipping the decomposition.
1474 Note that testing for transitive closedness on the result may
1475 be fairly expensive, so we may want to make this check
1476 configurable.
1477
1478 \begin{figure}
1479 \begin{center}
1480 \begin{tikzpicture}[x=0.5cm,y=0.5cm,>=stealth,shorten >=1pt]
1481 \foreach \x in {1,...,10}{
1482     \foreach \y in {1,...,10}{
1483         \draw[->] (\x,\y) -- (\x,\y+1);
1484     }
1485 }
1486 \foreach \x in {1,...,20}{
1487     \foreach \y in {5,...,15}{
1488         \draw[->] (\x,\y) -- (\x+1,\y);
1489     }
1490 }
1491 \end{tikzpicture}
1492 \end{center}
1493 \caption{The relation from \autoref{ex:closure4}}
1494 \label{f:closure4}
1495 \end{figure}
1496 \begin{example}
1497 \label{ex:closure4}
1498 Consider the relation in example {\tt closure4} that comes with
1499 the Omega calculator~\shortcite{Omega_calc}, $R = R_1 \cup R_2$,
1500 with
1501 $$
1502 \begin{aligned}
1503 R_1 & = \{\, (x,y) \to (x,y+1) \mid 1 \le x,y \le 10 \,\}
1504 \\
1505 R_2 & = \{\, (x,y) \to (x+1,y) \mid 1 \le x \le 20 \wedge 5 \le y \le 15 \,\}
1506 .
1507 \end{aligned}
1508 $$
1509 This relation is shown graphically in \autoref{f:closure4}.
1510 We have
1511 $$
1512 \begin{aligned}
1513 R_1 \circ R_2 &=
1514 \{\, (x,y) \to (x+1,y+1) \mid 1 \le x \le 9 \wedge 5 \le y \le 10 \,\}
1515 \\
1516 R_2 \circ R_1 &=
1517 \{\, (x,y) \to (x+1,y+1) \mid 1 \le x \le 10 \wedge 4 \le y \le 10 \,\}
1518 .
1519 \end{aligned}
1520 $$
1521 Clearly, $R_1 \circ R_2 \subseteq R_2 \circ R_1$ and so
1522 $$
1523 \left(
1524 R_1 \cup R_2
1525 \right)^+
1526 =
1527 \left(R_2^+ \circ R_1^+\right)
1528 \cup R_1^+
1529 \cup R_2^+
1530 .
1531 $$
1532 \end{example}
1533
1534 \begin{figure}
1535 \newcounter{n}
1536 \newcounter{t1}
1537 \newcounter{t2}
1538 \newcounter{t3}
1539 \newcounter{t4}
1540 \begin{center}
1541 \begin{tikzpicture}[>=stealth,shorten >=1pt]
1542 \setcounter{n}{7}
1543 \foreach \i in {1,...,\value{n}}{
1544     \foreach \j in {1,...,\value{n}}{
1545         \setcounter{t1}{2 * \j - 4 - \i + 1}
1546         \setcounter{t2}{\value{n} - 3 - \i + 1}
1547         \setcounter{t3}{2 * \i - 1 - \j + 1}
1548         \setcounter{t4}{\value{n} - \j + 1}
1549         \ifnum\value{t1}>0\ifnum\value{t2}>0
1550         \ifnum\value{t3}>0\ifnum\value{t4}>0
1551             \draw[thick,->] (\i,\j) to[out=20] (\i+3,\j);
1552         \fi\fi\fi\fi
1553         \setcounter{t1}{2 * \j - 1 - \i + 1}
1554         \setcounter{t2}{\value{n} - \i + 1}
1555         \setcounter{t3}{2 * \i - 4 - \j + 1}
1556         \setcounter{t4}{\value{n} - 3 - \j + 1}
1557         \ifnum\value{t1}>0\ifnum\value{t2}>0
1558         \ifnum\value{t3}>0\ifnum\value{t4}>0
1559             \draw[thick,->] (\i,\j) to[in=-20,out=20] (\i,\j+3);
1560         \fi\fi\fi\fi
1561         \setcounter{t1}{2 * \j - 1 - \i + 1}
1562         \setcounter{t2}{\value{n} - 1 - \i + 1}
1563         \setcounter{t3}{2 * \i - 1 - \j + 1}
1564         \setcounter{t4}{\value{n} - 1 - \j + 1}
1565         \ifnum\value{t1}>0\ifnum\value{t2}>0
1566         \ifnum\value{t3}>0\ifnum\value{t4}>0
1567             \draw[thick,->] (\i,\j) to (\i+1,\j+1);
1568         \fi\fi\fi\fi
1569     }
1570 }
1571 \end{tikzpicture}
1572 \end{center}
1573 \caption{The relation from \autoref{ex:decomposition}}
1574 \label{f:decomposition}
1575 \end{figure}
1576 \begin{example}
1577 \label{ex:decomposition}
1578 Consider the relation on the right of \shortciteN[Figure~2]{Beletska2009},
1579 reproduced in \autoref{f:decomposition}.
1580 The relation can be described as $R = R_1 \cup R_2 \cup R_3$,
1581 with
1582 $$
1583 \begin{aligned}
1584 R_1 &= n \mapsto \{\, (i,j) \to (i+3,j) \mid
1585 i \le 2 j - 4 \wedge
1586 i \le n - 3 \wedge
1587 j \le 2 i - 1 \wedge
1588 j \le n \,\}
1589 \\
1590 R_2 &= n \mapsto \{\, (i,j) \to (i,j+3) \mid
1591 i \le 2 j - 1 \wedge
1592 i \le n \wedge
1593 j \le 2 i - 4 \wedge
1594 j \le n - 3 \,\}
1595 \\
1596 R_3 &= n \mapsto \{\, (i,j) \to (i+1,j+1) \mid
1597 i \le 2 j - 1 \wedge
1598 i \le n - 1 \wedge
1599 j \le 2 i - 1 \wedge
1600 j \le n - 1\,\}
1601 .
1602 \end{aligned}
1603 $$
1604 The figure shows this relation for $n = 7$.
1605 Both
1606 $R_3 \circ R_1 \subseteq R_1 \circ R_3$
1607 and
1608 $R_3 \circ R_2 \subseteq R_2 \circ R_3$,
1609 which the reader can verify using the {\tt iscc} calculator:
1610 \begin{verbatim}
1611 R1 := [n] -> { [i,j] -> [i+3,j] : i <= 2 j - 4 and i <= n - 3 and
1612                                   j <= 2 i - 1 and j <= n };
1613 R2 := [n] -> { [i,j] -> [i,j+3] : i <= 2 j - 1 and i <= n and
1614                                   j <= 2 i - 4 and j <= n - 3 };
1615 R3 := [n] -> { [i,j] -> [i+1,j+1] : i <= 2 j - 1 and i <= n - 1 and
1616                                     j <= 2 i - 1 and j <= n - 1 };
1617 (R1 . R3) - (R3 . R1);
1618 (R2 . R3) - (R3 . R2);
1619 \end{verbatim}
1620 $R_3$ can therefore be moved forward in any path.
1621 For the other two basic relations, we have both
1622 $R_2 \circ R_1 \not\subseteq R_1 \circ R_2$
1623 and
1624 $R_1 \circ R_2 \not\subseteq R_2 \circ R_1$
1625 and so $R_1$ and $R_2$ form a strongly connected component.
1626 By computing the power of $R_3$ and $R_1 \cup R_2$ separately
1627 and composing the results, the power of $R$ can be computed exactly
1628 using \eqref{eq:transitive:singleton}.
1629 As explained by \shortciteN{Beletska2009}, applying the same formula
1630 to $R$ directly, without a decomposition, would result in
1631 an overapproximation of the power.
1632 \end{example}
1633
1634 \subsection{Partitioning the domains and ranges of $R$}
1635
1636 The algorithm of \autoref{s:power} assumes that the input relation $R$
1637 can be treated as a union of translations.
1638 This is a reasonable assumption if $R$ maps elements of a given
1639 abstract domain to the same domain.
1640 However, if $R$ is a union of relations that map between different
1641 domains, then this assumption no longer holds.
1642 In particular, when an entire dependence graph is encoded
1643 in a single relation, as is done by, e.g.,
1644 \shortciteN[Section~6.1]{Barthou2000MSE}, then it does not make
1645 sense to look at differences between iterations of different domains.
1646 Now, arguably, a modified Floyd-Warshall algorithm should
1647 be applied to the dependence graph, as advocated by
1648 \shortciteN{Kelly1996closure}, with the transitive closure operation
1649 only being applied to relations from a given domain to itself.
1650 However, it is also possible to detect disjoint domains and ranges
1651 and to apply Floyd-Warshall internally.
1652
1653 \linesnumbered
1654 \begin{algorithm}
1655 \caption{The modified Floyd-Warshall algorithm of
1656 \protect\shortciteN{Kelly1996closure}}
1657 \label{a:Floyd}
1658 \SetKwInput{Input}{Input}
1659 \SetKwInput{Output}{Output}
1660 \Input{Relations $R_{pq}$, $0 \le p, q < n$}
1661 \Output{Updated relations $R_{pq}$ such that each relation
1662 $R_{pq}$ contains all indirect paths from $p$ to $q$ in the input graph}
1663 %
1664 \BlankLine
1665 \SetVline
1666 \dontprintsemicolon
1667 %
1668 \For{$r \in [0, n-1]$}{
1669     $R_{rr} \coloneqq R_{rr}^+$ \nllabel{l:Floyd:closure}\;
1670     \For{$p \in [0, n-1]$}{
1671         \For{$q \in [0, n-1]$}{
1672             \If{$p \ne r$ or $q \ne r$}{
1673                 $R_{pq} \coloneqq R_{pq} \cup \left(R_{rq} \circ R_{pr}\right)
1674                              \cup \left(R_{rq} \circ R_{rr} \circ R_{pr}\right)$
1675              \nllabel{l:Floyd:update}
1676              }
1677         }
1678     }
1679 }
1680 \end{algorithm}
1681
1682 Let the input relation $R$ be a union of $m$ basic relations $R_i$.
1683 Let $D_{2i}$ be the domains of $R_i$ and $D_{2i+1}$ the ranges of $R_i$.
1684 The first step is to group overlapping $D_j$ until a partition is
1685 obtained.  If the resulting partition consists of a single part,
1686 then we continue with the algorithm of \autoref{s:power}.
1687 Otherwise, we apply Floyd-Warshall on the graph with as vertices
1688 the parts of the partition and as edges the $R_i$ attached to
1689 the appropriate pairs of vertices.
1690 In particular, let there be $n$ parts $P_k$ in the partition.
1691 We construct $n^2$ relations
1692 $$
1693 R_{pq} \coloneqq \bigcup_{i \text{ s.t. } \domain R_i \subseteq P_p \wedge
1694                                  \range R_i \subseteq P_q} R_i
1695 ,
1696 $$
1697 apply \autoref{a:Floyd} and return the union of all resulting
1698 $R_{pq}$ as the transitive closure of $R$.
1699 Each iteration of the $r$-loop in \autoref{a:Floyd} updates
1700 all relations $R_{pq}$ to include paths that go from $p$ to $r$,
1701 possibly stay there for a while, and then go from $r$ to $q$.
1702 Note that paths that ``stay in $r$'' include all paths that
1703 pass through earlier vertices since $R_{rr}$ itself has been updated
1704 accordingly in previous iterations of the outer loop.
1705 In principle, it would be sufficient to use the $R_{pr}$
1706 and $R_{rq}$ computed in the previous iteration of the
1707 $r$-loop in Line~\ref{l:Floyd:update}.
1708 However, from an implementation perspective, it is easier
1709 to allow either or both of these to have been updated
1710 in the same iteration of the $r$-loop.
1711 This may result in duplicate paths, but these can usually
1712 be removed by coalescing (\autoref{s:coalescing}) the result of the union
1713 in Line~\ref{l:Floyd:update}, which should be done in any case.
1714 The transitive closure in Line~\ref{l:Floyd:closure}
1715 is performed using a recursive call.  This recursive call
1716 includes the partitioning step, but the resulting partition will
1717 usually be a singleton.
1718 The result of the recursive call will either be exact or an
1719 overapproximation.  The final result of Floyd-Warshall is therefore
1720 also exact or an overapproximation.
1721
1722 \begin{figure}
1723 \begin{center}
1724 \begin{tikzpicture}[x=1cm,y=1cm,>=stealth,shorten >=3pt]
1725 \foreach \x/\y in {0/0,1/1,3/2} {
1726     \fill (\x,\y) circle (2pt);
1727 }
1728 \foreach \x/\y in {0/1,2/2,3/3} {
1729     \draw (\x,\y) circle (2pt);
1730 }
1731 \draw[->] (0,0) -- (0,1);
1732 \draw[->] (0,1) -- (1,1);
1733 \draw[->] (2,2) -- (3,2);
1734 \draw[->] (3,2) -- (3,3);
1735 \draw[->,dashed] (2,2) -- (3,3);
1736 \draw[->,dotted] (0,0) -- (1,1);
1737 \end{tikzpicture}
1738 \end{center}
1739 \caption{The relation (solid arrows) on the right of Figure~1 of
1740 \protect\shortciteN{Beletska2009} and its transitive closure}
1741 \label{f:COCOA:1}
1742 \end{figure}
1743 \begin{example}
1744 Consider the relation on the right of Figure~1 of
1745 \shortciteN{Beletska2009},
1746 reproduced in \autoref{f:COCOA:1}.
1747 This relation can be described as
1748 $$
1749 \begin{aligned}
1750 \{\, (x, y) \to (x_2, y_2) \mid {} & (3y = 2x \wedge x_2 = x \wedge 3y_2 = 3 + 2x \wedge x \ge 0 \wedge x \le 3) \vee {} \\
1751 & (x_2 = 1 + x \wedge y_2 = y \wedge x \ge 0 \wedge 3y \ge 2 + 2x \wedge x \le 2 \wedge 3y \le 3 + 2x) \,\}
1752 .
1753 \end{aligned}
1754 $$
1755 Note that the domain of the upward relation overlaps with the range
1756 of the rightward relation and vice versa, but that the domain
1757 of neither relation overlaps with its own range or the domain of
1758 the other relation.
1759 The domains and ranges can therefore be partitioned into two parts,
1760 $P_0$ and $P_1$, shown as the white and black dots in \autoref{f:COCOA:1},
1761 respectively.
1762 Initially, we have
1763 $$
1764 \begin{aligned}
1765 R_{00} & = \emptyset
1766 \\
1767 R_{01} & = 
1768 \{\, (x, y) \to (x+1, y) \mid 
1769 (x \ge 0 \wedge 3y \ge 2 + 2x \wedge x \le 2 \wedge 3y \le 3 + 2x) \,\}
1770 \\
1771 R_{10} & =
1772 \{\, (x, y) \to (x_2, y_2) \mid (3y = 2x \wedge x_2 = x \wedge 3y_2 = 3 + 2x \wedge x \ge 0 \wedge x \le 3) \,\}
1773 \\
1774 R_{11} & = \emptyset
1775 .
1776 \end{aligned}
1777 $$
1778 In the first iteration, $R_{00}$ remains the same ($\emptyset^+ = \emptyset$).
1779 $R_{01}$ and $R_{10}$ are therefore also unaffected, but
1780 $R_{11}$ is updated to include $R_{01} \circ R_{10}$, i.e.,
1781 the dashed arrow in the figure.
1782 This new $R_{11}$ is obviously transitively closed, so it is not
1783 changed in the second iteration and it does not have an effect
1784 on $R_{01}$ and $R_{10}$.  However, $R_{00}$ is updated to
1785 include $R_{10} \circ R_{01}$, i.e., the dotted arrow in the figure.
1786 The transitive closure of the original relation is then equal to
1787 $R_{00} \cup R_{01} \cup R_{10} \cup R_{11}$.
1788 \end{example}
1789
1790 \subsection{Incremental Computation}
1791 \label{s:incremental}
1792
1793 In some cases it is possible and useful to compute the transitive closure
1794 of union of basic relations incrementally.  In particular,
1795 if $R$ is a union of $m$ basic maps,
1796 $$
1797 R = \bigcup_j R_j
1798 ,
1799 $$
1800 then we can pick some $R_i$ and compute the transitive closure of $R$ as
1801 \begin{equation}
1802 \label{eq:transitive:incremental}
1803 R^+ = R_i^+ \cup
1804 \left(
1805 \bigcup_{j \ne i}
1806 R_i^* \circ R_j \circ R_i^*
1807 \right)^+
1808 .
1809 \end{equation}
1810 For this approach to be successful, it is crucial that each
1811 of the disjuncts in the argument of the second transitive
1812 closure in \eqref{eq:transitive:incremental} be representable
1813 as a single basic relation, i.e., without a union.
1814 If this condition holds, then by using \eqref{eq:transitive:incremental},
1815 the number of disjuncts in the argument of the transitive closure
1816 can be reduced by one.
1817 Now, $R_i^* = R_i^+ \cup \identity$, but in some cases it is possible
1818 to relax the constraints of $R_i^+$ to include part of the identity relation,
1819 say on domain $D$.  We will use the notation
1820 ${\cal C}(R_i,D) = R_i^+ \cup \identity_D$ to represent
1821 this relaxed version of $R^+$.
1822 \shortciteN{Kelly1996closure} use the notation $R_i^?$.
1823 ${\cal C}(R_i,D)$ can be computed by allowing $k$ to attain
1824 the value $0$ in \eqref{eq:transitive:Q} and by using
1825 $$
1826 P \cap \left(D \to D\right)
1827 $$
1828 instead of \eqref{eq:transitive:approx}.
1829 Typically, $D$ will be a strict superset of both $\domain R_i$
1830 and $\range R_i$.  We therefore need to check that domain
1831 and range of the transitive closure are part of ${\cal C}(R_i,D)$,
1832 i.e., the part that results from the paths of positive length ($k \ge 1$),
1833 are equal to the domain and range of $R_i$.
1834 If not, then the incremental approach cannot be applied for
1835 the given choice of $R_i$ and $D$.
1836
1837 In order to be able to replace $R^*$ by ${\cal C}(R_i,D)$
1838 in \eqref{eq:transitive:incremental}, $D$ should be chosen
1839 to include both $\domain R$ and $\range R$, i.e., such
1840 that $\identity_D \circ R_j \circ \identity_D = R_j$ for all $j\ne i$.
1841 \shortciteN{Kelly1996closure} say that they use
1842 $D = \domain R_i \cup \range R_i$, but presumably they mean that
1843 they use $D = \domain R \cup \range R$.
1844 Now, this expression of $D$ contains a union, so it not directly usable.
1845 \shortciteN{Kelly1996closure} do not explain how they avoid this union.
1846 Apparently, in their implementation,
1847 they are using the convex hull of $\domain R \cup \range R$
1848 or at least an approximation of this convex hull.
1849 We use the simple hull (\autoref{s:simple hull}) of $\domain R \cup \range R$.
1850
1851 It is also possible to use a domain $D$ that does {\em not\/}
1852 include $\domain R \cup \range R$, but then we have to
1853 compose with ${\cal C}(R_i,D)$ more selectively.
1854 In particular, if we have
1855 \begin{equation}
1856 \label{eq:transitive:right}
1857 \text{for each $j \ne i$ either }
1858 \domain R_j \subseteq D \text{ or } \domain R_j \cap \range R_i = \emptyset
1859 \end{equation}
1860 and, similarly,
1861 \begin{equation}
1862 \label{eq:transitive:left}
1863 \text{for each $j \ne i$ either }
1864 \range R_j \subseteq D \text{ or } \range R_j \cap \domain R_i = \emptyset
1865 \end{equation}
1866 then we can refine \eqref{eq:transitive:incremental} to
1867 $$
1868 R_i^+ \cup
1869 \left(
1870 \left(
1871 \bigcup_{\shortstack{$\scriptstyle\domain R_j \subseteq D $\\
1872                      $\scriptstyle\range R_j \subseteq D$}}
1873 {\cal C} \circ R_j \circ {\cal C}
1874 \right)
1875 \cup
1876 \left(
1877 \bigcup_{\shortstack{$\scriptstyle\domain R_j \cap \range R_i = \emptyset$\\
1878                      $\scriptstyle\range R_j \subseteq D$}}
1879 \!\!\!\!\!
1880 {\cal C} \circ R_j
1881 \right)
1882 \cup
1883 \left(
1884 \bigcup_{\shortstack{$\scriptstyle\domain R_j \subseteq D $\\
1885                      $\scriptstyle\range R_j \cap \domain R_i = \emptyset$}}
1886 \!\!\!\!\!
1887 R_j \circ {\cal C}
1888 \right)
1889 \cup
1890 \left(
1891 \bigcup_{\shortstack{$\scriptstyle\domain R_j \cap \range R_i = \emptyset$\\
1892                      $\scriptstyle\range R_j \cap \domain R_i = \emptyset$}}
1893 \!\!\!\!\!
1894 R_j
1895 \right)
1896 \right)^+
1897 .
1898 $$
1899 If only property~\eqref{eq:transitive:right} holds,
1900 we can use
1901 $$
1902 R_i^+ \cup
1903 \left(
1904 \left(
1905 R_i^+ \cup \identity
1906 \right)
1907 \circ
1908 \left(
1909 \left(
1910 \bigcup_{\shortstack{$\scriptstyle\domain R_j \subseteq D $}}
1911 R_j \circ {\cal C}
1912 \right)
1913 \cup
1914 \left(
1915 \bigcup_{\shortstack{$\scriptstyle\domain R_j \cap \range R_i = \emptyset$}}
1916 \!\!\!\!\!
1917 R_j
1918 \right)
1919 \right)^+
1920 \right)
1921 ,
1922 $$
1923 while if only property~\eqref{eq:transitive:left} holds,
1924 we can use
1925 $$
1926 R_i^+ \cup
1927 \left(
1928 \left(
1929 \left(
1930 \bigcup_{\shortstack{$\scriptstyle\range R_j \subseteq D $}}
1931 {\cal C} \circ R_j
1932 \right)
1933 \cup
1934 \left(
1935 \bigcup_{\shortstack{$\scriptstyle\range R_j \cap \domain R_i = \emptyset$}}
1936 \!\!\!\!\!
1937 R_j
1938 \right)
1939 \right)^+
1940 \circ
1941 \left(
1942 R_i^+ \cup \identity
1943 \right)
1944 \right)
1945 .
1946 $$
1947
1948 It should be noted that if we want the result of the incremental
1949 approach to be transitively closed, then we can only apply it
1950 if all of the transitive closure operations involved are exact.
1951 If, say, the second transitive closure in \eqref{eq:transitive:incremental}
1952 contains extra elements, then the result does not necessarily contain
1953 the composition of these extra elements with powers of $R_i$.
1954
1955 \subsection{An {\tt Omega}-like implementation}
1956
1957 While the main algorithm of \shortciteN{Kelly1996closure} is
1958 designed to compute and underapproximation of the transitive closure,
1959 the authors mention that they could also compute overapproximations.
1960 In this section, we describe our implementation of an algorithm
1961 that is based on their ideas.
1962 Note that the {\tt Omega} library computes underapproximations
1963 \shortcite[Section 6.4]{Omega_lib}.
1964
1965 The main tool is Equation~(2) of \shortciteN{Kelly1996closure}.
1966 The input relation $R$ is first overapproximated by a ``d-form'' relation
1967 $$
1968 \{\, \vec i \to \vec j \mid \exists \vec \alpha :
1969 \vec L \le \vec j - \vec i \le \vec U
1970 \wedge
1971 (\forall p : j_p - i_p = M_p \alpha_p)
1972 \,\}
1973 ,
1974 $$
1975 where $p$ ranges over the dimensions and $\vec L$, $\vec U$ and
1976 $\vec M$ are constant integer vectors.  The elements of $\vec U$
1977 may be $\infty$, meaning that there is no upper bound corresponding
1978 to that element, and similarly for $\vec L$.
1979 Such an overapproximation can be obtained by computing strides,
1980 lower and upper bounds on the difference set $\Delta \, R$.
1981 The transitive closure of such a ``d-form'' relation is
1982 \begin{equation}
1983 \label{eq:omega}
1984 \{\, \vec i \to \vec j \mid \exists \vec \alpha, k :
1985 k \ge 1 \wedge
1986 k \, \vec L \le \vec j - \vec i \le k \, \vec U
1987 \wedge
1988 (\forall p : j_p - i_p = M_p \alpha_p)
1989 \,\}
1990 .
1991 \end{equation}
1992 The domain and range of this transitive closure are then
1993 intersected with those of the input relation.
1994 This is a special case of the algorithm in \autoref{s:power}.
1995
1996 In their algorithm for computing lower bounds, the authors
1997 use the above algorithm as a substep on the disjuncts in the relation.
1998 At the end, they say
1999 \begin{quote}
2000 If an upper bound is required, it can be calculated in a manner
2001 similar to that of a single conjunct [sic] relation.
2002 \end{quote}
2003 Presumably, the authors mean that a ``d-form'' approximation
2004 of the whole input relation should be used.
2005 However, the accuracy can be improved by also trying to
2006 apply the incremental technique from the same paper,
2007 which is explained in more detail in \autoref{s:incremental}.
2008 In this case, ${\cal C}(R_i,D)$ can be obtained by
2009 allowing the value zero for $k$ in \eqref{eq:omega},
2010 i.e., by computing
2011 $$
2012 \{\, \vec i \to \vec j \mid \exists \vec \alpha, k :
2013 k \ge 0 \wedge
2014 k \, \vec L \le \vec j - \vec i \le k \, \vec U
2015 \wedge
2016 (\forall p : j_p - i_p = M_p \alpha_p)
2017 \,\}
2018 .
2019 $$
2020 In our implementation we take as $D$ the simple hull
2021 (\autoref{s:simple hull}) of $\domain R \cup \range R$.
2022 To determine whether it is safe to use ${\cal C}(R_i,D)$,
2023 we check the following conditions, as proposed by
2024 \shortciteN{Kelly1996closure}:
2025 ${\cal C}(R_i,D) - R_i^+$ is not a union and for each $j \ne i$
2026 the condition
2027 $$
2028 \left({\cal C}(R_i,D) - R_i^+\right)
2029 \circ
2030 R_j
2031 \circ
2032 \left({\cal C}(R_i,D) - R_i^+\right)
2033 =
2034 R_j
2035 $$
2036 holds.