2050 lines · plain
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 sets5$S = \bigcup_i S_i$, each of which can be represented using affine6constraints7$$8S_i : \Z^n \to 2^{\Z^d} : \vec s \mapsto9S_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 set23$$\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 relations29$R = \bigcup_i R_i$ of type30$\Z^n \to 2^{\Z^{d_1+d_2}}$,31each of which can be represented using affine32constraints33$$34R_i = \vec s \mapsto35R_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 set50$$\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 set56$$\domain R \coloneqq \vec s \mapsto57\{\, \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 set66$$67\range R \coloneqq \vec s \mapsto68\{\, \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}}$ and76$S \in \Z^n \to 2^{\Z^{d_2+d_3}}$ be two relations,77then the composition of78$R$ and $S$ is defined as79$$80S \circ R \coloneqq81\vec s \mapsto82\{\, \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) \wedge85\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 set94of differences between image elements and the corresponding95domain elements,96$$97\diff R \coloneqq98\vec s \mapsto99\{\, \vec \delta \in \Z^{d} \mid \exists \vec x \to \vec y \in R :100\vec \delta = \vec y - \vec x101\,\}102$$103\end{definition}104 105\section{Simple Hull}\label{s:simple hull}106 107It is sometimes useful to have a single108basic set or basic relation that contains a given set or relation.109For rational sets, the obvious choice would be to compute the110(rational) convex hull. For integer sets, the obvious choice111would be the integer hull.112However, {\tt isl} currently does not support an integer hull operation113and even if it did, it would be fairly expensive to compute.114The convex hull operation is supported, but it is also fairly115expensive 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 overapproximation120and it is defined as the (inclusion-wise) smallest basic set121that is described by constraints that are translates of122the constraints in the input set.123This means that the simple hull is relatively cheap to compute124and that the number of constraints in the simple hull is no125larger than the number of constraints in the input.126\begin{definition}[Simple Hull of a Set]127The {\em simple hull} of a set128$S = \bigcup_{1 \le i \le v} S_i$, with129$$130S : \Z^n \to 2^{\Z^d} : \vec s \mapsto131S(\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 set137$$138H : \Z^n \to 2^{\Z^d} : \vec s \mapsto139S(\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 0143\,\right\}144,145$$146with $\vec K_i$ the (component-wise) smallest non-negative integer vectors147such that $S \subseteq H$.148\end{definition}149The $\vec K_i$ can be obtained by solving a number of150LP problems, one for each element of each $\vec K_i$.151If any LP problem is unbounded, then the corresponding constraint152is 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 quantified162variables. The latter operation has been used for counting elements163in sets involving such variables164\parencite{BouletRe98,Verdoolaege2005experiences} and lies at the core165of the internal representation of {\tt isl}.166 167Parametric integer programming was first implemented in \texttt{PipLib}.168An alternative method for parametric integer programming169was later implemented in {\tt barvinok} \cite{barvinok-0.22}.170This method is not based on Feautrier's algorithm, but on rational171generating functions \cite{Woods2003short} and was inspired by the172``digging'' technique of \textcite{DeLoera2004Three} for solving173non-parametric integer programming problems.174 175In the following sections, we briefly recall the dual simplex176method combined with Gomory cuts and describe some extensions177and optimizations. The main algorithm is applied to a matrix178data structure known as a tableau. In case of parametric problems,179there are two tableaus, one for the main problem and one for180the 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 representation187as that used by its incremental LP solver based on the \emph{primal}188simplex method. The implementation of this LP solver is based189on 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 implemented192as a sparse matrix, but neither {\tt Simplify} nor the current193implementation 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 relationship197between the variables $\vec x$ and non-negative variables198$\vec y = A \vec x + \vec b$ corresponding to the constraints.199The initial tableau contains $\begin{pmatrix}200\vec b & A201\end{pmatrix}$ and expresses the constraints $\vec y$ in the rows in terms202of the variables $\vec x$ in the columns. The main operation defined203on 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 assigned208zero and each row variable is assigned the constant term of the row.209This sample value represents a valid solution if each constraint variable210is assigned a non-negative value, i.e., if the constant terms of211rows corresponding to constraints are all non-negative.212 213The dual simplex method starts from an initial sample value that214may be invalid, but that is known to be (lexicographically) no215greater than any solution, and gradually increments this sample value216through pivoting until a valid solution is obtained.217In particular, each pivot exchanges a row variable218$r = -n + \sum_i a_i \, c_i$ with negative219sample 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 tableau228is (implicitly) extended with extra rows defining the original229variables in terms of the column variables.230If we assume that all variables are non-negative, then we know231that the zero vector is no greater than the minimal solution and232then 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 positive247and 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 row251is 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 term255may depend on the parameters. Each time the constant term of a constraint row256changes, we therefore need to check whether the new term can attain257negative and/or positive values over the current set of possible258parameter values, i.e., the context.259If all these terms can only attain non-negative values, the current260state of the tableau represents a solution. If one of the terms261can 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 positive264and negative values and split the context into a part where265it only attains non-negative values and a part where it only attains266negative values.267 268\subsection{Gomory Cuts}269 270The solution found by the dual simplex method may have271non-integral coordinates. If so, some rational solutions272(including the current sample value), can be cut off by273applying a (parametric) Gomory cut.274Let $r = b(\vec p) + \sp {\vec a} {\vec c}$ be the row275corresponding to the first non-integral coordinate of $\vec x$,276with $b(\vec p)$ the constant term, an affine expression in the277parameters $\vec p$, i.e., $b(\vec p) = \sp {\vec f} {\vec p} + g$.278Note that only row variables can attain279non-integral values as the sample value of the column variables is zero.280Consider the expression281$b(\vec p) - \ceil{b(\vec p)} + \sp {\fract{\vec a}} {\vec c}$,282with $\ceil\cdot$ the ceiling function and $\fract\cdot$ the283fractional part. This expression is negative at the sample value284since $\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 integral286value of $r$ and $\vec c \ge 0$, the expression is non-negative287because $b(\vec p) - \ceil{b(\vec p)} > -1$.288Imposing this expression to be non-negative therefore does not289invalidate any integral solutions, while it does cut away the current290fractional sample value. To be able to formulate this constraint,291a new variable $q = \floor{-b(\vec p)} = - \ceil{b(\vec p)}$ is added292to the context. This integral variable is uniquely defined by the constraints293$0 \le -d \, b(\vec p) - d \, q \le d - 1$, with $d$ the common294denominator of $\vec f$ and $g$. In practice, the variable295$q' = \floor{\sp {\fract{-f}} {\vec p} + \fract{-g}}$ is used instead296and the coefficients of the new constraint are adjusted accordingly.297The sign of the constant term of this new constraint need not be determined298as it is non-positive by construction.299When several of these extra context variables are added, it is important300to 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 from307sample value $\vec x = \vec 0$ and $\vec c$ is assumed to be non-negative308during 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 be311an arbitrarily large positive number. Instead of looking for the312lexicographically minimal value of $\vec x$, we search instead313for the lexicographically minimal value of $\vec x' = \vec M + \vec x$.314The sample value $\vec x' = \vec 0$ of the initial tableau then315corresponds to $\vec x = -\vec M$, which is clearly not greater than316any potential solution. The sign of the constant term of a row317is determined lexicographically, with the coefficient of $M$ considered318first. That is, if the coefficient of $M$ is not zero, then its sign319is the sign of the entire term. Otherwise, the sign is determined320by the remaining affine expression in the parameters.321If the original problem has a bounded optimum, then the final sample322value will be of the form $\vec M + \vec v$ and the optimal value323of the original problem is then $\vec v$.324Maximization problems can be handled in a similar way by computing325the minimum of $\vec M - \vec x$.326 327When the optimum is unbounded, the optimal value computed for328the original problem will involve the big parameter.329In the original implementation of {\tt PipLib}, the big parameter could330even appear in some of the extra variables $\vec q$ created during331the application of a Gomory cut. The final result could then contain332implicit conditions on the big parameter through conditions on such333$\vec q$ variables. This problem was resolved in later versions334of {\tt PipLib} by taking $M$ to be divisible by any positive number.335The big parameter can then never appear in any $\vec q$ because336$\fract {\alpha M } = 0$. It should be noted, though, that an unbounded337problem 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 result342\parencite{Feautrier02}. Recent versions allow the user to simply343specify that the unknowns may be negative or that the maximum should344be 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 applications349of the big parameter and we believe that the user should not be bothered350with such implementation issues.351The current version of {\tt isl} therefore does not352provide any interface for specifying big parameters. Instead, the user353can specify whether a maximum needs to be computed and no assumptions354are made on the sign of the unknowns. Instead, the sign of the unknowns355is checked internally and a big parameter is automatically introduced when356needed. For compatibility with {\tt PipLib}, the {\tt isl\_pip} tool357does explicitly add non-negativity constraints on the unknowns unless358the \verb+Urs_unknowns+ option is specified.359Currently, there is also no way in {\tt isl} of expressing a big360parameter in the output. Even though361{\tt isl} makes the same divisibility assumption on the big parameter362as recent versions of {\tt PipLib}, it will therefore eventually363produce an error if the problem turns out to be unbounded.364 365\subsection{Preprocessing}366 367In this section, we describe some transformations that are368or can be applied in advance to reduce the running time369of 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 cuts374do 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 on377the original problem considered as a non-parametric problem378over the combined space of unknowns and parameters.379In fact, we do not simply check the feasibility, but we also380check for implicit equalities among the integer points by computing381the integer affine hull. The algorithm used is the same as that382described in \autoref{s:GBR} below.383Computing the affine hull is fairly expensive, but it can384bring huge benefits if any equalities can be found or if the problem385turns out to be empty.386 387\subsubsection{Constraint Simplification}388 389If the coefficients of the unknown and parameters in a constraint390have a common factor, then this factor should be removed, possibly391rounding down the constant term. For example, the constraint392$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 simplification395on 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$ inequalities402\parencite{Feautrier02}, but it is even better to \emph{exploit} the403equalities to reduce the dimensionality of the problem.404Given an equality involving at least one unknown, we pivot405the row corresponding to the equality with the column corresponding406to the last unknown with non-zero coefficient. The new column variable407can 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 initial410tableau remain lexicographically positive. In particular, if411the equality is of the form $b + \sum_{i \le j} a_i \, x_i = 0$ with412$a_j \ne 0$, then the (implicit) top rows of the initial tableau413are changed as follows414$$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 parameters444in 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, say452$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 parts455of the context and this context may have to be split for each456of the constraints. In the worst case, the basic algorithm may457have to consider all possible orderings of the constant terms.458Instead, {\tt isl} introduces a new parameter, say $u$, and459replaces the collection of constraints by the single460constraint $u + \sp {\vec a} {\vec x} \ge 0$ along with461context constraints $u \le b_i(\vec p)$.462Any solution to the new system is also a solution463to the original system since464$\sp {\vec a} {\vec x} \ge -u \ge -b_i(\vec p)$.465Conversely, $m = \min_i b_i(\vec p)$ satisfies the constraints466on $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 performed469in {\tt isl}.470The method described in this section can only detect symmetries471that are explicitly available in the input.472See \autoref{s:online} for the detection473and exploitation of symmetries that appear during the course of474the dual simplex method.475 476Note that the replacement of the $b_i(\vec p)$ by $u$ may lose477information if the parameters that occur in $b_i(\vec p)$ also478occur in other constraints. The replacement is therefore currently479only applied when all the parameters in all of the $b_i(\vec p)$480only occur in a single constraint, i.e., the one in which481the 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 experiments485of \autoref{s:pip:experiments} did not take into account486this single-occurrence constraint.487 488\subsubsection{Parameter Compression}\label{s:compression}489 490It may in some cases be apparent from the equalities in the problem491description that there can only be a solution for a sublattice492of the parameters. In such cases ``parameter compression''493\parencite{Meister2004PhD,Meister2008} can be used to replace494the parameters by alternative ``dense'' parameters.495For example, if there is a constraint $2x = n$, then the system496will only have solutions for even values of $n$ and $n$ can be replaced497by $2n'$. Similarly, the parameters $n$ and $m$ in a system with498the 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 to502preserve the lexicographical order. Moreover, due to our handling503of equalities described above there should be504no 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 context512based on a parametric constant term in the main tableau with indeterminate513sign. Each of these nodes may introduce extra variables in the context514corresponding to integer divisions. Each leaf of the tree prescribes515the solution in that part of the context that satisfies all the conditions516on the path leading to the leaf.517Such a quast is a very economical way of representing the solution, but518it would not be suitable as the (only) internal representation of519sets and relations in {\tt isl}. Instead, {\tt isl} represents520the constraints of a set or relation in disjunctive normal form.521The result of a parametric integer programming problem is then also522converted to this internal representation. Unfortunately, the conversion523to disjunctive normal form can lead to an explosion of the size524of the representation.525In some cases, this overhead would have to be paid anyway in subsequent526operations, but in other cases, especially for outside users that just527want to solve parametric integer programming problems, we would like528to avoid this overhead in future. That is, we are planning on introducing529quasts or a related representation as one of several possible internal530representations and on allowing the output of {\tt isl\_pip} to optionally531be printed as a quast.532 533Currently, {\tt isl} also does not have an internal representation534for expressions such as $\min_i b_i(\vec p)$ from the offline535symmetry detection of \autoref{s:offline}.536Assume that one of these expressions has $n$ bounds $b_i(\vec p)$.537If the expression538does not appear in the affine expression describing the solution,539but only in the constraints, and if moreover, the expression540only appears with a positive coefficient, i.e.,541$\min_i b_i(\vec p) \ge f_j(\vec p)$, then each of these constraints542can simply be reduplicated $n$ times, once for each of the bounds.543Otherwise, a conversion to disjunctive normal form544leads to $n$ cases, each described as $u = b_i(\vec p)$ with constraints545$b_i(\vec p) \le b_j(\vec p)$ for $j > i$546and547$b_i(\vec p) < b_j(\vec p)$ for $j < i$.548Note that even though this conversion leads to a size increase549by a factor of $n$, not detecting the symmetry could lead to550an increase by a factor of $n!$ if all possible orderings end up being551considered.552 553\subsection{Context Tableau}\label{s:context}554 555The main operation that a context tableau needs to provide is a test556on the sign of an affine expression over the elements of the context.557This sign can be determined by solving two integer linear feasibility558problems, one with a constraint added to the context that enforces559the expression to be non-negative and one where the expression is560negative. 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 Gomory563cuts 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, that566uses generalized basis reduction.567We start with some optimizations that are shared between the two568implementations 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 that573they will not only say whether a set is empty or not, but if the set574is 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 such576witnesses, we can avoid many feasibility tests during the determination577of the signs of affine expressions. In particular, if the expression578evaluates to a positive number on some of these points and to a negative579number on some others, then no feasibility test needs to be performed.580If all the evaluations are non-negative, we only need to check for the581possibility of a negative value and similarly in case of all582non-positive evaluations. Finally, in the rare case that all points583evaluate 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 result585of the first test.586 587When a new constraint is added to the context, the points that588violate the constraint are temporarily removed. They are reconsidered589when we backtrack over the addition of the constraint, as they will590satisfy the negation of the constraint. It is only when we backtrack591over 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 the594witnesses 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 context601needs 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 row603to split on first. {\tt PipLib} uses a heuristic based on the (absolute)604sizes of the coefficients. In particular, it takes the largest coefficient605of each row and then selects the row where this largest coefficient is smaller606than those of the other rows.607 608In {\tt isl}, we take that row for which non-negativity of its constant609term implies non-negativity of as many of the constant terms of the other610rows as possible. The intuition behind this heuristic is that on the611positive side, we will have fewer negative and indeterminate signs,612while on the negative side, we need to perform a pivot, which may613affect any number of rows meaning that the effect on the signs614is difficult to predict. This heuristic is of course much more615expensive 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 steps621of the dual simplex method applied to this new context will be the same622or at least very similar to those taken on the original context, i.e.,623before the constraint was added. In {\tt isl}, we therefore apply624the dual simplex method incrementally on the context and backtrack625to a previous state when a constraint is removed again.626An initial implementation that was never made public would also627keep the Gomory cuts, but the current implementation backtracks628to before the point where Gomory cuts are added before adding629an extra constraint to the context.630Keeping the Gomory cuts has the advantage that the sample value631is always an integer point and that this point may also satisfy632the new constraint. However, due to the technique of maintaining633witnesses explained above,634we would not perform a feasibility test in such cases and then635the previously added cuts may be redundant, possibly resulting636in an accumulation of a large number of cuts.637 638If the parameters may be negative, then the same big parameter trick639used in the main tableau is applied to the context. This big parameter640is 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 because644the big parameter of the main tableau also appears645in the context tableau.646 647Finally, it was reported by \textcite{Galea2009personal}, who648worked on a parametric integer programming implementation649in {\tt PPL} \parencite{PPL},650that it is beneficial to add cuts for \emph{all} rational coordinates651in 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 checking657is 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 which662the entire unit hypercube positioned at that point lies in the context.663This set is described by translates of the constraints of the context664and if (rationally) non-empty, any rational point665in 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 projected669out first. As soon as the algorithm is invoked, we then also670incrementally keep track of this recession cone. The reduced basis671found by one call of the algorithm is also reused as initial basis672for the next call.673 674Some problems lead to the675introduction of many integer divisions. Within a given context,676some of these integer divisions may be equal to each other, even677if the expressions are not identical, or they may be equal to some678affine combination of other variables.679To detect such cases, we compute the affine hull of the context680each time a new integer division is added. The algorithm used681for computing this affine hull is that of \textcite{Karr1976affine},682while the points used in this algorithm are obtained by performing683integer feasibility checks on that part of the context outside684the current approximation of the affine hull.685The list of witnesses is used to construct an initial approximation686of the hull, while any extra points found during the construction687of the hull is added to this list.688Any equality found in this way that expresses an integer division689as an \emph{integer} affine combination of other variables is690propagated to the main tableau, where it is used to eliminate that691integer 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} directory700of the {\tt isl} distribution.701Easier problems such as the702test cases distributed with {\tt Pip\-Lib} can be solved so quickly703that we would only be measuring overhead such as input/output and conversions704and 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 problem712originating 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 inspiration718for 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} and722were used to drive the first, Gomory cuts based, implementation723in {\tt isl}.724The third and final group of test cases are borrowed from725\textcite{Bygde2010licentiate} and inspired the offline symmetry detection726of \autoref{s:offline}. Without symmetry detection, the running times727are 11s and 5.9s.728All running times of {\tt barvinok} and {\tt isl} include a conversion729to disjunctive normal form. Without this conversion, the final two730cases can be solved in 0.07s and 0.21s.731The {\tt PipLib} implementation has some fixed limits and will732sometimes report the problem to be too complex (TC), while on some other733problems it will run out of memory (OOM).734The {\tt barvinok} implementation does not support problems735with a non-trivial lineality space (line) nor maximization problems (max).736The Gomory cuts based {\tt isl} implementation was terminated after 1000737minutes on the first problem. The gbr version introduces some738overhead 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\hline745\hline746% bart.pip747Phideo & TC & 793m & $>$999m & 2.7s & 372m \\748\hline749e1 & 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\hline761cnt/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 of772\textcite{Bygde2010licentiate} and an analysis of the results773by the approximate MPA method developed by \textcite{Bygde2010licentiate}774have revealed that these problems contain many more symmetries775than can be detected using the offline method of \autoref{s:offline}.776In this section, we present an online detection mechanism that has777not been implemented yet, but that has shown promising results778in manual applications.779 780Let us first consider what happens when we do not perform offline781symmetry detection. At some point, one of the782$b_i(\vec p) + \sp {\vec a} {\vec x} \ge 0$ constraints,783say the $j$th constraint, appears as a column784variable, say $c_1$, while the other constraints are represented785as rows of the form $b_i(\vec p) - b_j(\vec p) + c$.786The context is then split according to the relative order of787$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 minimum790of these $b_i(\vec p)$.791In the online method the split is similarly avoided by the introduction792of a new parameter. In particular, a new parameter is introduced793that represents794$\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 row798of the tableau such that the sign of $b(\vec p)$ is indeterminate799and 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$ with803context constraints $t \ge -b(\vec p)$ and $t \ge 0$ and replace804the column variable $c_j$ by $c' + t$. The row $r$ is now equal805to $b(\vec p) + t + c' - f$. The constant term of this row is always806non-negative because any negative value of $b(\vec p)$ is compensated807by $t \ge -b(\vec p)$ while and non-negative value remains non-negative808because $t \ge 0$.809 810We need to show that this transformation does not eliminate any valid811solutions and that it does not introduce any spurious solutions.812Given a valid solution for the original problem, we need to find813a non-negative value of $c'$ satisfying the constraints.814If $b(\vec p) \ge 0$, we can take $t = 0$ so that815$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 in820$t = \left|-b(\vec p)\right|_+ = \max(-b(\vec p), 0)$.821Conversely, given a solution to the new problem, we need to find822a 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 in826\autoref{s:post}, but, as in the case of offline symmetry detection,827it may be better to provide a direct representation for such828expressions in the internal representation of sets and relations829or 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 and841$k \in \Z_{\ge 1}$842a positive number, then power $k$ of relation $R$ is defined as843\begin{equation}844\label{eq:transitive:power}845R^k \coloneqq846\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 union858of all positive powers of $R$,859$$860R^+ \coloneqq \bigcup_{k \ge 1} R^k861.862$$863\end{definition}864Alternatively, the transitive closure may be defined865inductively as866\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 relation873may no longer be a polyhedral relation \parencite{Kelly1996closure},874we can, in the general case, only compute an approximation875of 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 approximation880to be as close as possible to the actual transitive closure881$R^+$ and we want to detect the cases where the approximation is882exact, 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 project887out the parameter $k$ from the resulting relation.888 889\begin{example}890As a trivial example, consider the relation891$R = \{\, x \to x + 1 \,\}$. The $k$th power of this map892for arbitrary $k$ is893$$894R^k = k \mapsto \{\, x \to x + k \mid k \ge 1 \,\}895.896$$897The transitive closure is then898$$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 equal915to $R$ itself.916 917In general, it is impossible to construct a closed form918of $R^k$ as a polyhedral relation.919We will therefore need to make some approximations.920As a first approximations, we will consider each of the basic921relations in $R$ as simply adding one or more offsets to a domain element922to arrive at an image element and ignore the fact that some of these923offsets 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 paths926that move through927a total of $k$ offsets and then intersect domain and range of this928collection 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}935with936\begin{equation}937\label{eq:transitive:path}938P = \vec s \mapsto \{\, \vec x \to \vec y \mid939\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_i941\wedge942\sum_i k_i = k > 0943\,\}944\end{equation}945and with $\Delta_i$ the basic sets that compose946the difference set $\Delta\,R$.947Note that the number of basic sets $\Delta_i$ need not be948the same as the number of basic relations in $R$.949Also note that since addition is commutative, it does not950matter in which order we add the offsets and so we are allowed951to group them as we did in \eqref{eq:transitive:path}.952 953If all the $\Delta_i$s are singleton sets954$\Delta_i = \{\, \vec \delta_i \,\}$ with $\vec \delta_i \in \Z^d$,955then \eqref{eq:transitive:path} simplifies to956\begin{equation}957\label{eq:transitive:singleton}958P = \{\, \vec x \to \vec y \mid959\exists k_i \in \Z_{\ge 0} :960\vec y = \vec x + \sum_i k_i \, \vec \delta_i961\wedge962\sum_i k_i = k > 0963\,\}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 if968some of $\vec \delta_i$s are parametric, then we need969to resort to further approximations.970 971To ease both the exposition and the implementation, we will for972the remainder of this section work with extended offsets973$\Delta_i' = \Delta_i \times \{\, 1 \,\}$.974That is, each offset is extended with an extra coordinate that is975set equal to one. The paths constructed by summing such extended976offsets have the length encoded as the difference of their977final coordinates. The path $P'$ can then be decomposed into978paths $P_i'$, one for each $\Delta_i$,979\begin{equation}980\label{eq:transitive:decompose}981P' = \left(982(P_m' \cup \identity) \circ \cdots \circ983(P_2' \cup \identity) \circ984(P_1' \cup \identity)985\right) \cap986\{\,987\vec x' \to \vec y' \mid y_{d+1} - x_{d+1} = k > 0988\,\}989,990\end{equation}991with992$$993P_i' = \vec s \mapsto \{\, \vec x' \to \vec y' \mid994\exists k \in \Z_{\ge 1}, \vec \delta \in k \, \Delta_i'(\vec s) :995\vec y' = \vec x' + \vec \delta996\,\}997.998$$999Note that each $P_i'$ contains paths of length at least one.1000We therefore need to take the union with the identity relation1001when composing the $P_i'$s to allow for paths that do not contain1002any offsets from one or more $\Delta_i'$.1003The path that consists of only identity relations is removed1004by imposing the constraint $y_{d+1} - x_{d+1} > 0$.1005Taking the union with the identity relation means that1006that 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 application1009of the composition operation may therefore result in a relation1010with $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 together1016and handled as in \eqref{eq:transitive:singleton}.1017Note that this is just an optimization. The procedure described1018below 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 end1022of this section.1023Without existentially quantified variables, we can classify1024the constraints of $\Delta_i'$ as follows1025\begin{enumerate}1026\item non-parametric constraints1027\begin{equation}1028\label{eq:transitive:non-parametric}1029A_1 \vec x + \vec c_1 \geq \vec 01030\end{equation}1031\item purely parametric constraints1032\begin{equation}1033\label{eq:transitive:parametric}1034B_2 \vec s + \vec c_2 \geq \vec 01035\end{equation}1036\item negative mixed constraints1037\begin{equation}1038\label{eq:transitive:mixed}1039A_3 \vec x + B_3 \vec s + \vec c_3 \geq \vec 01040\end{equation}1041such that for each row $j$ and for all $\vec s$,1042$$1043\Delta_i'(\vec s) \cap1044\{\, \vec \delta' \mid B_{3,j} \vec s + c_{3,j} > 0 \,\}1045= \emptyset1046$$1047\item positive mixed constraints1048$$1049A_4 \vec x + B_4 \vec s + \vec c_4 \geq \vec 01050$$1051such that for each row $j$, there is at least one $\vec s$ such that1052$$1053\Delta_i'(\vec s) \cap1054\{\, \vec \delta' \mid B_{4,j} \vec s + c_{4,j} > 0 \,\}1055\ne \emptyset1056$$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 \mapsto1063\{\,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 01071\wedge1072B_2 \vec s + \vec c_2 \geq \vec 01073\wedge1074A_3 \vec f + B_3 \vec s + \vec c_3 \geq \vec 01075\,\}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 every1081$k \in \Z_{\ge 1}$ and for every $\vec f \in k \, \Delta_i(\vec s)$1082we have that1083$(\vec f, k)$ satisfies the constraints in \eqref{eq:transitive:Q}.1084If $\Delta_i(\vec s)$ is non-empty, then $\vec s$ must satisfy1085the constraints in \eqref{eq:transitive:parametric}.1086Each element $(\vec f, k) \in k \, \Delta_i'(\vec s)$ is a sum1087of $k$ elements $(\vec f_j, 1)$ in $\Delta_i'(\vec s)$.1088Each of these elements satisfies the constraints in1089\eqref{eq:transitive:non-parametric}, i.e.,1090$$1091\left[1092\begin{matrix}1093A_1 & \vec c_11094\end{matrix}1095\right]1096\left[1097\begin{matrix}1098\vec f_j \\ 11099\end{matrix}1100\right]1101\ge \vec 01102.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 such1107that 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 the1112rational 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 cone1117$\{\, \vec x' \in \Q^{d+1} \mid \left[1118\begin{matrix}1119A_1 & \vec c_11120\end{matrix}1121\right] \vec x' \,\}$ and therefore $\Delta'_i(\vec s)$ is1122a 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 may1126not compute the most accurate affine approximation of1127$k \, \Delta_i(\vec s)$ with $k \ge 1$.1128In particular, we only consider the negative mixed constraints that1129happen to appear in the description of $\Delta_i(\vec s)$, while we1130should instead consider \emph{all} valid such constraints.1131It is also sufficient to consider those constraints because any1132constraint that is valid for $k \, \Delta_i(\vec s)$ is also1133valid for $1 \, \Delta_i(\vec s) = \Delta_i(\vec s)$.1134Take therefore any constraint1135$\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)$ iff1137$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 constraint1141is not valid for $k \, \Delta_i(\vec s)$.1142We therefore need to impose $\spv b s + c \le 0$ for all values1143of $\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 valid1145constraint of $\Delta_i(\vec s)$. That is, $(\vec b, c)$ are the opposites1146of the coefficients of a valid constraint of $\Delta_i(\vec s)$.1147The approximation of $k \, \Delta_i(\vec s)$ can therefore be obtained1148using three applications of Farkas' lemma. The first obtains the coefficients1149of constraints valid for $\Delta_i(\vec s)$. The second obtains1150the coefficients of constraints valid for the projection of $\Delta_i(\vec s)$1151onto the parameters. The opposite of the second set is then computed1152and intersected with the first set. The result is the set of coefficients1153of constraints valid for $k \, \Delta_i(\vec s)$. A final application1154of Farkas' lemma is needed to obtain the approximation of1155$k \, \Delta_i(\vec s)$ itself.1156 1157\begin{example}1158Consider the relation1159$$1160n \to \{\, (x, y) \to (1 + x, 1 - n + y) \mid n \ge 2 \,\}1161.1162$$1163The difference set of this relation is1164$$1165\Delta = n \to \{\, (1, 1 - n) \mid n \ge 2 \,\}1166.1167$$1168Using our approach, we would only consider the mixed constraint1169$y - 1 + n \ge 0$, leading to the following approximation of the1170transitive 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 obtain1182\begin{verbatim}1183{ rat: coefficients[[c_cst, c_n] -> [i2, i3]] : i3 <= c_n and1184 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}1192are1193\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 obtain1204\begin{verbatim}1205{ rat: [[c_cst, c_n] -> [i2, i3]] : i3 <= c_n and1206 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 then1214\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}1224is1225\begin{verbatim}1226[n] -> { [x, y] -> [o0, o1] : o1 <= x + y - o0 and1227 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 by1232classifying them into variables that are uniquely1233determined by the parameters, variables that are independent1234of the parameters and others. The first set can be treated1235as parameters and the second as variables. Constraints involving1236the other existentially quantified variables are removed.1237 1238\begin{example}1239Consider the relation1240$$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 is1246$$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 terms1252of the parameters and variables as1253$$1254\alpha_0 = \floor{\frac{-2 + n}7}1255\qquad1256\text{and}1257\qquad1258\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 as1264a purely parametric constraint, while the other two constraints are1265non-parametric.1266The corresponding $Q$~\eqref{eq:transitive:Q} is therefore1267$$1268\begin{aligned}1269n \to \{\, (x,z) \to (y,w) \mid1270\exists\, \alpha_0, \alpha_1, k, f : {} &1271k \ge 1 \wedge1272y = x + f \wedge1273w = z + k \wedge {} \\1274&12757\alpha_0 = -2 + n \wedge12765\alpha_1 = -k + x \wedge1277x \ge 6 k1278\,\}1279.1280\end{aligned}1281$$1282Projecting out the final coordinates encoding the length of the paths,1283results in the exact transitive closure1284$$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 leads1292to a loss of accuracy. In some cases, some of this loss can be recovered1293by not considering the parameters in a special way.1294That is, instead of considering the set1295$$1296\Delta = \diff R =1297\vec s \mapsto1298\{\, \vec \delta \in \Z^{d} \mid \exists \vec x \to \vec y \in R :1299\vec \delta = \vec y - \vec x1300\,\}1301$$1302we consider the set1303$$1304\Delta' = \diff R' =1305\{\, \vec \delta \in \Z^{n+d} \mid \exists1306(\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 equivalent1313to projecting out the parameters in $\Delta$.1314The result is obviously a superset of $\Delta$, but all its constraints1315are of type \eqref{eq:transitive:non-parametric} and they can therefore1316all be used in the construction of $Q_i$.1317 1318\begin{example}1319Consider the relation1320$$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 have1326$$1327\diff R = n \to \{\, (1, 1 - n) \mid n \ge 2 \,\}1328$$1329and so, by treating the parameters in a special way, we obtain1330the 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 instead1336$$1337R' = \{\, (n, x, y) \to (n, 1 + x, 1 - n + y) \mid n \ge 2 \,\}1338$$1339then1340$$1341\diff R' = \{\, (0, 1, y) \mid y \le -1 \,\}1342$$1343and we obtain the approximation1344$$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 obtain1349$$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 that1354can be obtained. That would be1355$$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 obtained1364by 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 an1367overapproximation of $R^+$.1368To check whether the results are exact, we need to consider two1369cases depending on whether $R$ is {\em cyclic}, where $R$ is defined1370to 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 of1373\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 check1380whether1381$$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 lexicographically1387forward relations, a special case of acyclic relations.1388 1389If, on the other hand, $R$ is cyclic, then we have to resort1390to checking whether the approximation $K$ of the power is exact.1391Note that $T$ may be exact even if $K$ is not exact, so the check1392is sound, but incomplete.1393To check exactness of the power, we simply need to check1394\eqref{eq:transitive:power}. Since again $K$ is known1395to be an overapproximation, we only need to check whether1396$$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 R1405\wedge y_{d+1} - x_{d+1} = 1\,\}$, i.e., $R$ extended with path1406lengths 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, even1410in the acyclic case, so we only need to be careful that we find1411all cyclic cases. Now, if $R$ is cyclic, i.e.,1412$R^+ \cap \identity \ne \emptyset$, then, since $T$ is1413an overapproximation of $R^+$, also1414$T \cap \identity \ne \emptyset$. This in turn means1415that $\Delta \, K'$ contains a point whose first $d$ coordinates1416are 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 approximation1419is 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 relations1425that can be partially ordered1426then the accuracy of the approximation may be improved by computing1427an 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 later1430pass 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 closures1437of $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 segments1442in any path that moves through both $R_1$ and $R_2$ to1443first move through $R_1$ and then through $R_2$.1444 1445This idea can be generalized to relations that are unions1446of more than two basic relations by constructing the1447strongly connected components in the graph with as vertices1448the basic relations and an edge between two basic relations1449$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$ iff1451\begin{equation}1452\label{eq:transitive:edge}1453R_i \circ R_j1454\not\subseteq1455R_j \circ R_i1456.1457\end{equation}1458The components can be obtained from the graph by applying1459Tarjan's algorithm \parencite{Tarjan1972}.1460 1461In practice, we compute the (extended) powers $K_i'$ of each component1462separately and then compose them as in \eqref{eq:transitive:decompose}.1463Note, however, that in this case the order in which we apply them is1464important and should correspond to a topological ordering of the1465strongly connected components. Simply applying Tarjan's1466algorithm 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 between1469two 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 check1473is 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 will1478still be an overapproximation of the left hand side, but this result1479may not be transitively closed. If we only separate components based1480on the condition $R_i \circ R_j = \emptyset$, then there is no problem,1481as this condition will still hold on the computed approximations1482of the transitive closures. If, however, we have exploited1483\eqref{eq:transitive:edge} during the decomposition and if the1484result turns out not to be exact, then we check whether1485the result is transitively closed. If not, we recompute1486the transitive closure, skipping the decomposition.1487Note that testing for transitive closedness on the result may1488be fairly expensive, so we may want to make this check1489configurable.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 with1512the Omega calculator~\parencite{Omega_calc}, $R = R_1 \cup R_2$,1513with1514$$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 have1524$$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 so1535$$1536\left(1537R_1 \cup R_21538\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}>01563 \ifnum\value{t3}>0\ifnum\value{t4}>01564 \draw[thick,->] (\i,\j) to[out=20] (\i+3,\j);1565 \fi\fi\fi\fi1566 \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}>01571 \ifnum\value{t3}>0\ifnum\value{t4}>01572 \draw[thick,->] (\i,\j) to[in=-20,out=20] (\i,\j+3);1573 \fi\fi\fi\fi1574 \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}>01579 \ifnum\value{t3}>0\ifnum\value{t4}>01580 \draw[thick,->] (\i,\j) to (\i+1,\j+1);1581 \fi\fi\fi\fi1582 }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$,1594with1595$$1596\begin{aligned}1597R_1 &= n \mapsto \{\, (i,j) \to (i+3,j) \mid1598i \le 2 j - 4 \wedge1599i \le n - 3 \wedge1600j \le 2 i - 1 \wedge1601j \le n \,\}1602\\1603R_2 &= n \mapsto \{\, (i,j) \to (i,j+3) \mid1604i \le 2 j - 1 \wedge1605i \le n \wedge1606j \le 2 i - 4 \wedge1607j \le n - 3 \,\}1608\\1609R_3 &= n \mapsto \{\, (i,j) \to (i+1,j+1) \mid1610i \le 2 j - 1 \wedge1611i \le n - 1 \wedge1612j \le 2 i - 1 \wedge1613j \le n - 1\,\}1614.1615\end{aligned}1616$$1617The figure shows this relation for $n = 7$.1618Both1619$R_3 \circ R_1 \subseteq R_1 \circ R_3$1620and1621$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 and1625 j <= 2 i - 1 and j <= n };1626R2 := [n] -> { [i,j] -> [i,j+3] : i <= 2 j - 1 and i <= n and1627 j <= 2 i - 4 and j <= n - 3 };1628R3 := [n] -> { [i,j] -> [i+1,j+1] : i <= 2 j - 1 and i <= n - 1 and1629 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 both1635$R_2 \circ R_1 \not\subseteq R_1 \circ R_2$1636and1637$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$ separately1640and composing the results, the power of $R$ can be computed exactly1641using \eqref{eq:transitive:singleton}.1642As explained by \textcite{Beletska2009}, applying the same formula1643to $R$ directly, without a decomposition, would result in1644an 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 given1652abstract domain to the same domain.1653However, if $R$ is a union of relations that map between different1654domains, then this assumption no longer holds.1655In particular, when an entire dependence graph is encoded1656in a single relation, as is done by, e.g.,1657\textcite[Section~6.1]{Barthou2000MSE}, then it does not make1658sense to look at differences between iterations of different domains.1659Now, arguably, a modified Floyd-Warshall algorithm should1660be applied to the dependence graph, as advocated by1661\textcite{Kelly1996closure}, with the transitive closure operation1662only being applied to relations from a given domain to itself.1663However, it is also possible to detect disjoint domains and ranges1664and to apply Floyd-Warshall internally.1665 1666\LinesNumbered1667\begin{algorithm}1668\caption{The modified Floyd-Warshall algorithm of1669\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 relation1675$R_{pq}$ contains all indirect paths from $p$ to $q$ in the input graph}1676%1677\BlankLine1678\SetAlgoVlined1679\DontPrintSemicolon1680%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 is1698obtained. 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 vertices1701the parts of the partition and as edges the $R_i$ attached to1702the appropriate pairs of vertices.1703In particular, let there be $n$ parts $P_k$ in the partition.1704We construct $n^2$ relations1705$$1706R_{pq} \coloneqq \bigcup_{i \text{ s.t. } \domain R_i \subseteq P_p \wedge1707 \range R_i \subseteq P_q} R_i1708,1709$$1710apply \autoref{a:Floyd} and return the union of all resulting1711$R_{pq}$ as the transitive closure of $R$.1712Each iteration of the $r$-loop in \autoref{a:Floyd} updates1713all 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 that1716pass through earlier vertices since $R_{rr}$ itself has been updated1717accordingly 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 the1720$r$-loop in Line~\ref{l:Floyd:update}.1721However, from an implementation perspective, it is easier1722to allow either or both of these to have been updated1723in the same iteration of the $r$-loop.1724This may result in duplicate paths, but these can usually1725be removed by coalescing (\autoref{s:coalescing}) the result of the union1726in 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 call1729includes the partitioning step, but the resulting partition will1730usually be a singleton.1731The result of the recursive call will either be exact or an1732overapproximation. The final result of Floyd-Warshall is therefore1733also 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 of1753\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 of1758\textcite{Beletska2009},1759reproduced in \autoref{f:COCOA:1}.1760This relation can be described as1761$$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 range1769of the rightward relation and vice versa, but that the domain1770of neither relation overlaps with its own range or the domain of1771the 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 have1776$$1777\begin{aligned}1778R_{00} & = \emptyset1779\\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} & = \emptyset1788.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, but1793$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 not1796changed in the second iteration and it does not have an effect1797on $R_{01}$ and $R_{10}$. However, $R_{00}$ is updated to1798include $R_{10} \circ R_{01}$, i.e., the dotted arrow in the figure.1799The transitive closure of the original relation is then equal to1800$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 closure1807of union of basic relations incrementally. In particular,1808if $R$ is a union of $m$ basic maps,1809$$1810R = \bigcup_j R_j1811,1812$$1813then we can pick some $R_i$ and compute the transitive closure of $R$ as1814\begin{equation}1815\label{eq:transitive:incremental}1816R^+ = R_i^+ \cup1817\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 each1824of the disjuncts in the argument of the second transitive1825closure in \eqref{eq:transitive:incremental} be representable1826as 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 closure1829can be reduced by one.1830Now, $R_i^* = R_i^+ \cup \identity$, but in some cases it is possible1831to relax the constraints of $R_i^+$ to include part of the identity relation,1832say on domain $D$. We will use the notation1833${\cal C}(R_i,D) = R_i^+ \cup \identity_D$ to represent1834this 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 attain1837the value $0$ in \eqref{eq:transitive:Q} and by using1838$$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 domain1844and 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 for1848the 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 chosen1852to include both $\domain R$ and $\range R$, i.e., such1853that $\identity_D \circ R_j \circ \identity_D = R_j$ for all $j\ne i$.1854\textcite{Kelly1996closure} say that they use1855$D = \domain R_i \cup \range R_i$, but presumably they mean that1856they 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 to1866compose with ${\cal C}(R_i,D)$ more selectively.1867In particular, if we have1868\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 = \emptyset1872\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 = \emptyset1878\end{equation}1879then we can refine \eqref{eq:transitive:incremental} to1880$$1881R_i^+ \cup1882\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\cup1889\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_j1894\right)1895\cup1896\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\cup1903\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_j1908\right)1909\right)^+1910.1911$$1912If only property~\eqref{eq:transitive:right} holds,1913we can use1914$$1915R_i^+ \cup1916\left(1917\left(1918R_i^+ \cup \identity1919\right)1920\circ1921\left(1922\left(1923\bigcup_{\shortstack{$\scriptstyle\domain R_j \subseteq D $}}1924R_j \circ {\cal C}1925\right)1926\cup1927\left(1928\bigcup_{\shortstack{$\scriptstyle\domain R_j \cap \range R_i = \emptyset$}}1929\!\!\!\!\!1930R_j1931\right)1932\right)^+1933\right)1934,1935$$1936while if only property~\eqref{eq:transitive:left} holds,1937we can use1938$$1939R_i^+ \cup1940\left(1941\left(1942\left(1943\bigcup_{\shortstack{$\scriptstyle\range R_j \subseteq D $}}1944{\cal C} \circ R_j1945\right)1946\cup1947\left(1948\bigcup_{\shortstack{$\scriptstyle\range R_j \cap \domain R_i = \emptyset$}}1949\!\!\!\!\!1950R_j1951\right)1952\right)^+1953\circ1954\left(1955R_i^+ \cup \identity1956\right)1957\right)1958.1959$$1960 1961It should be noted that if we want the result of the incremental1962approach to be transitively closed, then we can only apply it1963if 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 contain1966the 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} is1971designed 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 algorithm1974that is based on their ideas.1975Note that the {\tt Omega} library computes underapproximations1976\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'' relation1980$$1981\{\, \vec i \to \vec j \mid \exists \vec \alpha :1982\vec L \le \vec j - \vec i \le \vec U1983\wedge1984(\forall p : j_p - i_p = M_p \alpha_p)1985\,\}1986,1987$$1988where $p$ ranges over the dimensions and $\vec L$, $\vec U$ and1989$\vec M$ are constant integer vectors. The elements of $\vec U$1990may be $\infty$, meaning that there is no upper bound corresponding1991to 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 is1995\begin{equation}1996\label{eq:omega}1997\{\, \vec i \to \vec j \mid \exists \vec \alpha, k :1998k \ge 1 \wedge1999k \, \vec L \le \vec j - \vec i \le k \, \vec U2000\wedge2001(\forall p : j_p - i_p = M_p \alpha_p)2002\,\}2003.2004\end{equation}2005The domain and range of this transitive closure are then2006intersected 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 authors2010use the above algorithm as a substep on the disjuncts in the relation.2011At the end, they say2012\begin{quote}2013If an upper bound is required, it can be calculated in a manner2014similar to that of a single conjunct [sic] relation.2015\end{quote}2016Presumably, the authors mean that a ``d-form'' approximation2017of the whole input relation should be used.2018However, the accuracy can be improved by also trying to2019apply 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 by2022allowing the value zero for $k$ in \eqref{eq:omega},2023i.e., by computing2024$$2025\{\, \vec i \to \vec j \mid \exists \vec \alpha, k :2026k \ge 0 \wedge2027k \, \vec L \le \vec j - \vec i \le k \, \vec U2028\wedge2029(\forall p : j_p - i_p = M_p \alpha_p)2030\,\}2031.2032$$2033In our implementation we take as $D$ the simple hull2034(\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 by2037\textcite{Kelly1996closure}:2038${\cal C}(R_i,D) - R_i^+$ is not a union and for each $j \ne i$2039the condition2040$$2041\left({\cal C}(R_i,D) - R_i^+\right)2042\circ2043R_j2044\circ2045\left({\cal C}(R_i,D) - R_i^+\right)2046=2047R_j2048$$2049holds.2050