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 \shortcite{Feautrier88parametric} 159is used to solve many problems within the context of the polyhedral model. 160Here, we are mainly interested in dependence analysis \shortcite{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\shortcite{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 \shortciteN{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} \shortcite{Detlefs2005simplify}, which, in turn, 190was derived from the work of \shortciteN{Nelson1980phd}. 191In the original \shortcite{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, \shortciteN[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\shortcite{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 \shortciteN{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\shortcite{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 \shortciteN{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 476\subsubsection{Parameter Compression}\label{s:compression} 477 478It may in some cases be apparent from the equalities in the problem 479description that there can only be a solution for a sublattice 480of the parameters. In such cases ``parameter compression'' 481\shortcite{Meister2004PhD,Meister2008} can be used to replace 482the parameters by alternative ``dense'' parameters. 483For example, if there is a constraint $2x = n$, then the system 484will only have solutions for even values of $n$ and $n$ can be replaced 485by $2n'$. Similarly, the parameters $n$ and $m$ in a system with 486the constraint $2n = 3m$ can be replaced by a single parameter $n'$ 487with $n=3n'$ and $m=2n'$. 488It is also possible to perform a similar compression on the unknowns, 489but it would be more complicated as the compression would have to 490preserve the lexicographical order. Moreover, due to our handling 491of equalities described above there should be 492no need for such variable compression. 493Although parameter compression has been implemented in {\tt isl}, 494it is currently not yet used during parametric integer programming. 495 496\subsection{Postprocessing}\label{s:post} 497 498The output of {\tt PipLib} is a quast (quasi-affine selection tree). 499Each internal node in this tree corresponds to a split of the context 500based on a parametric constant term in the main tableau with indeterminate 501sign. Each of these nodes may introduce extra variables in the context 502corresponding to integer divisions. Each leaf of the tree prescribes 503the solution in that part of the context that satisfies all the conditions 504on the path leading to the leaf. 505Such a quast is a very economical way of representing the solution, but 506it would not be suitable as the (only) internal representation of 507sets and relations in {\tt isl}. Instead, {\tt isl} represents 508the constraints of a set or relation in disjunctive normal form. 509The result of a parametric integer programming problem is then also 510converted to this internal representation. Unfortunately, the conversion 511to disjunctive normal form can lead to an explosion of the size 512of the representation. 513In some cases, this overhead would have to be paid anyway in subsequent 514operations, but in other cases, especially for outside users that just 515want to solve parametric integer programming problems, we would like 516to avoid this overhead in future. That is, we are planning on introducing 517quasts or a related representation as one of several possible internal 518representations and on allowing the output of {\tt isl\_pip} to optionally 519be printed as a quast. 520 521Currently, {\tt isl} also does not have an internal representation 522for expressions such as $\min_i b_i(\vec p)$ from the offline 523symmetry detection of \autoref{s:offline}. 524Assume that one of these expressions has $n$ bounds $b_i(\vec p)$. 525If the expression 526does not appear in the affine expression describing the solution, 527but only in the constraints, and if moreover, the expression 528only appears with a positive coefficient, i.e., 529$\min_i b_i(\vec p) \ge f_j(\vec p)$, then each of these constraints 530can simply be reduplicated $n$ times, once for each of the bounds. 531Otherwise, a conversion to disjunctive normal form 532leads 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$ 534and 535$b_i(\vec p) < b_j(\vec p)$ for $j < i$. 536Note that even though this conversion leads to a size increase 537by a factor of $n$, not detecting the symmetry could lead to 538an increase by a factor of $n!$ if all possible orderings end up being 539considered. 540 541\subsection{Context Tableau}\label{s:context} 542 543The main operation that a context tableau needs to provide is a test 544on the sign of an affine expression over the elements of the context. 545This sign can be determined by solving two integer linear feasibility 546problems, one with a constraint added to the context that enforces 547the expression to be non-negative and one where the expression is 548negative. As already mentioned by \shortciteN{Feautrier88parametric}, 549any integer linear feasibility solver could be used, but the {\tt PipLib} 550implementation uses a recursive call to the dual simplex with Gomory 551cuts algorithm to determine the feasibility of a context. 552In {\tt isl}, two ways of handling the context have been implemented, 553one that performs the recursive call and one, used by default, that 554uses generalized basis reduction. 555We start with some optimizations that are shared between the two 556implementations and then discuss additional details of each of them. 557 558\subsubsection{Maintaining Witnesses}\label{s:witness} 559 560A common feature of both integer linear feasibility solvers is that 561they will not only say whether a set is empty or not, but if the set 562is non-empty, they will also provide a \emph{witness} for this result, 563i.e., a point that belongs to the set. By maintaining a list of such 564witnesses, we can avoid many feasibility tests during the determination 565of the signs of affine expressions. In particular, if the expression 566evaluates to a positive number on some of these points and to a negative 567number on some others, then no feasibility test needs to be performed. 568If all the evaluations are non-negative, we only need to check for the 569possibility of a negative value and similarly in case of all 570non-positive evaluations. Finally, in the rare case that all points 571evaluate to zero or at the start, when no points have been collected yet, 572one or two feasibility tests need to be performed depending on the result 573of the first test. 574 575When a new constraint is added to the context, the points that 576violate the constraint are temporarily removed. They are reconsidered 577when we backtrack over the addition of the constraint, as they will 578satisfy the negation of the constraint. It is only when we backtrack 579over the addition of the points that they are finally removed completely. 580When an extra integer division is added to the context, 581the new coordinates of the 582witnesses can easily be computed by evaluating the integer division. 583The idea of keeping track of witnesses was first used in {\tt barvinok}. 584 585\subsubsection{Choice of Constant Term on which to Split} 586 587Recall that if there are no rows with a non-positive constant term, 588but there are rows with an indeterminate sign, then the context 589needs to be split along the constant term of one of these rows. 590If there is more than one such row, then we need to choose which row 591to split on first. {\tt PipLib} uses a heuristic based on the (absolute) 592sizes of the coefficients. In particular, it takes the largest coefficient 593of each row and then selects the row where this largest coefficient is smaller 594than those of the other rows. 595 596In {\tt isl}, we take that row for which non-negativity of its constant 597term implies non-negativity of as many of the constant terms of the other 598rows as possible. The intuition behind this heuristic is that on the 599positive side, we will have fewer negative and indeterminate signs, 600while on the negative side, we need to perform a pivot, which may 601affect any number of rows meaning that the effect on the signs 602is difficult to predict. This heuristic is of course much more 603expensive to evaluate than the heuristic used by {\tt PipLib}. 604More extensive tests are needed to evaluate whether the heuristic is worthwhile. 605 606\subsubsection{Dual Simplex + Gomory Cuts} 607 608When a new constraint is added to the context, the first steps 609of the dual simplex method applied to this new context will be the same 610or at least very similar to those taken on the original context, i.e., 611before the constraint was added. In {\tt isl}, we therefore apply 612the dual simplex method incrementally on the context and backtrack 613to a previous state when a constraint is removed again. 614An initial implementation that was never made public would also 615keep the Gomory cuts, but the current implementation backtracks 616to before the point where Gomory cuts are added before adding 617an extra constraint to the context. 618Keeping the Gomory cuts has the advantage that the sample value 619is always an integer point and that this point may also satisfy 620the new constraint. However, due to the technique of maintaining 621witnesses explained above, 622we would not perform a feasibility test in such cases and then 623the previously added cuts may be redundant, possibly resulting 624in an accumulation of a large number of cuts. 625 626If the parameters may be negative, then the same big parameter trick 627used in the main tableau is applied to the context. This big parameter 628is of course unrelated to the big parameter from the main tableau. 629Note that it is not a requirement for this parameter to be ``big'', 630but it does allow for some code reuse in {\tt isl}. 631In {\tt PipLib}, the extra parameter is not ``big'', but this may be because 632the big parameter of the main tableau also appears 633in the context tableau. 634 635Finally, it was reported by \shortciteN{Galea2009personal}, who 636worked on a parametric integer programming implementation 637in {\tt PPL} \shortcite{PPL}, 638that it is beneficial to add cuts for \emph{all} rational coordinates 639in the context tableau. Based on this report, 640the initial {\tt isl} implementation was adapted accordingly. 641 642\subsubsection{Generalized Basis Reduction}\label{s:GBR} 643 644The default algorithm used in {\tt isl} for feasibility checking 645is generalized basis reduction \shortcite{Cook1991implementation}. 646This algorithm is also used in the {\tt barvinok} implementation. 647The algorithm is fairly robust, but it has some overhead. 648We therefore try to avoid calling the algorithm in easy cases. 649In particular, we incrementally keep track of points for which 650the entire unit hypercube positioned at that point lies in the context. 651This set is described by translates of the constraints of the context 652and if (rationally) non-empty, any rational point 653in the set can be rounded up to yield an integer point in the context. 654 655A restriction of the algorithm is that it only works on bounded sets. 656The affine hull of the recession cone therefore needs to be projected 657out first. As soon as the algorithm is invoked, we then also 658incrementally keep track of this recession cone. The reduced basis 659found by one call of the algorithm is also reused as initial basis 660for the next call. 661 662Some problems lead to the 663introduction of many integer divisions. Within a given context, 664some of these integer divisions may be equal to each other, even 665if the expressions are not identical, or they may be equal to some 666affine combination of other variables. 667To detect such cases, we compute the affine hull of the context 668each time a new integer division is added. The algorithm used 669for computing this affine hull is that of \shortciteN{Karr1976affine}, 670while the points used in this algorithm are obtained by performing 671integer feasibility checks on that part of the context outside 672the current approximation of the affine hull. 673The list of witnesses is used to construct an initial approximation 674of the hull, while any extra points found during the construction 675of the hull is added to this list. 676Any equality found in this way that expresses an integer division 677as an \emph{integer} affine combination of other variables is 678propagated to the main tableau, where it is used to eliminate that 679integer division. 680 681\subsection{Experiments} 682 683\autoref{t:comparison} compares the execution times of {\tt isl} 684(with both types of context tableau) 685on some more difficult instances to those of other tools, 686run on an Intel Xeon W3520 @ 2.66GHz. 687Easier problems such as the 688test cases distributed with {\tt Pip\-Lib} can be solved so quickly 689that we would only be measuring overhead such as input/output and conversions 690and not the running time of the actual algorithm. 691We 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} 695and {\tt PPL} version 0.11.2. 696 697The first test case is the following dependence analysis problem 698originating from the Phideo project \shortcite{Verhaegh1995PhD} 699that was communicated to us by Bart Kienhuis: 700\begin{lstlisting}[flexiblecolumns=true,breaklines=true]{} 701lexmax { [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} 703This problem was the main inspiration 704for some of the optimizations in \autoref{s:GBR}. 705The second group of test cases are projections used during counting. 706The first nine of these come from \shortciteN{Seghir2006minimizing}. 707The remaining two come from \shortciteN{Verdoolaege2005experiences} and 708were used to drive the first, Gomory cuts based, implementation 709in {\tt isl}. 710The third and final group of test cases are borrowed from 711\shortciteN{Bygde2010licentiate} and inspired the offline symmetry detection 712of \autoref{s:offline}. Without symmetry detection, the running times 713are 11s and 5.9s. 714All running times of {\tt barvinok} and {\tt isl} include a conversion 715to disjunctive normal form. Without this conversion, the final two 716cases can be solved in 0.07s and 0.21s. 717The {\tt PipLib} implementation has some fixed limits and will 718sometimes report the problem to be too complex (TC), while on some other 719problems it will run out of memory (OOM). 720The {\tt barvinok} implementation does not support problems 721with a non-trivial lineality space (line) nor maximization problems (max). 722The Gomory cuts based {\tt isl} implementation was terminated after 1000 723minutes on the first problem. The gbr version introduces some 724overhead 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 733Phideo & TC & 793m & $>$999m & 2.7s & 372m \\ 734\hline 735e1 & 0.33s & 3.5s & 0.08s & 0.11s & 0.18s \\ 736e3 & 0.14s & 0.13s & 0.10s & 0.10s & 0.17s \\ 737e4 & 0.24s & 9.1s & 0.09s & 0.11s & 0.70s \\ 738e5 & 0.12s & 6.0s & 0.06s & 0.14s & 0.17s \\ 739e6 & 0.10s & 6.8s & 0.17s & 0.08s & 0.21s \\ 740e7 & 0.03s & 0.27s & 0.04s & 0.04s & 0.03s \\ 741e8 & 0.03s & 0.18s & 0.03s & 0.04s & 0.01s \\ 742e9 & OOM & 70m & 2.6s & 0.94s & 22s \\ 743vd & 0.04s & 0.10s & 0.03s & 0.03s & 0.03s \\ 744bouleti & 0.25s & line & 0.06s & 0.06s & 0.15s \\ 745difficult & OOM & 1.3s & 1.7s & 0.33s & 1.4s \\ 746\hline 747cnt/sum & TC & max & 2.2s & 2.2s & OOM \\ 748jcomplex & 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 757Manual experiments on small instances of the problems of 758\shortciteN{Bygde2010licentiate} and an analysis of the results 759by the approximate MPA method developed by \shortciteN{Bygde2010licentiate} 760have revealed that these problems contain many more symmetries 761than can be detected using the offline method of \autoref{s:offline}. 762In this section, we present an online detection mechanism that has 763not been implemented yet, but that has shown promising results 764in manual applications. 765 766Let us first consider what happens when we do not perform offline 767symmetry detection. At some point, one of the 768$b_i(\vec p) + \sp {\vec a} {\vec x} \ge 0$ constraints, 769say the $j$th constraint, appears as a column 770variable, say $c_1$, while the other constraints are represented 771as rows of the form $b_i(\vec p) - b_j(\vec p) + c$. 772The context is then split according to the relative order of 773$b_j(\vec p)$ and one of the remaining $b_i(\vec p)$. 774The offline method avoids this split by replacing all $b_i(\vec p)$ 775by a single newly introduced parameter that represents the minimum 776of these $b_i(\vec p)$. 777In the online method the split is similarly avoided by the introduction 778of a new parameter. In particular, a new parameter is introduced 779that represents 780$\left| b_j(\vec p) - b_i(\vec p) \right|_+ = 781\max(b_j(\vec p) - b_i(\vec p), 0)$. 782 783In general, let $r = b(\vec p) + \sp {\vec a} {\vec c}$ be a row 784of the tableau such that the sign of $b(\vec p)$ is indeterminate 785and such that exactly one of the elements of $\vec a$ is a $1$, 786while all remaining elements are non-positive. 787That is, $r = b(\vec p) + c_j - f$ with $f = -\sum_{i\ne j} a_i c_i \ge 0$. 788We introduce a new parameter $t$ with 789context constraints $t \ge -b(\vec p)$ and $t \ge 0$ and replace 790the column variable $c_j$ by $c' + t$. The row $r$ is now equal 791to $b(\vec p) + t + c' - f$. The constant term of this row is always 792non-negative because any negative value of $b(\vec p)$ is compensated 793by $t \ge -b(\vec p)$ while and non-negative value remains non-negative 794because $t \ge 0$. 795 796We need to show that this transformation does not eliminate any valid 797solutions and that it does not introduce any spurious solutions. 798Given a valid solution for the original problem, we need to find 799a non-negative value of $c'$ satisfying the constraints. 800If $b(\vec p) \ge 0$, we can take $t = 0$ so that 801$c' = c_j - t = c_j \ge 0$. 802If $b(\vec p) < 0$, we can take $t = -b(\vec p)$. 803Since $r = b(\vec p) + c_j - f \ge 0$ and $f \ge 0$, we have 804$c' = c_j + b(\vec p) \ge 0$. 805Note that these choices amount to plugging in 806$t = \left|-b(\vec p)\right|_+ = \max(-b(\vec p), 0)$. 807Conversely, given a solution to the new problem, we need to find 808a non-negative value of $c_j$, but this is easy since $c_j = c' + t$ 809and both of these are non-negative. 810 811Plugging 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, 813it may be better to provide a direct representation for such 814expressions in the internal representation of sets and relations 815or at least in a quast-like output format. 816 817\section{Coalescing}\label{s:coalescing} 818 819See \shortciteN{Verdoolaege2009isl}, for now. 820More details will be added later. 821 822\section{Transitive Closure} 823 824\subsection{Introduction} 825 826\begin{definition}[Power of a Relation] 827Let $R \in \Z^n \to 2^{\Z^{d+d}}$ be a relation and 828$k \in \Z_{\ge 1}$ 829a positive number, then power $k$ of relation $R$ is defined as 830\begin{equation} 831\label{eq:transitive:power} 832R^k \coloneqq 833\begin{cases} 834R & \text{if $k = 1$} 835\\ 836R \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] 843Let $R \in \Z^n \to 2^{\Z^{d+d}}$ be a relation, 844then the transitive closure $R^+$ of $R$ is the union 845of all positive powers of $R$, 846$$ 847R^+ \coloneqq \bigcup_{k \ge 1} R^k 848. 849$$ 850\end{definition} 851Alternatively, the transitive closure may be defined 852inductively as 853\begin{equation} 854\label{eq:transitive:inductive} 855R^+ \coloneqq R \cup \left(R \circ R^+\right) 856. 857\end{equation} 858 859Since the transitive closure of a polyhedral relation 860may no longer be a polyhedral relation \shortcite{Kelly1996closure}, 861we can, in the general case, only compute an approximation 862of the transitive closure. 863Whereas \shortciteN{Kelly1996closure} compute underapproximations, 864we, like \shortciteN{Beletska2009}, compute overapproximations. 865That is, given a relation $R$, we will compute a relation $T$ 866such that $R^+ \subseteq T$. Of course, we want this approximation 867to be as close as possible to the actual transitive closure 868$R^+$ and we want to detect the cases where the approximation is 869exact, i.e., where $T = R^+$. 870 871For computing an approximation of the transitive closure of $R$, 872we follow the same general strategy as \shortciteN{Beletska2009} 873and first compute an approximation of $R^k$ for $k \ge 1$ and then project 874out the parameter $k$ from the resulting relation. 875 876\begin{example} 877As a trivial example, consider the relation 878$R = \{\, x \to x + 1 \,\}$. The $k$th power of this map 879for arbitrary $k$ is 880$$ 881R^k = k \mapsto \{\, x \to x + k \mid k \ge 1 \,\} 882. 883$$ 884The transitive closure is then 885$$ 886\begin{aligned} 887R^+ & = \{\, 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 898There are some special cases where the computation of $R^k$ is very easy. 899One such case is that where $R$ does not compose with itself, 900i.e., $R \circ R = \emptyset$ or $\domain R \cap \range R = \emptyset$. 901In this case, $R^k$ is only non-empty for $k=1$ where it is equal 902to $R$ itself. 903 904In general, it is impossible to construct a closed form 905of $R^k$ as a polyhedral relation. 906We will therefore need to make some approximations. 907As a first approximations, we will consider each of the basic 908relations in $R$ as simply adding one or more offsets to a domain element 909to arrive at an image element and ignore the fact that some of these 910offsets may only be applied to some of the domain elements. 911That is, we will only consider the difference set $\Delta\,R$ of the relation. 912In particular, we will first construct a collection $P$ of paths 913that move through 914a total of $k$ offsets and then intersect domain and range of this 915collection with those of $R$. 916That is, 917\begin{equation} 918\label{eq:transitive:approx} 919K = P \cap \left(\domain R \to \range R\right) 920, 921\end{equation} 922with 923\begin{equation} 924\label{eq:transitive:path} 925P = \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} 932and with $\Delta_i$ the basic sets that compose 933the difference set $\Delta\,R$. 934Note that the number of basic sets $\Delta_i$ need not be 935the same as the number of basic relations in $R$. 936Also note that since addition is commutative, it does not 937matter in which order we add the offsets and so we are allowed 938to group them as we did in \eqref{eq:transitive:path}. 939 940If all the $\Delta_i$s are singleton sets 941$\Delta_i = \{\, \vec \delta_i \,\}$ with $\vec \delta_i \in \Z^d$, 942then \eqref{eq:transitive:path} simplifies to 943\begin{equation} 944\label{eq:transitive:singleton} 945P = \{\, \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} 952and then the approximation computed in \eqref{eq:transitive:approx} 953is essentially the same as that of \shortciteN{Beletska2009}. 954If some of the $\Delta_i$s are not singleton sets or if 955some of $\vec \delta_i$s are parametric, then we need 956to resort to further approximations. 957 958To ease both the exposition and the implementation, we will for 959the remainder of this section work with extended offsets 960$\Delta_i' = \Delta_i \times \{\, 1 \,\}$. 961That is, each offset is extended with an extra coordinate that is 962set equal to one. The paths constructed by summing such extended 963offsets have the length encoded as the difference of their 964final coordinates. The path $P'$ can then be decomposed into 965paths $P_i'$, one for each $\Delta_i$, 966\begin{equation} 967\label{eq:transitive:decompose} 968P' = \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} 978with 979$$ 980P_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$$ 986Note that each $P_i'$ contains paths of length at least one. 987We therefore need to take the union with the identity relation 988when composing the $P_i'$s to allow for paths that do not contain 989any offsets from one or more $\Delta_i'$. 990The path that consists of only identity relations is removed 991by imposing the constraint $y_{d+1} - x_{d+1} > 0$. 992Taking the union with the identity relation means that 993that the relations we compose in \eqref{eq:transitive:decompose} 994each consist of two basic relations. If there are $m$ 995disjuncts in the input relation, then a direct application 996of the composition operation may therefore result in a relation 997with $2^m$ disjuncts, which is prohibitively expensive. 998It is therefore crucial to apply coalescing (\autoref{s:coalescing}) 999after each composition. 1000 1001Let us now consider how to compute an overapproximation of $P_i'$. 1002Those that correspond to singleton $\Delta_i$s are grouped together 1003and handled as in \eqref{eq:transitive:singleton}. 1004Note that this is just an optimization. The procedure described 1005below would produce results that are at least as accurate. 1006For simplicity, we first assume that no constraint in $\Delta_i'$ 1007involves any existentially quantified variables. 1008We will return to existentially quantified variables at the end 1009of this section. 1010Without existentially quantified variables, we can classify 1011the constraints of $\Delta_i'$ as follows 1012\begin{enumerate} 1013\item non-parametric constraints 1014\begin{equation} 1015\label{eq:transitive:non-parametric} 1016A_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} 1021B_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} 1026A_3 \vec x + B_3 \vec s + \vec c_3 \geq \vec 0 1027\end{equation} 1028such 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$$ 1036A_4 \vec x + B_4 \vec s + \vec c_4 \geq \vec 0 1037$$ 1038such 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} 1045We will use the following approximation $Q_i$ for $P_i'$: 1046\begin{equation} 1047\label{eq:transitive:Q} 1048\begin{aligned} 1049Q_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& 1057A_1 \vec f + k \vec c_1 \geq \vec 0 1058\wedge 1059B_2 \vec s + \vec c_2 \geq \vec 0 1060\wedge 1061A_3 \vec f + B_3 \vec s + \vec c_3 \geq \vec 0 1062\,\} 1063. 1064\end{aligned} 1065\end{equation} 1066To prove that $Q_i$ is indeed an overapproximation of $P_i'$, 1067we 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)$ 1069we have that 1070$(\vec f, k)$ satisfies the constraints in \eqref{eq:transitive:Q}. 1071If $\Delta_i(\vec s)$ is non-empty, then $\vec s$ must satisfy 1072the constraints in \eqref{eq:transitive:parametric}. 1073Each element $(\vec f, k) \in k \, \Delta_i'(\vec s)$ is a sum 1074of $k$ elements $(\vec f_j, 1)$ in $\Delta_i'(\vec s)$. 1075Each of these elements satisfies the constraints in 1076\eqref{eq:transitive:non-parametric}, i.e., 1077$$ 1078\left[ 1079\begin{matrix} 1080A_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$$ 1091The sum of these elements therefore satisfies the same set of inequalities, 1092i.e., $A_1 \vec f + k \vec c_1 \geq \vec 0$. 1093Finally, the constraints in \eqref{eq:transitive:mixed} are such 1094that for any $\vec s$ in the parameter domain of $\Delta$, 1095we have $-\vec r(\vec s) \coloneqq B_3 \vec s + \vec c_3 \le \vec 0$, 1096i.e., $A_3 \vec f_j \ge \vec r(\vec s) \ge \vec 0$ 1097and therefore also $A_3 \vec f \ge \vec r(\vec s)$. 1098Note that if there are no mixed constraints and if the 1099rational 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\,\}$, 1101has integer vertices, then the approximation is exact, i.e., 1102$Q_i = P_i'$. In this case, the vertices of $\Delta'_i(\vec s)$ 1103generate the rational cone 1104$\{\, \vec x' \in \Q^{d+1} \mid \left[ 1105\begin{matrix} 1106A_1 & \vec c_1 1107\end{matrix} 1108\right] \vec x' \,\}$ and therefore $\Delta'_i(\vec s)$ is 1109a Hilbert basis of this cone \shortcite[Theorem~16.4]{Schrijver1986}. 1110 1111Note however that, as pointed out by \shortciteN{DeSmet2010personal}, 1112if there \emph{are} any mixed constraints, then the above procedure may 1113not compute the most accurate affine approximation of 1114$k \, \Delta_i(\vec s)$ with $k \ge 1$. 1115In particular, we only consider the negative mixed constraints that 1116happen to appear in the description of $\Delta_i(\vec s)$, while we 1117should instead consider \emph{all} valid such constraints. 1118It is also sufficient to consider those constraints because any 1119constraint that is valid for $k \, \Delta_i(\vec s)$ is also 1120valid for $1 \, \Delta_i(\vec s) = \Delta_i(\vec s)$. 1121Take therefore any constraint 1122$\spv a x + \spv b s + c \ge 0$ valid for $\Delta_i(\vec s)$. 1123This constraint is also valid for $k \, \Delta_i(\vec s)$ iff 1124$k \, \spv a x + \spv b s + c \ge 0$. 1125If $\spv b s + c$ can attain any positive value, then $\spv a x$ 1126may be negative for some elements of $\Delta_i(\vec s)$. 1127We then have $k \, \spv a x < \spv a x$ for $k > 1$ and so the constraint 1128is not valid for $k \, \Delta_i(\vec s)$. 1129We therefore need to impose $\spv b s + c \le 0$ for all values 1130of $\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 1132constraint of $\Delta_i(\vec s)$. That is, $(\vec b, c)$ are the opposites 1133of the coefficients of a valid constraint of $\Delta_i(\vec s)$. 1134The approximation of $k \, \Delta_i(\vec s)$ can therefore be obtained 1135using three applications of Farkas' lemma. The first obtains the coefficients 1136of constraints valid for $\Delta_i(\vec s)$. The second obtains 1137the coefficients of constraints valid for the projection of $\Delta_i(\vec s)$ 1138onto the parameters. The opposite of the second set is then computed 1139and intersected with the first set. The result is the set of coefficients 1140of constraints valid for $k \, \Delta_i(\vec s)$. A final application 1141of Farkas' lemma is needed to obtain the approximation of 1142$k \, \Delta_i(\vec s)$ itself. 1143 1144\begin{example} 1145Consider the relation 1146$$ 1147n \to \{\, (x, y) \to (1 + x, 1 - n + y) \mid n \ge 2 \,\} 1148. 1149$$ 1150The difference set of this relation is 1151$$ 1152\Delta = n \to \{\, (1, 1 - n) \mid n \ge 2 \,\} 1153. 1154$$ 1155Using our approach, we would only consider the mixed constraint 1156$y - 1 + n \ge 0$, leading to the following approximation of the 1157transitive closure: 1158$$ 1159n \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$$ 1162If, instead, we apply Farkas's lemma to $\Delta$, i.e., 1163\begin{verbatim} 1164D := [n] -> { [1, 1 - n] : n >= 2 }; 1165CD := coefficients D; 1166CD; 1167\end{verbatim} 1168we 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} 1173The pure-parametric constraints valid for $\Delta$, 1174\begin{verbatim} 1175P := { [a,b] -> [] }(D); 1176CP := coefficients P; 1177CP; 1178\end{verbatim} 1179are 1180\begin{verbatim} 1181{ rat: coefficients[[c_cst, c_n] -> []] : c_n >= 0 and 2c_n >= -c_cst } 1182\end{verbatim} 1183Negating these coefficients and intersecting with \verb+CD+, 1184\begin{verbatim} 1185NCP := { rat: coefficients[[a,b] -> []] 1186 -> coefficients[[-a,-b] -> []] }(CP); 1187CK := wrap((unwrap CD) * (dom (unwrap NCP))); 1188CK; 1189\end{verbatim} 1190we 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} 1195The approximation for $k\,\Delta$, 1196\begin{verbatim} 1197K := solutions CK; 1198K; 1199\end{verbatim} 1200is then 1201\begin{verbatim} 1202[n] -> { rat: [i0, i1] : i1 <= -i0 and i0 >= 1 and i1 <= 2 - n - i0 } 1203\end{verbatim} 1204Finally, the computed approximation for $R^+$, 1205\begin{verbatim} 1206T := unwrap({ [dx,dy] -> [[x,y] -> [x+dx,y+dy]] }(K)); 1207R := [n] -> { [x,y] -> [x+1,y+1-n] : n >= 2 }; 1208T := T * ((dom R) -> (ran R)); 1209T; 1210\end{verbatim} 1211is 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 1218Existentially quantified variables can be handled by 1219classifying them into variables that are uniquely 1220determined by the parameters, variables that are independent 1221of the parameters and others. The first set can be treated 1222as parameters and the second as variables. Constraints involving 1223the other existentially quantified variables are removed. 1224 1225\begin{example} 1226Consider the relation 1227$$ 1228R = 1229n \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$$ 1232The difference set of this relation is 1233$$ 1234\Delta = \Delta \, R = 1235n \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$$ 1238The existentially quantified variables can be defined in terms 1239of 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, 1249while $\alpha_1$ can be treated as a variable. 1250This in turn means that $7\alpha_0 = -2 + n$ can be treated as 1251a purely parametric constraint, while the other two constraints are 1252non-parametric. 1253The corresponding $Q$~\eqref{eq:transitive:Q} is therefore 1254$$ 1255\begin{aligned} 1256n \to \{\, (x,z) \to (y,w) \mid 1257\exists\, \alpha_0, \alpha_1, k, f : {} & 1258k \ge 1 \wedge 1259y = x + f \wedge 1260w = z + k \wedge {} \\ 1261& 12627\alpha_0 = -2 + n \wedge 12635\alpha_1 = -k + x \wedge 1264x \ge 6 k 1265\,\} 1266. 1267\end{aligned} 1268$$ 1269Projecting out the final coordinates encoding the length of the paths, 1270results in the exact transitive closure 1271$$ 1272R^+ = 1273n \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 1278The fact that we ignore some impure constraints clearly leads 1279to a loss of accuracy. In some cases, some of this loss can be recovered 1280by not considering the parameters in a special way. 1281That 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$$ 1289we 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$$ 1298The first $n$ coordinates of every element in $\Delta'$ are zero. 1299Projecting out these zero coordinates from $\Delta'$ is equivalent 1300to projecting out the parameters in $\Delta$. 1301The result is obviously a superset of $\Delta$, but all its constraints 1302are of type \eqref{eq:transitive:non-parametric} and they can therefore 1303all be used in the construction of $Q_i$. 1304 1305\begin{example} 1306Consider the relation 1307$$ 1308% [n] -> { [x, y] -> [1 + x, 1 - n + y] | n >= 2 } 1309R = n \to \{\, (x, y) \to (1 + x, 1 - n + y) \mid n \ge 2 \,\} 1310. 1311$$ 1312We have 1313$$ 1314\diff R = n \to \{\, (1, 1 - n) \mid n \ge 2 \,\} 1315$$ 1316and so, by treating the parameters in a special way, we obtain 1317the following approximation for $R^+$: 1318$$ 1319n \to \{\, (x, y) \to (x', y') \mid n \ge 2 \wedge y' \le 1 - n + y \wedge x' \ge 1 + x \,\} 1320. 1321$$ 1322If we consider instead 1323$$ 1324R' = \{\, (n, x, y) \to (n, 1 + x, 1 - n + y) \mid n \ge 2 \,\} 1325$$ 1326then 1327$$ 1328\diff R' = \{\, (0, 1, y) \mid y \le -1 \,\} 1329$$ 1330and we obtain the approximation 1331$$ 1332n \to \{\, (x, y) \to (x', y') \mid n \ge 2 \wedge x' \ge 1 + x \wedge y' \le x + y - x' \,\} 1333. 1334$$ 1335If we consider both $\diff R$ and $\diff R'$, then we obtain 1336$$ 1337n \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$$ 1340Note, however, that this is not the most accurate affine approximation that 1341can be obtained. That would be 1342$$ 1343n \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 1350The approximation $T$ for the transitive closure $R^+$ can be obtained 1351by projecting out the parameter $k$ from the approximation $K$ 1352\eqref{eq:transitive:approx} of the power $R^k$. 1353Since $K$ is an overapproximation of $R^k$, $T$ will also be an 1354overapproximation of $R^+$. 1355To check whether the results are exact, we need to consider two 1356cases depending on whether $R$ is {\em cyclic}, where $R$ is defined 1357to be cyclic if $R^+$ maps any element to itself, i.e., 1358$R^+ \cap \identity \ne \emptyset$. 1359If $R$ is acyclic, then the inductive definition of 1360\eqref{eq:transitive:inductive} is equivalent to its completion, 1361i.e., 1362$$ 1363R^+ = R \cup \left(R \circ R^+\right) 1364$$ 1365is a defining property. 1366Since $T$ is known to be an overapproximation, we only need to check 1367whether 1368$$ 1369T \subseteq R \cup \left(R \circ T\right) 1370. 1371$$ 1372This is essentially Theorem~5 of \shortciteN{Kelly1996closure}. 1373The only difference is that they only consider lexicographically 1374forward relations, a special case of acyclic relations. 1375 1376If, on the other hand, $R$ is cyclic, then we have to resort 1377to checking whether the approximation $K$ of the power is exact. 1378Note that $T$ may be exact even if $K$ is not exact, so the check 1379is sound, but incomplete. 1380To check exactness of the power, we simply need to check 1381\eqref{eq:transitive:power}. Since again $K$ is known 1382to be an overapproximation, we only need to check whether 1383$$ 1384\begin{aligned} 1385K'|_{y_{d+1} - x_{d+1} = 1} & \subseteq R' 1386\\ 1387K'|_{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$$ 1391where $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 1393lengths equal to 1. 1394 1395All that remains is to explain how to check the cyclicity of $R$. 1396Note that the exactness on the power is always sound, even 1397in the acyclic case, so we only need to be careful that we find 1398all cyclic cases. Now, if $R$ is cyclic, i.e., 1399$R^+ \cap \identity \ne \emptyset$, then, since $T$ is 1400an overapproximation of $R^+$, also 1401$T \cap \identity \ne \emptyset$. This in turn means 1402that $\Delta \, K'$ contains a point whose first $d$ coordinates 1403are zero and whose final coordinate is positive. 1404In the implementation we currently perform this test on $P'$ instead of $K'$. 1405Note that if $R^+$ is acyclic and $T$ is not, then the approximation 1406is clearly not exact and the approximation of the power $K$ 1407will not be exact either. 1408 1409\subsection{Decomposing $R$ into strongly connected components} 1410 1411If the input relation $R$ is a union of several basic relations 1412that can be partially ordered 1413then the accuracy of the approximation may be improved by computing 1414an approximation of each strongly connected components separately. 1415For example, if $R = R_1 \cup R_2$ and $R_1 \circ R_2 = \emptyset$, 1416then we know that any path that passes through $R_2$ cannot later 1417pass through $R_1$, i.e., 1418\begin{equation} 1419\label{eq:transitive:components} 1420R^+ = R_1^+ \cup R_2^+ \cup \left(R_2^+ \circ R_1^+\right) 1421. 1422\end{equation} 1423We can therefore compute (approximations of) transitive closures 1424of $R_1$ and $R_2$ separately. 1425Note, however, that the condition $R_1 \circ R_2 = \emptyset$ 1426is actually too strong. 1427If $R_1 \circ R_2$ is a subset of $R_2 \circ R_1$ 1428then we can reorder the segments 1429in any path that moves through both $R_1$ and $R_2$ to 1430first move through $R_1$ and then through $R_2$. 1431 1432This idea can be generalized to relations that are unions 1433of more than two basic relations by constructing the 1434strongly connected components in the graph with as vertices 1435the 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. 1437That is, there is an edge from $R_i$ to $R_j$ iff 1438\begin{equation} 1439\label{eq:transitive:edge} 1440R_i \circ R_j 1441\not\subseteq 1442R_j \circ R_i 1443. 1444\end{equation} 1445The components can be obtained from the graph by applying 1446Tarjan's algorithm \shortcite{Tarjan1972}. 1447 1448In practice, we compute the (extended) powers $K_i'$ of each component 1449separately and then compose them as in \eqref{eq:transitive:decompose}. 1450Note, however, that in this case the order in which we apply them is 1451important and should correspond to a topological ordering of the 1452strongly connected components. Simply applying Tarjan's 1453algorithm will produce topologically sorted strongly connected components. 1454The graph on which Tarjan's algorithm is applied is constructed on-the-fly. 1455That is, whenever the algorithm checks if there is an edge between 1456two vertices, we evaluate \eqref{eq:transitive:edge}. 1457The exactness check is performed on each component separately. 1458If the approximation turns out to be inexact for any of the components, 1459then the entire result is marked inexact and the exactness check 1460is skipped on the components that still need to be handled. 1461 1462It should be noted that \eqref{eq:transitive:components} 1463is only valid for exact transitive closures. 1464If overapproximations are computed in the right hand side, then the result will 1465still be an overapproximation of the left hand side, but this result 1466may not be transitively closed. If we only separate components based 1467on the condition $R_i \circ R_j = \emptyset$, then there is no problem, 1468as this condition will still hold on the computed approximations 1469of the transitive closures. If, however, we have exploited 1470\eqref{eq:transitive:edge} during the decomposition and if the 1471result turns out not to be exact, then we check whether 1472the result is transitively closed. If not, we recompute 1473the transitive closure, skipping the decomposition. 1474Note that testing for transitive closedness on the result may 1475be fairly expensive, so we may want to make this check 1476configurable. 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} 1498Consider the relation in example {\tt closure4} that comes with 1499the Omega calculator~\shortcite{Omega_calc}, $R = R_1 \cup R_2$, 1500with 1501$$ 1502\begin{aligned} 1503R_1 & = \{\, (x,y) \to (x,y+1) \mid 1 \le x,y \le 10 \,\} 1504\\ 1505R_2 & = \{\, (x,y) \to (x+1,y) \mid 1 \le x \le 20 \wedge 5 \le y \le 15 \,\} 1506. 1507\end{aligned} 1508$$ 1509This relation is shown graphically in \autoref{f:closure4}. 1510We have 1511$$ 1512\begin{aligned} 1513R_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\\ 1516R_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$$ 1521Clearly, $R_1 \circ R_2 \subseteq R_2 \circ R_1$ and so 1522$$ 1523\left( 1524R_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} 1578Consider the relation on the right of \shortciteN[Figure~2]{Beletska2009}, 1579reproduced in \autoref{f:decomposition}. 1580The relation can be described as $R = R_1 \cup R_2 \cup R_3$, 1581with 1582$$ 1583\begin{aligned} 1584R_1 &= n \mapsto \{\, (i,j) \to (i+3,j) \mid 1585i \le 2 j - 4 \wedge 1586i \le n - 3 \wedge 1587j \le 2 i - 1 \wedge 1588j \le n \,\} 1589\\ 1590R_2 &= n \mapsto \{\, (i,j) \to (i,j+3) \mid 1591i \le 2 j - 1 \wedge 1592i \le n \wedge 1593j \le 2 i - 4 \wedge 1594j \le n - 3 \,\} 1595\\ 1596R_3 &= n \mapsto \{\, (i,j) \to (i+1,j+1) \mid 1597i \le 2 j - 1 \wedge 1598i \le n - 1 \wedge 1599j \le 2 i - 1 \wedge 1600j \le n - 1\,\} 1601. 1602\end{aligned} 1603$$ 1604The figure shows this relation for $n = 7$. 1605Both 1606$R_3 \circ R_1 \subseteq R_1 \circ R_3$ 1607and 1608$R_3 \circ R_2 \subseteq R_2 \circ R_3$, 1609which the reader can verify using the {\tt iscc} calculator: 1610\begin{verbatim} 1611R1 := [n] -> { [i,j] -> [i+3,j] : i <= 2 j - 4 and i <= n - 3 and 1612 j <= 2 i - 1 and j <= n }; 1613R2 := [n] -> { [i,j] -> [i,j+3] : i <= 2 j - 1 and i <= n and 1614 j <= 2 i - 4 and j <= n - 3 }; 1615R3 := [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. 1621For the other two basic relations, we have both 1622$R_2 \circ R_1 \not\subseteq R_1 \circ R_2$ 1623and 1624$R_1 \circ R_2 \not\subseteq R_2 \circ R_1$ 1625and so $R_1$ and $R_2$ form a strongly connected component. 1626By computing the power of $R_3$ and $R_1 \cup R_2$ separately 1627and composing the results, the power of $R$ can be computed exactly 1628using \eqref{eq:transitive:singleton}. 1629As explained by \shortciteN{Beletska2009}, applying the same formula 1630to $R$ directly, without a decomposition, would result in 1631an overapproximation of the power. 1632\end{example} 1633 1634\subsection{Partitioning the domains and ranges of $R$} 1635 1636The algorithm of \autoref{s:power} assumes that the input relation $R$ 1637can be treated as a union of translations. 1638This is a reasonable assumption if $R$ maps elements of a given 1639abstract domain to the same domain. 1640However, if $R$ is a union of relations that map between different 1641domains, then this assumption no longer holds. 1642In particular, when an entire dependence graph is encoded 1643in a single relation, as is done by, e.g., 1644\shortciteN[Section~6.1]{Barthou2000MSE}, then it does not make 1645sense to look at differences between iterations of different domains. 1646Now, arguably, a modified Floyd-Warshall algorithm should 1647be applied to the dependence graph, as advocated by 1648\shortciteN{Kelly1996closure}, with the transitive closure operation 1649only being applied to relations from a given domain to itself. 1650However, it is also possible to detect disjoint domains and ranges 1651and 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\SetAlgoVlined 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 1682Let the input relation $R$ be a union of $m$ basic relations $R_i$. 1683Let $D_{2i}$ be the domains of $R_i$ and $D_{2i+1}$ the ranges of $R_i$. 1684The first step is to group overlapping $D_j$ until a partition is 1685obtained. If the resulting partition consists of a single part, 1686then we continue with the algorithm of \autoref{s:power}. 1687Otherwise, we apply Floyd-Warshall on the graph with as vertices 1688the parts of the partition and as edges the $R_i$ attached to 1689the appropriate pairs of vertices. 1690In particular, let there be $n$ parts $P_k$ in the partition. 1691We construct $n^2$ relations 1692$$ 1693R_{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$$ 1697apply \autoref{a:Floyd} and return the union of all resulting 1698$R_{pq}$ as the transitive closure of $R$. 1699Each iteration of the $r$-loop in \autoref{a:Floyd} updates 1700all relations $R_{pq}$ to include paths that go from $p$ to $r$, 1701possibly stay there for a while, and then go from $r$ to $q$. 1702Note that paths that ``stay in $r$'' include all paths that 1703pass through earlier vertices since $R_{rr}$ itself has been updated 1704accordingly in previous iterations of the outer loop. 1705In principle, it would be sufficient to use the $R_{pr}$ 1706and $R_{rq}$ computed in the previous iteration of the 1707$r$-loop in Line~\ref{l:Floyd:update}. 1708However, from an implementation perspective, it is easier 1709to allow either or both of these to have been updated 1710in the same iteration of the $r$-loop. 1711This may result in duplicate paths, but these can usually 1712be removed by coalescing (\autoref{s:coalescing}) the result of the union 1713in Line~\ref{l:Floyd:update}, which should be done in any case. 1714The transitive closure in Line~\ref{l:Floyd:closure} 1715is performed using a recursive call. This recursive call 1716includes the partitioning step, but the resulting partition will 1717usually be a singleton. 1718The result of the recursive call will either be exact or an 1719overapproximation. The final result of Floyd-Warshall is therefore 1720also 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} 1744Consider the relation on the right of Figure~1 of 1745\shortciteN{Beletska2009}, 1746reproduced in \autoref{f:COCOA:1}. 1747This 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$$ 1755Note that the domain of the upward relation overlaps with the range 1756of the rightward relation and vice versa, but that the domain 1757of neither relation overlaps with its own range or the domain of 1758the other relation. 1759The 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}, 1761respectively. 1762Initially, we have 1763$$ 1764\begin{aligned} 1765R_{00} & = \emptyset 1766\\ 1767R_{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\\ 1771R_{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\\ 1774R_{11} & = \emptyset 1775. 1776\end{aligned} 1777$$ 1778In 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., 1781the dashed arrow in the figure. 1782This new $R_{11}$ is obviously transitively closed, so it is not 1783changed in the second iteration and it does not have an effect 1784on $R_{01}$ and $R_{10}$. However, $R_{00}$ is updated to 1785include $R_{10} \circ R_{01}$, i.e., the dotted arrow in the figure. 1786The 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 1793In some cases it is possible and useful to compute the transitive closure 1794of union of basic relations incrementally. In particular, 1795if $R$ is a union of $m$ basic maps, 1796$$ 1797R = \bigcup_j R_j 1798, 1799$$ 1800then we can pick some $R_i$ and compute the transitive closure of $R$ as 1801\begin{equation} 1802\label{eq:transitive:incremental} 1803R^+ = R_i^+ \cup 1804\left( 1805\bigcup_{j \ne i} 1806R_i^* \circ R_j \circ R_i^* 1807\right)^+ 1808. 1809\end{equation} 1810For this approach to be successful, it is crucial that each 1811of the disjuncts in the argument of the second transitive 1812closure in \eqref{eq:transitive:incremental} be representable 1813as a single basic relation, i.e., without a union. 1814If this condition holds, then by using \eqref{eq:transitive:incremental}, 1815the number of disjuncts in the argument of the transitive closure 1816can be reduced by one. 1817Now, $R_i^* = R_i^+ \cup \identity$, but in some cases it is possible 1818to relax the constraints of $R_i^+$ to include part of the identity relation, 1819say on domain $D$. We will use the notation 1820${\cal C}(R_i,D) = R_i^+ \cup \identity_D$ to represent 1821this 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 1824the value $0$ in \eqref{eq:transitive:Q} and by using 1825$$ 1826P \cap \left(D \to D\right) 1827$$ 1828instead of \eqref{eq:transitive:approx}. 1829Typically, $D$ will be a strict superset of both $\domain R_i$ 1830and $\range R_i$. We therefore need to check that domain 1831and range of the transitive closure are part of ${\cal C}(R_i,D)$, 1832i.e., the part that results from the paths of positive length ($k \ge 1$), 1833are equal to the domain and range of $R_i$. 1834If not, then the incremental approach cannot be applied for 1835the given choice of $R_i$ and $D$. 1836 1837In order to be able to replace $R^*$ by ${\cal C}(R_i,D)$ 1838in \eqref{eq:transitive:incremental}, $D$ should be chosen 1839to include both $\domain R$ and $\range R$, i.e., such 1840that $\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 1843they use $D = \domain R \cup \range R$. 1844Now, this expression of $D$ contains a union, so it not directly usable. 1845\shortciteN{Kelly1996closure} do not explain how they avoid this union. 1846Apparently, in their implementation, 1847they are using the convex hull of $\domain R \cup \range R$ 1848or at least an approximation of this convex hull. 1849We use the simple hull (\autoref{s:simple hull}) of $\domain R \cup \range R$. 1850 1851It is also possible to use a domain $D$ that does {\em not\/} 1852include $\domain R \cup \range R$, but then we have to 1853compose with ${\cal C}(R_i,D)$ more selectively. 1854In 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} 1860and, 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} 1866then we can refine \eqref{eq:transitive:incremental} to 1867$$ 1868R_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\!\!\!\!\! 1887R_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\!\!\!\!\! 1894R_j 1895\right) 1896\right)^+ 1897. 1898$$ 1899If only property~\eqref{eq:transitive:right} holds, 1900we can use 1901$$ 1902R_i^+ \cup 1903\left( 1904\left( 1905R_i^+ \cup \identity 1906\right) 1907\circ 1908\left( 1909\left( 1910\bigcup_{\shortstack{$\scriptstyle\domain R_j \subseteq D $}} 1911R_j \circ {\cal C} 1912\right) 1913\cup 1914\left( 1915\bigcup_{\shortstack{$\scriptstyle\domain R_j \cap \range R_i = \emptyset$}} 1916\!\!\!\!\! 1917R_j 1918\right) 1919\right)^+ 1920\right) 1921, 1922$$ 1923while if only property~\eqref{eq:transitive:left} holds, 1924we can use 1925$$ 1926R_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\!\!\!\!\! 1937R_j 1938\right) 1939\right)^+ 1940\circ 1941\left( 1942R_i^+ \cup \identity 1943\right) 1944\right) 1945. 1946$$ 1947 1948It should be noted that if we want the result of the incremental 1949approach to be transitively closed, then we can only apply it 1950if all of the transitive closure operations involved are exact. 1951If, say, the second transitive closure in \eqref{eq:transitive:incremental} 1952contains extra elements, then the result does not necessarily contain 1953the composition of these extra elements with powers of $R_i$. 1954 1955\subsection{An {\tt Omega}-like implementation} 1956 1957While the main algorithm of \shortciteN{Kelly1996closure} is 1958designed to compute and underapproximation of the transitive closure, 1959the authors mention that they could also compute overapproximations. 1960In this section, we describe our implementation of an algorithm 1961that is based on their ideas. 1962Note that the {\tt Omega} library computes underapproximations 1963\shortcite[Section 6.4]{Omega_lib}. 1964 1965The main tool is Equation~(2) of \shortciteN{Kelly1996closure}. 1966The 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$$ 1975where $p$ ranges over the dimensions and $\vec L$, $\vec U$ and 1976$\vec M$ are constant integer vectors. The elements of $\vec U$ 1977may be $\infty$, meaning that there is no upper bound corresponding 1978to that element, and similarly for $\vec L$. 1979Such an overapproximation can be obtained by computing strides, 1980lower and upper bounds on the difference set $\Delta \, R$. 1981The 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 : 1985k \ge 1 \wedge 1986k \, \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} 1992The domain and range of this transitive closure are then 1993intersected with those of the input relation. 1994This is a special case of the algorithm in \autoref{s:power}. 1995 1996In their algorithm for computing lower bounds, the authors 1997use the above algorithm as a substep on the disjuncts in the relation. 1998At the end, they say 1999\begin{quote} 2000If an upper bound is required, it can be calculated in a manner 2001similar to that of a single conjunct [sic] relation. 2002\end{quote} 2003Presumably, the authors mean that a ``d-form'' approximation 2004of the whole input relation should be used. 2005However, the accuracy can be improved by also trying to 2006apply the incremental technique from the same paper, 2007which is explained in more detail in \autoref{s:incremental}. 2008In this case, ${\cal C}(R_i,D)$ can be obtained by 2009allowing the value zero for $k$ in \eqref{eq:omega}, 2010i.e., by computing 2011$$ 2012\{\, \vec i \to \vec j \mid \exists \vec \alpha, k : 2013k \ge 0 \wedge 2014k \, \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$$ 2020In our implementation we take as $D$ the simple hull 2021(\autoref{s:simple hull}) of $\domain R \cup \range R$. 2022To determine whether it is safe to use ${\cal C}(R_i,D)$, 2023we 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$ 2026the condition 2027$$ 2028\left({\cal C}(R_i,D) - R_i^+\right) 2029\circ 2030R_j 2031\circ 2032\left({\cal C}(R_i,D) - R_i^+\right) 2033= 2034R_j 2035$$ 2036holds. 2037