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