isl_map_transitive_closure: break early if input map doesn't compose with itself
[platform/upstream/isl.git] / isl_transitive_closure.c
index 295c2da..ba7153b 100644 (file)
 #include "isl_map_private.h"
 #include "isl_seq.h"
  
+/* Given a map that represents a path with the length of the path
+ * encoded as the difference between the last output coordindate
+ * and the last input coordinate, set this length to either
+ * exactly "length" (if "exactly" is set) or at least "length"
+ * (if "exactly" is not set).
+ */
+static __isl_give isl_map *set_path_length(__isl_take isl_map *map,
+       int exactly, int length)
+{
+       struct isl_dim *dim;
+       struct isl_basic_map *bmap;
+       unsigned d;
+       unsigned nparam;
+       int k;
+       isl_int *c;
+
+       if (!map)
+               return NULL;
+
+       dim = isl_map_get_dim(map);
+       d = isl_dim_size(dim, isl_dim_in);
+       nparam = isl_dim_size(dim, isl_dim_param);
+       bmap = isl_basic_map_alloc_dim(dim, 0, 1, 1);
+       if (exactly) {
+               k = isl_basic_map_alloc_equality(bmap);
+               c = bmap->eq[k];
+       } else {
+               k = isl_basic_map_alloc_inequality(bmap);
+               c = bmap->ineq[k];
+       }
+       if (k < 0)
+               goto error;
+       isl_seq_clr(c, 1 + isl_basic_map_total_dim(bmap));
+       isl_int_set_si(c[0], -length);
+       isl_int_set_si(c[1 + nparam + d - 1], -1);
+       isl_int_set_si(c[1 + nparam + d + d - 1], 1);
+
+       bmap = isl_basic_map_finalize(bmap);
+       map = isl_map_intersect(map, isl_map_from_basic_map(bmap));
+
+       return map;
+error:
+       isl_basic_map_free(bmap);
+       isl_map_free(map);
+       return NULL;
+}
+
+/* Check whether the overapproximation of the power of "map" is exactly
+ * the power of "map".  Let R be "map" and A_k the overapproximation.
+ * The approximation is exact if
+ *
+ *     A_1 = R
+ *     A_k = A_{k-1} \circ R                   k >= 2
+ *
+ * Since A_k is known to be an overapproximation, we only need to check
+ *
+ *     A_1 \subset R
+ *     A_k \subset A_{k-1} \circ R             k >= 2
+ *
+ * In practice, "app" has an extra input and output coordinate
+ * to encode the length of the path.  So, we first need to add
+ * this coordinate to "map" and set the length of the path to
+ * one.
+ */
+static int check_power_exactness(__isl_take isl_map *map,
+       __isl_take isl_map *app)
+{
+       int exact;
+       isl_map *app_1;
+       isl_map *app_2;
+
+       map = isl_map_add(map, isl_dim_in, 1);
+       map = isl_map_add(map, isl_dim_out, 1);
+       map = set_path_length(map, 1, 1);
+
+       app_1 = set_path_length(isl_map_copy(app), 1, 1);
+
+       exact = isl_map_is_subset(app_1, map);
+       isl_map_free(app_1);
+
+       if (!exact || exact < 0) {
+               isl_map_free(app);
+               isl_map_free(map);
+               return exact;
+       }
+
+       app_1 = set_path_length(isl_map_copy(app), 0, 1);
+       app_2 = set_path_length(app, 0, 2);
+       app_1 = isl_map_apply_range(map, app_1);
+
+       exact = isl_map_is_subset(app_2, app_1);
+
+       isl_map_free(app_1);
+       isl_map_free(app_2);
+
+       return exact;
+}
+
+/* Check whether the overapproximation of the power of "map" is exactly
+ * the power of "map", possibly after projecting out the power (if "project"
+ * is set).
+ *
+ * If "project" is set and if "steps" can only result in acyclic paths,
+ * then we check
+ *
+ *     A = R \cup (A \circ R)
+ *
+ * where A is the overapproximation with the power projected out, i.e.,
+ * an overapproximation of the transitive closure.
+ * More specifically, since A is known to be an overapproximation, we check
+ *
+ *     A \subset R \cup (A \circ R)
+ *
+ * Otherwise, we check if the power is exact.
+ *
+ * Note that "app" has an extra input and output coordinate to encode
+ * the length of the part.  If we are only interested in the transitive
+ * closure, then we can simply project out these coordinates first.
+ */
+static int check_exactness(__isl_take isl_map *map, __isl_take isl_map *app,
+       int project)
+{
+       isl_map *test;
+       int exact;
+       unsigned d;
+
+       if (!project)
+               return check_power_exactness(map, app);
+
+       d = isl_map_dim(map, isl_dim_in);
+       app = set_path_length(app, 0, 1);
+       app = isl_map_project_out(app, isl_dim_in, d, 1);
+       app = isl_map_project_out(app, isl_dim_out, d, 1);
+
+       test = isl_map_apply_range(isl_map_copy(map), isl_map_copy(app));
+       test = isl_map_union(test, isl_map_copy(map));
+
+       exact = isl_map_is_subset(app, test);
+
+       isl_map_free(app);
+       isl_map_free(test);
+
+       isl_map_free(map);
+
+       return exact;
+error:
+       isl_map_free(app);
+       isl_map_free(map);
+       return -1;
+}
+
 /*
  * The transitive closure implementation is based on the paper
  * "Computing the Transitive Closure of a Union of Affine Integer
@@ -95,16 +246,22 @@ error:
 #define IMPURE         0
 #define PURE_PARAM     1
 #define PURE_VAR       2
+#define MIXED          3
 
 /* Return PURE_PARAM if only the coefficients of the parameters are non-zero.
  * Return PURE_VAR if only the coefficients of the set variables are non-zero.
+ * Return MIXED if only the coefficients of the parameters and the set
+ *     variables are non-zero and if moreover the parametric constant
+ *     can never attain positive values.
  * Return IMPURE otherwise.
  */
-static int purity(__isl_keep isl_basic_set *bset, isl_int *c)
+static int purity(__isl_keep isl_basic_set *bset, isl_int *c, int eq)
 {
        unsigned d;
        unsigned n_div;
        unsigned nparam;
+       int k;
+       int empty;
 
        n_div = isl_basic_set_dim(bset, isl_dim_div);
        d = isl_basic_set_dim(bset, isl_dim_set);
@@ -116,7 +273,25 @@ static int purity(__isl_keep isl_basic_set *bset, isl_int *c)
                return PURE_VAR;
        if (isl_seq_first_non_zero(c + 1 + nparam, d) == -1)
                return PURE_PARAM;
-       return IMPURE;
+       if (eq)
+               return IMPURE;
+
+       bset = isl_basic_set_copy(bset);
+       bset = isl_basic_set_cow(bset);
+       bset = isl_basic_set_extend_constraints(bset, 0, 1);
+       k = isl_basic_set_alloc_inequality(bset);
+       if (k < 0)
+               goto error;
+       isl_seq_clr(bset->ineq[k], 1 + isl_basic_set_total_dim(bset));
+       isl_seq_cpy(bset->ineq[k], c, 1 + nparam);
+       isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
+       empty = isl_basic_set_is_empty(bset);
+       isl_basic_set_free(bset);
+
+       return empty < 0 ? -1 : empty ? MIXED : IMPURE;
+error:
+       isl_basic_set_free(bset);
+       return -1;
 }
 
 /* Given a set of offsets "delta", construct a relation of the
@@ -136,12 +311,16 @@ static int purity(__isl_keep isl_basic_set *bset, isl_int *c)
  * In particular, let delta be defined as
  *
  *     \delta = [p] -> { [x] : A x + a >= and B p + b >= 0 and
- *                             C x + C'p + c >= 0 }
+ *                             C x + C'p + c >= 0 and
+ *                             D x + D'p + d >= 0 }
  *
- * then the relation is constructed as
+ * where the constraints C x + C'p + c >= 0 are such that the parametric
+ * constant term of each constraint j, "C_j x + C'_j p + c_j",
+ * can never attain positive values, then the relation is constructed as
  *
  *     { [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and
- *                     A f + k a >= 0 and B p + b >= 0 and k >= 1 }
+ *                     A f + k a >= 0 and B p + b >= 0 and
+ *                     C f + C'p + c >= 0 and k >= 1 }
  *     union { [x] -> [x] }
  *
  * Existentially quantified variables in \delta are currently ignored.
@@ -184,7 +363,9 @@ static __isl_give isl_map *path_along_delta(__isl_take isl_dim *dim,
        }
 
        for (i = 0; i < delta->n_eq; ++i) {
-               int p = purity(delta, delta->eq[i]);
+               int p = purity(delta, delta->eq[i], 1);
+               if (p < 0)
+                       goto error;
                if (p == IMPURE)
                        continue;
                k = isl_basic_map_alloc_equality(path);
@@ -200,7 +381,9 @@ static __isl_give isl_map *path_along_delta(__isl_take isl_dim *dim,
        }
 
        for (i = 0; i < delta->n_ineq; ++i) {
-               int p = purity(delta, delta->ineq[i]);
+               int p = purity(delta, delta->ineq[i], 0);
+               if (p < 0)
+                       goto error;
                if (p == IMPURE)
                        continue;
                k = isl_basic_map_alloc_inequality(path);
@@ -211,8 +394,13 @@ static __isl_give isl_map *path_along_delta(__isl_take isl_dim *dim,
                        isl_seq_cpy(path->ineq[k] + off,
                                    delta->ineq[i] + 1 + nparam, d);
                        isl_int_set(path->ineq[k][off + d], delta->ineq[i][0]);
-               } else
+               } else if (p == PURE_PARAM) {
+                       isl_seq_cpy(path->ineq[k], delta->ineq[i], 1 + nparam);
+               } else {
+                       isl_seq_cpy(path->ineq[k] + off,
+                                   delta->ineq[i] + 1 + nparam, d);
                        isl_seq_cpy(path->ineq[k], delta->ineq[i], 1 + nparam);
+               }
        }
 
        k = isl_basic_map_alloc_inequality(path);
@@ -368,6 +556,7 @@ static __isl_give isl_map *construct_extended_path(__isl_take isl_dim *dim,
                if (j < d) {
                        path = isl_map_apply_range(path,
                                path_along_delta(isl_dim_copy(dim), delta));
+                       path = isl_map_coalesce(path);
                } else {
                        isl_basic_set_free(delta);
                        ++n;
@@ -417,10 +606,11 @@ error:
  *     { (x) -> (x) }
  */
 static __isl_give isl_map *construct_component(__isl_take isl_dim *dim,
-       __isl_keep isl_map *map, int *project)
+       __isl_keep isl_map *map, int *exact, int project)
 {
        struct isl_set *domain = NULL;
        struct isl_set *range = NULL;
+       struct isl_set *overlap;
        struct isl_map *app = NULL;
        struct isl_map *path = NULL;
 
@@ -428,14 +618,38 @@ static __isl_give isl_map *construct_component(__isl_take isl_dim *dim,
        domain = isl_set_coalesce(domain);
        range = isl_map_range(isl_map_copy(map));
        range = isl_set_coalesce(range);
+       overlap = isl_set_intersect(isl_set_copy(domain), isl_set_copy(range));
+       if (isl_set_is_empty(overlap) == 1) {
+               isl_set_free(domain);
+               isl_set_free(range);
+               isl_set_free(overlap);
+               isl_dim_free(dim);
+
+               map = isl_map_copy(map);
+               map = isl_map_add(map, isl_dim_in, 1);
+               map = isl_map_add(map, isl_dim_out, 1);
+               map = set_path_length(map, 1, 1);
+               return map;
+       }
+       isl_set_free(overlap);
        app = isl_map_from_domain_and_range(domain, range);
        app = isl_map_add(app, isl_dim_in, 1);
        app = isl_map_add(app, isl_dim_out, 1);
 
-       path = construct_extended_path(isl_dim_copy(dim), map, project);
+       path = construct_extended_path(isl_dim_copy(dim), map,
+                                       exact && *exact ? &project : NULL);
        app = isl_map_intersect(app, path);
 
+       if (exact && *exact &&
+           (*exact = check_exactness(isl_map_copy(map), isl_map_copy(app),
+                                     project)) < 0)
+               goto error;
+
        return isl_map_union(app, isl_map_identity(isl_dim_domain(dim)));
+error:
+       isl_dim_free(dim);
+       isl_map_free(app);
+       return NULL;
 }
 
 /* Structure for representing the nodes in the graph being traversed
@@ -518,31 +732,28 @@ error:
  *
  *     R_1 \circ R_2
  *
- * is non-empty and that moreover, it is non-empty on the set
- * of elements that do not get mapped to the same set of elements
- * by both "R_1 \circ R_2" and "R_2 \circ R_1".
- * For elements that do get mapped to the same elements by these
- * two compositions, R_1 and R_2 are commutative, so if these
- * elements are the only ones for which R_1 \circ R_2 is non-empty,
- * then you may just as well apply R_1 first.
+ * is a subset of
+ *
+ *     R_2 \circ R_1
+ *
+ * If so, then there is no reason for R_1 to immediately follow R_2
+ * in any path.
  */
 static int basic_map_follows(__isl_keep isl_basic_map *bmap1,
        __isl_keep isl_basic_map *bmap2)
 {
        struct isl_map *map12 = NULL;
        struct isl_map *map21 = NULL;
-       struct isl_map *d = NULL;
-       struct isl_set *dom = NULL;
-       int empty;
+       int subset;
 
        map21 = isl_map_from_basic_map(
                        isl_basic_map_apply_range(
                                isl_basic_map_copy(bmap2),
                                isl_basic_map_copy(bmap1)));
-       empty = isl_map_is_empty(map21);
-       if (empty < 0)
+       subset = isl_map_is_empty(map21);
+       if (subset < 0)
                goto error;
-       if (empty) {
+       if (subset) {
                isl_map_free(map21);
                return 0;
        }
@@ -551,18 +762,13 @@ static int basic_map_follows(__isl_keep isl_basic_map *bmap1,
                        isl_basic_map_apply_range(
                                isl_basic_map_copy(bmap1),
                                isl_basic_map_copy(bmap2)));
-       d = isl_map_subtract(isl_map_copy(map12), isl_map_copy(map21));
-       d = isl_map_union(d,
-               isl_map_subtract(isl_map_copy(map21), isl_map_copy(map12)));
-       dom = isl_map_domain(d);
 
-       map21 = isl_map_intersect_domain(map21, dom);
-       empty = isl_map_is_empty(map21);
+       subset = isl_map_is_subset(map21, map12);
 
        isl_map_free(map12);
        isl_map_free(map21);
 
-       return empty < 0 ? -1 : !empty;
+       return subset < 0 ? -1 : !subset;
 error:
        isl_map_free(map21);
        return -1;
@@ -651,7 +857,7 @@ static int power_components_tarjan(struct basic_map_sort *s,
  * to be strictly positive.
  */
 static __isl_give isl_map *construct_power_components(__isl_take isl_dim *dim,
-       __isl_keep isl_map *map, int *project)
+       __isl_keep isl_map *map, int *exact, int project)
 {
        int i, n;
        struct isl_map *path = NULL;
@@ -660,7 +866,7 @@ static __isl_give isl_map *construct_power_components(__isl_take isl_dim *dim,
        if (!map)
                goto error;
        if (map->n <= 1)
-               return construct_component(dim, map, project);
+               return construct_component(dim, map, exact, project);
 
        s = basic_map_sort_alloc(map->ctx, map->n);
        if (!s)
@@ -685,7 +891,8 @@ static __isl_give isl_map *construct_power_components(__isl_take isl_dim *dim,
                        ++i;
                }
                path = isl_map_apply_range(path,
-                           construct_component(isl_dim_copy(dim), comp, project));
+                           construct_component(isl_dim_copy(dim), comp,
+                                               exact, project));
                isl_map_free(comp);
                ++i;
        }
@@ -725,7 +932,7 @@ error:
  * and made positive.  The extra coordinates are subsequently projected out.
  */
 static __isl_give isl_map *construct_power(__isl_keep isl_map *map,
-       unsigned param, int *project)
+       unsigned param, int *exact, int project)
 {
        struct isl_map *app = NULL;
        struct isl_map *diff;
@@ -741,7 +948,8 @@ static __isl_give isl_map *construct_power(__isl_keep isl_map *map,
        dim = isl_dim_add(dim, isl_dim_in, 1);
        dim = isl_dim_add(dim, isl_dim_out, 1);
 
-       app = construct_power_components(isl_dim_copy(dim), map, project);
+       app = construct_power_components(isl_dim_copy(dim), map,
+                                       exact, project);
 
        diff = equate_parameter_to_length(dim, param);
        app = isl_map_intersect(app, diff);
@@ -751,148 +959,6 @@ static __isl_give isl_map *construct_power(__isl_keep isl_map *map,
        return app;
 }
 
-/* Shift variable at position "pos" up by one.
- * That is, replace the corresponding variable v by v - 1.
- */
-static __isl_give isl_basic_map *basic_map_shift_pos(
-       __isl_take isl_basic_map *bmap, unsigned pos)
-{
-       int i;
-
-       bmap = isl_basic_map_cow(bmap);
-       if (!bmap)
-               return NULL;
-
-       for (i = 0; i < bmap->n_eq; ++i)
-               isl_int_sub(bmap->eq[i][0], bmap->eq[i][0], bmap->eq[i][pos]);
-
-       for (i = 0; i < bmap->n_ineq; ++i)
-               isl_int_sub(bmap->ineq[i][0],
-                               bmap->ineq[i][0], bmap->ineq[i][pos]);
-
-       for (i = 0; i < bmap->n_div; ++i) {
-               if (isl_int_is_zero(bmap->div[i][0]))
-                       continue;
-               isl_int_sub(bmap->div[i][1],
-                               bmap->div[i][1], bmap->div[i][1 + pos]);
-       }
-
-       return bmap;
-}
-
-/* Shift variable at position "pos" up by one.
- * That is, replace the corresponding variable v by v - 1.
- */
-static __isl_give isl_map *map_shift_pos(__isl_take isl_map *map, unsigned pos)
-{
-       int i;
-
-       map = isl_map_cow(map);
-       if (!map)
-               return NULL;
-
-       for (i = 0; i < map->n; ++i) {
-               map->p[i] = basic_map_shift_pos(map->p[i], pos);
-               if (!map->p[i])
-                       goto error;
-       }
-       ISL_F_CLR(map, ISL_MAP_NORMALIZED);
-       return map;
-error:
-       isl_map_free(map);
-       return NULL;
-}
-
-/* Check whether the overapproximation of the power of "map" is exactly
- * the power of "map".  Let R be "map" and A_k the overapproximation.
- * The approximation is exact if
- *
- *     A_1 = R
- *     A_k = A_{k-1} \circ R                   k >= 2
- *
- * Since A_k is known to be an overapproximation, we only need to check
- *
- *     A_1 \subset R
- *     A_k \subset A_{k-1} \circ R             k >= 2
- *
- */
-static int check_power_exactness(__isl_take isl_map *map,
-       __isl_take isl_map *app, unsigned param)
-{
-       int exact;
-       isl_map *app_1;
-       isl_map *app_2;
-
-       app_1 = isl_map_fix_si(isl_map_copy(app), isl_dim_param, param, 1);
-
-       exact = isl_map_is_subset(app_1, map);
-       isl_map_free(app_1);
-
-       if (!exact || exact < 0) {
-               isl_map_free(app);
-               isl_map_free(map);
-               return exact;
-       }
-
-       app_2 = isl_map_lower_bound_si(isl_map_copy(app),
-                                       isl_dim_param, param, 2);
-       app_1 = map_shift_pos(app, 1 + param);
-       app_1 = isl_map_apply_range(map, app_1);
-
-       exact = isl_map_is_subset(app_2, app_1);
-
-       isl_map_free(app_1);
-       isl_map_free(app_2);
-
-       return exact;
-}
-
-/* Check whether the overapproximation of the power of "map" is exactly
- * the power of "map", possibly after projecting out the power (if "project"
- * is set).
- *
- * If "project" is set and if "steps" can only result in acyclic paths,
- * then we check
- *
- *     A = R \cup (A \circ R)
- *
- * where A is the overapproximation with the power projected out, i.e.,
- * an overapproximation of the transitive closure.
- * More specifically, since A is known to be an overapproximation, we check
- *
- *     A \subset R \cup (A \circ R)
- *
- * Otherwise, we check if the power is exact.
- */
-static int check_exactness(__isl_take isl_map *map, __isl_take isl_map *app,
-       unsigned param, int project)
-{
-       isl_map *test;
-       int exact;
-
-       if (!project)
-               return check_power_exactness(map, app, param);
-
-       map = isl_map_project_out(map, isl_dim_param, param, 1);
-       app = isl_map_project_out(app, isl_dim_param, param, 1);
-
-       test = isl_map_apply_range(isl_map_copy(map), isl_map_copy(app));
-       test = isl_map_union(test, isl_map_copy(map));
-
-       exact = isl_map_is_subset(app, test);
-
-       isl_map_free(app);
-       isl_map_free(test);
-
-       isl_map_free(map);
-
-       return exact;
-error:
-       isl_map_free(app);
-       isl_map_free(map);
-       return -1;
-}
-
 /* Compute the positive powers of "map", or an overapproximation.
  * The power is given by parameter "param".  If the result is exact,
  * then *exact is set to 1.
@@ -919,12 +985,7 @@ static __isl_give isl_map *map_power(__isl_take isl_map *map, unsigned param,
                isl_map_dim(map, isl_dim_in) == isl_map_dim(map, isl_dim_out),
                goto error);
 
-       app = construct_power(map, param, exact ? &project : NULL);
-
-       if (exact &&
-           (*exact = check_exactness(isl_map_copy(map), isl_map_copy(app),
-                                     param, project)) < 0)
-               goto error;
+       app = construct_power(map, param, exact, project);
 
        isl_map_free(map);
        return app;