Català-Valencià – Catalan中文 – Chinese (Simplified)中文 – Chinese (Traditional)Česky – CzechDansk – DanishNederlands – DutchEnglish – EnglishSuomi – FinnishFrançais – FrenchDeutsch – Germanעברית – Hebrewहिंदी – HindiMagyar – HungarianBahasa Indonesia – IndonesianItaliano – Italian日本語 – Japanese한국어 – KoreanМакедонски – Macedonianमराठी – MarathiNorsk – NorwegianPolski – PolishPortuguês – PortuguesePortuguês – Portuguese (Brazil)Русский – RussianSlovenčina – SlovakSlovenščina – SlovenianEspañol – SpanishSvenska – SwedishTürkçe – TurkishУкраїнська – UkrainianOëzbekcha – Uzbek

nexmon – Blame information for rev 1

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