transitive closure: project out parameters when any constraints are impure

The resulting delta set can lead to extra constraints on the path,
resulting in a possibly more accurate approximation.

Signed-off-by: Sven Verdoolaege <skimo@kotnet.org>

3 files changed, 117 insertions, 5 deletions

diff --git a/doc/implementation.tex b/doc/implementation.tex
index 758d27ee..d79f6550 100644
--- a/doc/implementation.tex
+++ b/doc/implementation.tex
@@ -94,7 +94,7 @@ The difference set ($\Delta \, R$) of $R$ is the set
 of differences between image elements and the corresponding
 domain elements,
 $$
-\Delta \, R \coloneqq
+\diff R \coloneqq
 \vec s \mapsto
 \{\, \vec \delta \in \Z^{d} \mid \exists \vec x \to \vec y \in R :
 \vec \delta = \vec y - \vec x
@@ -505,6 +505,76 @@ n \to \{\, x \to y \mid \exists \, \alpha_0, \alpha_1: 7\alpha_1 = -2 + n \wedge
 $$
 \end{example}
 
+The fact that we ignore some impure constraints clearly leads
+to a loss of accuracy.  In some cases, some of this loss can be recovered
+by not considering the parameters in a special way.
+That is, instead of considering the set
+$$
+\Delta = \diff R =
+\vec s \mapsto
+\{\, \vec \delta \in \Z^{d} \mid \exists \vec x \to \vec y \in R :
+\vec \delta = \vec y - \vec x
+\,\}
+$$
+we consider the set
+$$
+\Delta' = \diff R' =
+\{\, \vec \delta \in \Z^{n+d} \mid \exists
+(\vec s, \vec x) \to (\vec s, \vec y) \in R' :
+\vec \delta = (\vec s - \vec s, \vec y - \vec x)
+\,\}
+.
+$$
+The first $n$ coordinates of every element in $\Delta'$ are zero.
+Projecting out these zero coordinates from $\Delta'$ is equivalent
+to projecting out the parameters in $\Delta$.
+The result is obviously a superset of $\Delta$, but all its constraints
+are of type \eqref{eq:transitive:non-parametric} and they can therefore
+all be used in the construction of $Q_i$.
+
+\begin{example}
+Consider the relation
+$$
+% [n] -> { [x, y] -> [1 + x, 1 - n + y] | n >= 2 }
+R = n \to \{\, (x, y) \to (1 + x, 1 - n + y) \mid n \ge 2 \,\}
+.
+$$
+We have
+$$
+\diff R = n \to \{\, (1, 1 - n) \mid n \ge 2 \,\}
+$$
+and so, by treating the parameters in a special way, we obtain
+the following approximation for $R^+$:
+$$
+n \to \{\, (x, y) \to (x', y') \mid n \ge 2 \wedge y' \le 1 - n + y \wedge x' \ge 1 + x \,\}
+.
+$$
+If we consider instead
+$$
+R' = \{\, (n, x, y) \to (n, 1 + x, 1 - n + y) \mid n \ge 2 \,\}
+$$
+then
+$$
+\diff R' = \{\, (0, 1, y) \mid y \le -1 \,\}
+$$
+and we obtain the approximation
+$$
+n \to \{\, (x, y) \to (x', y') \mid n \ge 2 \wedge x' \ge 1 + x \wedge y' \le x + y - x' \,\}
+.
+$$
+If we consider both $\diff R$ and $\diff R'$, then we obtain
+$$
+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' \,\}
+.
+$$
+Note, however, that this is not the most accurate affine approximation that
+can be obtained.  That would be
+$$
+n \to \{\, (x, y) \to (x', y') \mid y' \le 2 - n + x + y - x' \wedge n \ge 2 \wedge x' \ge 1 + x \,\}
+.
+$$
+\end{example}
+
 \subsection{Checking Exactness}
 
 The approximation $T$ for the transitive closure $R^+$ can be obtained
diff --git a/doc/manual.tex b/doc/manual.tex
index 191f4d19..eae014e8 100644
--- a/doc/manual.tex
+++ b/doc/manual.tex
@@ -38,6 +38,7 @@
 \def\domain{\mathop{\rm dom}\nolimits}
 \def\range{\mathop{\rm ran}\nolimits}
 \def\identity{\mathop{\rm Id}\nolimits}
+\def\diff{\mathop{\Delta}\nolimits}
 
 \providecommand{\floor}[1]{\left\lfloor#1\right\rfloor}
 
diff --git a/isl_transitive_closure.c b/isl_transitive_closure.c
index 1452fa99..bd0bd7a3 100644
--- a/isl_transitive_closure.c
+++ b/isl_transitive_closure.c
@@ -322,6 +322,9 @@ error:
  * 	variables are non-zero and if moreover the parametric constant
  * 	can never attain positive values.
  * Return IMPURE otherwise.
+ *
+ * If div_purity is NULL then we are dealing with a non-parametric set
+ * and so the constraint is obviously PURE_VAR.
  */
 static int purity(__isl_keep isl_basic_set *bset, isl_int *c, int *div_purity,
 	int eq)
@@ -333,6 +336,9 @@ static int purity(__isl_keep isl_basic_set *bset, isl_int *c, int *div_purity,
 	int i;
 	int p = 0, v = 0;
 
+	if (!div_purity)
+		return PURE_VAR;
+
 	n_div = isl_basic_set_dim(bset, isl_dim_div);
 	d = isl_basic_set_dim(bset, isl_dim_set);
 	nparam = isl_basic_set_dim(bset, isl_dim_param);
@@ -439,10 +445,13 @@ error:
 	return -1;
 }
 
+/* If any of the constraints is found to be impure then this function
+ * sets *impurity to 1.
+ */
 static __isl_give isl_basic_map *add_delta_constraints(
 	__isl_take isl_basic_map *path,
 	__isl_keep isl_basic_set *delta, unsigned off, unsigned nparam,
-	unsigned d, int *div_purity, int eq)
+	unsigned d, int *div_purity, int eq, int *impurity)
 {
 	int i, k;
 	int n = eq ? delta->n_eq : delta->n_ineq;
@@ -456,6 +465,8 @@ static __isl_give isl_basic_map *add_delta_constraints(
 		int p = purity(delta, delta_c[i], div_purity, eq);
 		if (p < 0)
 			goto error;
+		if (p != PURE_VAR && p != PURE_PARAM && !*impurity)
+			*impurity = 1;
 		if (p == IMPURE)
 			continue;
 		if (eq && p != MIXED) {
@@ -505,7 +516,7 @@ error:
  *
  * In particular, let delta be defined as
  *
- *	\delta = [p] -> { [x] : A x + a >= and B p + b >= 0 and
+ *	\delta = [p] -> { [x] : A x + a >= 0 and B p + b >= 0 and
  *				C x + C'p + c >= 0 and
  *				D x + D'p + d >= 0 }
  *
@@ -532,6 +543,17 @@ error:
  * parameter dependent and others.  Constraints containing
  * any of the other existentially quantified variables are removed.
  * This is safe, but leads to an additional overapproximation.
+ *
+ * If there are any impure constraints, then we also eliminate
+ * the parameters from \delta, resulting in a set
+ *
+ *	\delta' = { [x] : E x + e >= 0 }
+ *
+ * and add the constraints
+ *
+ *			E f + k e >= 0
+ *
+ * to the constructed relation.
  */
 static __isl_give isl_map *path_along_delta(__isl_take isl_dim *dim,
 	__isl_take isl_basic_set *delta)
@@ -544,6 +566,7 @@ static __isl_give isl_map *path_along_delta(__isl_take isl_dim *dim,
 	int i, k;
 	int is_id;
 	int *div_purity = NULL;
+	int impurity = 0;
 
 	if (!delta)
 		goto error;
@@ -575,8 +598,26 @@ static __isl_give isl_map *path_along_delta(__isl_take isl_dim *dim,
 	if (!div_purity)
 		goto error;
 
-	path = add_delta_constraints(path, delta, off, nparam, d, div_purity, 1);
-	path = add_delta_constraints(path, delta, off, nparam, d, div_purity, 0);
+	path = add_delta_constraints(path, delta, off, nparam, d,
+				     div_purity, 1, &impurity);
+	path = add_delta_constraints(path, delta, off, nparam, d,
+				     div_purity, 0, &impurity);
+	if (impurity) {
+		isl_dim *dim = isl_basic_set_get_dim(delta);
+		delta = isl_basic_set_project_out(delta,
+						  isl_dim_param, 0, nparam);
+		delta = isl_basic_set_add(delta, isl_dim_param, nparam);
+		delta = isl_basic_set_reset_dim(delta, dim);
+		if (!delta)
+			goto error;
+		path = isl_basic_map_extend_constraints(path, delta->n_eq,
+							delta->n_ineq + 1);
+		path = add_delta_constraints(path, delta, off, nparam, d,
+					     NULL, 1, &impurity);
+		path = add_delta_constraints(path, delta, off, nparam, d,
+					     NULL, 0, &impurity);
+		path = isl_basic_map_gauss(path, NULL);
+	}
 
 	is_id = empty_path_is_identity(path, off + d);
 	if (is_id < 0)