change isl_map_power interface
[platform/upstream/isl.git] / isl_transitive_closure.c
index 7e65915..37b5142 100644 (file)
@@ -8,10 +8,13 @@
  * 91893 Orsay, France 
  */
 
-#include "isl_map.h"
-#include "isl_map_private.h"
-#include "isl_seq.h"
-#include <isl_lp.h>
+#include <isl_map_private.h>
+#include <isl/map.h>
+#include <isl/seq.h>
+#include <isl_dim_private.h>
+#include <isl/lp.h>
+#include <isl/union_map.h>
+#include <isl_mat_private.h>
 
 int isl_map_is_transitively_closed(__isl_keep isl_map *map)
 {
@@ -24,6 +27,19 @@ int isl_map_is_transitively_closed(__isl_keep isl_map *map)
 
        return closed;
 }
+
+int isl_union_map_is_transitively_closed(__isl_keep isl_union_map *umap)
+{
+       isl_union_map *umap2;
+       int closed;
+
+       umap2 = isl_union_map_apply_range(isl_union_map_copy(umap),
+                                         isl_union_map_copy(umap));
+       closed = isl_union_map_is_subset(umap2, umap);
+       isl_union_map_free(umap2);
+
+       return closed;
+}
  
 /* Given a map that represents a path with the length of the path
  * encoded as the difference between the last output coordindate
@@ -96,8 +112,8 @@ static int check_power_exactness(__isl_take isl_map *map,
        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 = isl_map_add_dims(map, isl_dim_in, 1);
+       map = isl_map_add_dims(map, isl_dim_out, 1);
        map = set_path_length(map, 1, 1);
 
        app_1 = set_path_length(isl_map_copy(app), 1, 1);
@@ -159,6 +175,8 @@ static int check_exactness(__isl_take isl_map *map, __isl_take isl_map *app,
        app = isl_map_project_out(app, isl_dim_in, d, 1);
        app = isl_map_project_out(app, isl_dim_out, d, 1);
 
+       app = isl_map_reset_dim(app, isl_map_get_dim(map));
+
        test = isl_map_apply_range(isl_map_copy(map), isl_map_copy(app));
        test = isl_map_union(test, isl_map_copy(map));
 
@@ -304,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)
@@ -315,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);
@@ -397,7 +421,7 @@ static __isl_give int *get_div_purity(__isl_keep isl_basic_set *bset)
  * check if setting the length to zero results in only the identity
  * mapping.
  */
-int empty_path_is_identity(__isl_keep isl_basic_map *path, unsigned pos)
+static int empty_path_is_identity(__isl_keep isl_basic_map *path, unsigned pos)
 {
        isl_basic_map *test = NULL;
        isl_basic_map *id = NULL;
@@ -411,7 +435,7 @@ int empty_path_is_identity(__isl_keep isl_basic_map *path, unsigned pos)
                goto error;
        isl_seq_clr(test->eq[k], 1 + isl_basic_map_total_dim(test));
        isl_int_set_si(test->eq[k][pos], 1);
-       id = isl_basic_map_identity(isl_dim_domain(isl_basic_map_get_dim(path)));
+       id = isl_basic_map_identity(isl_basic_map_get_dim(path));
        is_id = isl_basic_map_is_equal(test, id);
        isl_basic_map_free(test);
        isl_basic_map_free(id);
@@ -421,9 +445,13 @@ error:
        return -1;
 }
 
-__isl_give isl_basic_map *add_delta_constraints(__isl_take isl_basic_map *path,
+/* 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;
@@ -437,6 +465,8 @@ __isl_give isl_basic_map *add_delta_constraints(__isl_take isl_basic_map *path,
                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) {
@@ -486,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 >= and B p + b >= 0 and
  *                             C x + C'p + c >= 0 and
  *                             D x + D'p + d >= 0 }
  *
@@ -513,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)
@@ -525,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;
@@ -556,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)
@@ -578,8 +638,7 @@ static __isl_give isl_map *path_along_delta(__isl_take isl_dim *dim,
                isl_dim_free(dim);
                return isl_map_from_basic_map(path);
        }
-       return isl_basic_map_union(path,
-                               isl_basic_map_identity(isl_dim_domain(dim)));
+       return isl_basic_map_union(path, isl_basic_map_identity(dim));
 error:
        free(div_purity);
        isl_dim_free(dim);
@@ -694,7 +753,7 @@ static __isl_give isl_map *construct_extended_path(__isl_take isl_dim *dim,
 
        d = isl_map_dim(map, isl_dim_in);
 
-       path = isl_map_identity(isl_dim_domain(isl_dim_copy(dim)));
+       path = isl_map_identity(isl_dim_copy(dim));
 
        steps = isl_mat_alloc(map->ctx, map->n, d);
        if (!steps)
@@ -757,6 +816,9 @@ static int isl_set_overlaps(__isl_keep isl_set *set1, __isl_keep isl_set *set2)
        isl_set *i;
        int no_overlap;
 
+       if (!isl_dim_tuple_match(set1->dim, isl_dim_set, set2->dim, isl_dim_set))
+               return 0;
+
        i = isl_set_intersect(isl_set_copy(set1), isl_set_copy(set2));
        no_overlap = isl_set_is_empty(i);
        isl_set_free(i);
@@ -801,14 +863,14 @@ static __isl_give isl_map *construct_component(__isl_take isl_dim *dim,
                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 = isl_map_add_dims(map, isl_dim_in, 1);
+               map = isl_map_add_dims(map, isl_dim_out, 1);
                map = set_path_length(map, 1, 1);
                return map;
        }
        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);
+       app = isl_map_add_dims(app, isl_dim_in, 1);
+       app = isl_map_add_dims(app, isl_dim_out, 1);
 
        path = construct_extended_path(isl_dim_copy(dim), map,
                                        exact && *exact ? &project : NULL);
@@ -863,7 +925,7 @@ static __isl_give isl_map *q_closure(__isl_take isl_dim *dim,
        isl_map *map;
        isl_map *app;
 
-       dom = isl_set_add(dom, isl_dim_set, 1);
+       dom = isl_set_add_dims(dom, isl_dim_set, 1);
        app = isl_map_from_domain_and_range(dom, isl_set_copy(dom));
        map = isl_map_from_basic_map(isl_basic_map_copy(bmap));
        path = construct_extended_path(dim, map, &project);
@@ -1013,6 +1075,13 @@ static int composability(__isl_keep isl_set *C, int i,
        return ok;
 }
 
+static __isl_give isl_map *anonymize(__isl_take isl_map *map)
+{
+       map = isl_map_reset(map, isl_dim_in);
+       map = isl_map_reset(map, isl_dim_out);
+       return map;
+}
+
 /* Return a map that is a union of the basic maps in "map", except i,
  * composed to left and right with qc based on the entries of "left"
  * and "right".
@@ -1031,6 +1100,7 @@ static __isl_give isl_map *compose(__isl_keep isl_map *map, int i,
                        continue;
 
                map_j = isl_map_from_basic_map(isl_basic_map_copy(map->p[j]));
+               map_j = anonymize(map_j);
                if (left && left[j])
                        map_j = isl_map_apply_range(map_j, isl_map_copy(qc));
                if (right && right[j])
@@ -1090,7 +1160,7 @@ static __isl_give isl_map *compute_incremental(
 
        if (!left || !right)
                rtc = isl_map_union(isl_map_copy(tc),
-                       isl_map_identity(isl_dim_domain(isl_map_get_dim(tc))));
+                                   isl_map_identity(isl_map_get_dim(tc)));
        if (!right)
                qc = isl_map_apply_range(rtc, qc);
        if (!left)
@@ -1363,9 +1433,9 @@ static int merge(isl_set **set, int *group, __isl_take isl_set *dom, int pos)
                        continue;
 
                set[i] = isl_set_union(set[i], set[group[pos]]);
+               set[group[pos]] = NULL;
                if (!set[i])
                        goto error;
-               set[group[pos]] = NULL;
                group[group[pos]] = i;
                group[pos] = i;
        }
@@ -1420,6 +1490,48 @@ static int add_length(__isl_keep isl_map *map, isl_map ***grid, int n)
        return 0;
 }
 
+/* The core of the Floyd-Warshall algorithm.
+ * Updates the given n x x matrix of relations in place.
+ *
+ * The algorithm iterates over all vertices.  In each step, the whole
+ * matrix is updated to include all paths that go to the current vertex,
+ * possibly stay there a while (including passing through earlier vertices)
+ * and then come back.  At the start of each iteration, the diagonal
+ * element corresponding to the current vertex is replaced by its
+ * transitive closure to account for all indirect paths that stay
+ * in the current vertex.
+ */
+static void floyd_warshall_iterate(isl_map ***grid, int n, int *exact)
+{
+       int r, p, q;
+
+       for (r = 0; r < n; ++r) {
+               int r_exact;
+               grid[r][r] = isl_map_transitive_closure(grid[r][r],
+                               (exact && *exact) ? &r_exact : NULL);
+               if (exact && *exact && !r_exact)
+                       *exact = 0;
+
+               for (p = 0; p < n; ++p)
+                       for (q = 0; q < n; ++q) {
+                               isl_map *loop;
+                               if (p == r && q == r)
+                                       continue;
+                               loop = isl_map_apply_range(
+                                               isl_map_copy(grid[p][r]),
+                                               isl_map_copy(grid[r][q]));
+                               grid[p][q] = isl_map_union(grid[p][q], loop);
+                               loop = isl_map_apply_range(
+                                               isl_map_copy(grid[p][r]),
+                                       isl_map_apply_range(
+                                               isl_map_copy(grid[r][r]),
+                                               isl_map_copy(grid[r][q])));
+                               grid[p][q] = isl_map_union(grid[p][q], loop);
+                               grid[p][q] = isl_map_coalesce(grid[p][q]);
+                       }
+       }
+}
+
 /* Given a partition of the domains and ranges of the basic maps in "map",
  * apply the Floyd-Warshall algorithm with the elements in the partition
  * as vertices.
@@ -1431,14 +1543,6 @@ static int add_length(__isl_keep isl_map *map, isl_map ***grid, int n)
  * apply Floyd-Warshall on this matrix of relations and then take the
  * union of all entries in the matrix as the final result.
  *
- * The algorithm iterates over all vertices.  In each step, the whole
- * matrix is updated to include all paths that go to the current vertex,
- * possibly stay there a while (including passing through earlier vertices)
- * and then come back.  At the start of each iteration, the diagonal
- * element corresponding to the current vertex is replaced by its
- * transitive closure to account for all indirect paths that stay
- * in the current vertex.
- *
  * If we are actually computing the power instead of the transitive closure,
  * i.e., when "project" is not set, then the result should have the
  * path lengths encoded as the difference between an extra pair of
@@ -1451,7 +1555,6 @@ static __isl_give isl_map *floyd_warshall_with_groups(__isl_take isl_dim *dim,
        __isl_keep isl_map *map, int *exact, int project, int *group, int n)
 {
        int i, j, k;
-       int r, p, q;
        isl_map ***grid = NULL;
        isl_map *app;
 
@@ -1485,31 +1588,7 @@ static __isl_give isl_map *floyd_warshall_with_groups(__isl_take isl_dim *dim,
        if (!project && add_length(map, grid, n) < 0)
                goto error;
 
-       for (r = 0; r < n; ++r) {
-               int r_exact;
-               grid[r][r] = isl_map_transitive_closure(grid[r][r],
-                               (exact && *exact) ? &r_exact : NULL);
-               if (exact && *exact && !r_exact)
-                       *exact = 0;
-
-               for (p = 0; p < n; ++p)
-                       for (q = 0; q < n; ++q) {
-                               isl_map *loop;
-                               if (p == r && q == r)
-                                       continue;
-                               loop = isl_map_apply_range(
-                                               isl_map_copy(grid[p][r]),
-                                               isl_map_copy(grid[r][q]));
-                               grid[p][q] = isl_map_union(grid[p][q], loop);
-                               loop = isl_map_apply_range(
-                                               isl_map_copy(grid[p][r]),
-                                       isl_map_apply_range(
-                                               isl_map_copy(grid[r][r]),
-                                               isl_map_copy(grid[r][q])));
-                               grid[p][q] = isl_map_union(grid[p][q], loop);
-                               grid[p][q] = isl_map_coalesce(grid[p][q]);
-                       }
-       }
+       floyd_warshall_iterate(grid, n, exact);
 
        app = isl_map_empty(isl_map_get_dim(map));
 
@@ -1539,13 +1618,8 @@ error:
        return NULL;
 }
 
-/* Check if the domains and ranges of the basic maps in "map" can
- * be partitioned, and if so, apply Floyd-Warshall on the elements
- * of the partition.  Note that we also apply this algorithm
- * if we want to compute the power, i.e., when "project" is not set.
- * However, the results are unlikely to be exact since the recursive
- * calls inside the Floyd-Warshall algorithm typically result in
- * non-linear path lengths quite quickly.
+/* Partition the domains and ranges of the n basic relations in list
+ * into disjoint cells.
  *
  * To find the partition, we simply consider all of the domains
  * and ranges in turn and combine those that overlap.
@@ -1558,44 +1632,81 @@ error:
  * ranges in the corresponding group, or is equal to some l < k,
  * with l another domain or range in the same group.
  */
-static __isl_give isl_map *floyd_warshall(__isl_take isl_dim *dim,
-       __isl_keep isl_map *map, int *exact, int project)
+static int *setup_groups(isl_ctx *ctx, __isl_keep isl_basic_map **list, int n,
+       isl_set ***set, int *n_group)
 {
        int i;
-       isl_set **set = NULL;
        int *group = NULL;
-       int n;
+       int g;
 
-       if (!map)
-               goto error;
-       if (map->n <= 1)
-               return incremental_closure(dim, map, exact, project);
-
-       set = isl_calloc_array(map->ctx, isl_set *, 2 * map->n);
-       group = isl_alloc_array(map->ctx, int, 2 * map->n);
+       *set = isl_calloc_array(ctx, isl_set *, 2 * n);
+       group = isl_alloc_array(ctx, int, 2 * n);
 
-       if (!set || !group)
+       if (!*set || !group)
                goto error;
 
-       for (i = 0; i < map->n; ++i) {
+       for (i = 0; i < n; ++i) {
                isl_set *dom;
                dom = isl_set_from_basic_set(isl_basic_map_domain(
-                               isl_basic_map_copy(map->p[i])));
-               if (merge(set, group, dom, 2 * i) < 0)
+                               isl_basic_map_copy(list[i])));
+               if (merge(*set, group, dom, 2 * i) < 0)
                        goto error;
                dom = isl_set_from_basic_set(isl_basic_map_range(
-                               isl_basic_map_copy(map->p[i])));
-               if (merge(set, group, dom, 2 * i + 1) < 0)
+                               isl_basic_map_copy(list[i])));
+               if (merge(*set, group, dom, 2 * i + 1) < 0)
                        goto error;
        }
 
-       n = 0;
-       for (i = 0; i < 2 * map->n; ++i)
-               if (group[i] == i)
-                       group[i] = n++;
-               else
+       g = 0;
+       for (i = 0; i < 2 * n; ++i)
+               if (group[i] == i) {
+                       if (g != i) {
+                               (*set)[g] = (*set)[i];
+                               (*set)[i] = NULL;
+                       }
+                       group[i] = g++;
+               } else
                        group[i] = group[group[i]];
 
+       *n_group = g;
+
+       return group;
+error:
+       if (*set) {
+               for (i = 0; i < 2 * n; ++i)
+                       isl_set_free((*set)[i]);
+               free(*set);
+               *set = NULL;
+       }
+       free(group);
+       return NULL;
+}
+
+/* Check if the domains and ranges of the basic maps in "map" can
+ * be partitioned, and if so, apply Floyd-Warshall on the elements
+ * of the partition.  Note that we also apply this algorithm
+ * if we want to compute the power, i.e., when "project" is not set.
+ * However, the results are unlikely to be exact since the recursive
+ * calls inside the Floyd-Warshall algorithm typically result in
+ * non-linear path lengths quite quickly.
+ */
+static __isl_give isl_map *floyd_warshall(__isl_take isl_dim *dim,
+       __isl_keep isl_map *map, int *exact, int project)
+{
+       int i;
+       isl_set **set = NULL;
+       int *group = NULL;
+       int n;
+
+       if (!map)
+               goto error;
+       if (map->n <= 1)
+               return incremental_closure(dim, map, exact, project);
+
+       group = setup_groups(map->ctx, map->p, map->n, &set, &n);
+       if (!group)
+               goto error;
+
        for (i = 0; i < 2 * map->n; ++i)
                isl_set_free(set[i]);
 
@@ -1603,10 +1714,6 @@ static __isl_give isl_map *floyd_warshall(__isl_take isl_dim *dim,
 
        return floyd_warshall_with_groups(dim, map, exact, project, group, n);
 error:
-       for (i = 0; i < 2 * map->n; ++i)
-               isl_set_free(set[i]);
-       free(set);
-       free(group);
        isl_dim_free(dim);
        return NULL;
 }
@@ -1714,6 +1821,9 @@ static int basic_map_follows(__isl_keep isl_basic_map *bmap1,
        struct isl_map *map21 = NULL;
        int subset;
 
+       if (!isl_dim_tuple_match(bmap1->dim, isl_dim_in, bmap2->dim, isl_dim_out))
+               return 0;
+
        map21 = isl_map_from_basic_map(
                        isl_basic_map_apply_range(
                                isl_basic_map_copy(bmap2),
@@ -1726,6 +1836,12 @@ static int basic_map_follows(__isl_keep isl_basic_map *bmap1,
                return 0;
        }
 
+       if (!isl_dim_tuple_match(bmap1->dim, isl_dim_in, bmap1->dim, isl_dim_out) ||
+           !isl_dim_tuple_match(bmap2->dim, isl_dim_in, bmap2->dim, isl_dim_out)) {
+               isl_map_free(map21);
+               return 1;
+       }
+
        map12 = isl_map_from_basic_map(
                        isl_basic_map_apply_range(
                                isl_basic_map_copy(bmap1),
@@ -1751,7 +1867,7 @@ error:
  * to be applied after the second.
  */
 static int power_components_tarjan(struct basic_map_sort *s,
-       __isl_keep isl_map *map, int i)
+       __isl_keep isl_basic_map **list, int i)
 {
        int j;
 
@@ -1771,14 +1887,14 @@ static int power_components_tarjan(struct basic_map_sort *s,
                         s->node[j].index > s->node[i].min_index))
                        continue;
 
-               f = basic_map_follows(map->p[i], map->p[j], &s->check_closed);
+               f = basic_map_follows(list[i], list[j], &s->check_closed);
                if (f < 0)
                        return -1;
                if (!f)
                        continue;
 
                if (s->node[j].index < 0) {
-                       power_components_tarjan(s, map, j);
+                       power_components_tarjan(s, list, j);
                        if (s->node[j].min_index < s->node[i].min_index)
                                s->node[i].min_index = s->node[j].min_index;
                } else if (s->node[j].index < s->node[i].min_index)
@@ -1798,6 +1914,31 @@ static int power_components_tarjan(struct basic_map_sort *s,
        return 0;
 }
 
+/* Decompose the "len" basic relations in "list" into strongly connected
+ * components.
+ */
+static struct basic_map_sort *basic_map_sort_init(isl_ctx *ctx, int len,
+       __isl_keep isl_basic_map **list)
+{
+       int i;
+       struct basic_map_sort *s = NULL;
+
+       s = basic_map_sort_alloc(ctx, len);
+       if (!s)
+               return NULL;
+       for (i = len - 1; i >= 0; --i) {
+               if (s->node[i].index >= 0)
+                       continue;
+               if (power_components_tarjan(s, list, i) < 0)
+                       goto error;
+       }
+
+       return s;
+error:
+       basic_map_sort_free(s);
+       return NULL;
+}
+
 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
  * and a dimension specification (Z^{n+1} -> Z^{n+1}),
  * construct a map that is an overapproximation of the map
@@ -1845,15 +1986,9 @@ static __isl_give isl_map *construct_power_components(__isl_take isl_dim *dim,
        if (map->n <= 1)
                return floyd_warshall(dim, map, exact, project);
 
-       s = basic_map_sort_alloc(map->ctx, map->n);
+       s = basic_map_sort_init(map->ctx, map->n, map->p);
        if (!s)
                goto error;
-       for (i = map->n - 1; i >= 0; --i) {
-               if (s->node[i].index >= 0)
-                       continue;
-               if (power_components_tarjan(s, map, i) < 0)
-                       goto error;
-       }
 
        orig_exact = exact;
        if (s->check_closed && !exact)
@@ -1866,6 +2001,7 @@ static __isl_give isl_map *construct_power_components(__isl_take isl_dim *dim,
                path = isl_map_empty(isl_map_get_dim(map));
        else
                path = isl_map_empty(isl_dim_copy(dim));
+       path = anonymize(path);
        while (n) {
                struct isl_map *comp;
                isl_map *path_comp, *path_comb;
@@ -1878,6 +2014,7 @@ static __isl_give isl_map *construct_power_components(__isl_take isl_dim *dim,
                }
                path_comp = floyd_warshall(isl_dim_copy(dim),
                                                comp, exact, project);
+               path_comp = anonymize(path_comp);
                path_comb = isl_map_apply_range(isl_map_copy(path),
                                                isl_map_copy(path_comp));
                path = isl_map_union(path, path_comp);
@@ -1998,46 +2135,105 @@ error:
 }
 
 /* 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.
+ * The result maps the exponent to a nested copy of the corresponding power.
+ * If the result is exact, then *exact is set to 1.
  * map_power constructs an extended relation with the path lengths
  * encoded as the difference between the final coordinates.
- * In the final step, this difference is equated to the parameter "param"
- * and made positive.  The extra coordinates are subsequently projected out.
+ * In the final step, this difference is equated to an extra parameter
+ * and made positive.  The extra coordinates are subsequently projected out
+ * and the parameter is turned into the domain of the result.
  */
-__isl_give isl_map *isl_map_power(__isl_take isl_map *map, unsigned param,
-       int *exact)
+__isl_give isl_map *isl_map_power(__isl_take isl_map *map, int *exact)
 {
+       isl_dim *target_dim;
        isl_dim *dim;
        isl_map *diff;
        unsigned d;
+       unsigned param;
 
        if (!map)
                return NULL;
 
-       isl_assert(map->ctx, param < isl_map_dim(map, isl_dim_param),
-               goto error);
-
        d = isl_map_dim(map, isl_dim_in);
+       param = isl_map_dim(map, isl_dim_param);
 
        map = isl_map_compute_divs(map);
        map = isl_map_coalesce(map);
 
-       if (isl_map_fast_is_empty(map))
+       if (isl_map_fast_is_empty(map)) {
+               map = isl_map_from_range(isl_map_wrap(map));
+               map = isl_map_add_dims(map, isl_dim_in, 1);
+               map = isl_map_set_dim_name(map, isl_dim_in, 0, "k");
                return map;
+       }
+
+       target_dim = isl_map_get_dim(map);
+       target_dim = isl_dim_from_range(isl_dim_wrap(target_dim));
+       target_dim = isl_dim_add(target_dim, isl_dim_in, 1);
+       target_dim = isl_dim_set_name(target_dim, isl_dim_in, 0, "k");
 
        map = map_power(map, exact, 0);
 
+       map = isl_map_add_dims(map, isl_dim_param, 1);
        dim = isl_map_get_dim(map);
        diff = equate_parameter_to_length(dim, param);
        map = isl_map_intersect(map, diff);
        map = isl_map_project_out(map, isl_dim_in, d, 1);
        map = isl_map_project_out(map, isl_dim_out, d, 1);
+       map = isl_map_from_range(isl_map_wrap(map));
+       map = isl_map_move_dims(map, isl_dim_in, 0, isl_dim_param, param, 1);
+
+       map = isl_map_reset_dim(map, target_dim);
+
+       return map;
+}
+
+/* Compute a relation that maps each element in the range of the input
+ * relation to the lengths of all paths composed of edges in the input
+ * relation that end up in the given range element.
+ * The result may be an overapproximation, in which case *exact is set to 0.
+ * The resulting relation is very similar to the power relation.
+ * The difference are that the domain has been projected out, the
+ * range has become the domain and the exponent is the range instead
+ * of a parameter.
+ */
+__isl_give isl_map *isl_map_reaching_path_lengths(__isl_take isl_map *map,
+       int *exact)
+{
+       isl_dim *dim;
+       isl_map *diff;
+       unsigned d;
+       unsigned param;
+
+       if (!map)
+               return NULL;
+
+       d = isl_map_dim(map, isl_dim_in);
+       param = isl_map_dim(map, isl_dim_param);
+
+       map = isl_map_compute_divs(map);
+       map = isl_map_coalesce(map);
+
+       if (isl_map_fast_is_empty(map)) {
+               if (exact)
+                       *exact = 1;
+               map = isl_map_project_out(map, isl_dim_out, 0, d);
+               map = isl_map_add_dims(map, isl_dim_out, 1);
+               return map;
+       }
+
+       map = map_power(map, exact, 0);
+
+       map = isl_map_add_dims(map, isl_dim_param, 1);
+       dim = isl_map_get_dim(map);
+       diff = equate_parameter_to_length(dim, param);
+       map = isl_map_intersect(map, diff);
+       map = isl_map_project_out(map, isl_dim_in, 0, d + 1);
+       map = isl_map_project_out(map, isl_dim_out, d, 1);
+       map = isl_map_reverse(map);
+       map = isl_map_move_dims(map, isl_dim_out, 0, isl_dim_param, param, 1);
 
        return map;
-error:
-       isl_map_free(map);
-       return NULL;
 }
 
 /* Check whether equality i of bset is a pure stride constraint
@@ -2507,12 +2703,13 @@ error:
 __isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
        int *exact)
 {
+       isl_dim *target_dim;
        int closed;
 
        if (!map)
                goto error;
 
-       if (map->ctx->opt->closure == ISL_CLOSURE_OMEGA)
+       if (map->ctx->opt->closure == ISL_CLOSURE_BOX)
                return transitive_closure_omega(map, exact);
 
        map = isl_map_compute_divs(map);
@@ -2526,10 +2723,297 @@ __isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
                return map;
        }
 
+       target_dim = isl_map_get_dim(map);
        map = map_power(map, exact, 1);
+       map = isl_map_reset_dim(map, target_dim);
 
        return map;
 error:
        isl_map_free(map);
        return NULL;
 }
+
+static int inc_count(__isl_take isl_map *map, void *user)
+{
+       int *n = user;
+
+       *n += map->n;
+
+       isl_map_free(map);
+
+       return 0;
+}
+
+static int collect_basic_map(__isl_take isl_map *map, void *user)
+{
+       int i;
+       isl_basic_map ***next = user;
+
+       for (i = 0; i < map->n; ++i) {
+               **next = isl_basic_map_copy(map->p[i]);
+               if (!**next)
+                       goto error;
+               (*next)++;
+       }
+
+       isl_map_free(map);
+       return 0;
+error:
+       isl_map_free(map);
+       return -1;
+}
+
+/* Perform Floyd-Warshall on the given list of basic relations.
+ * The basic relations may live in different dimensions,
+ * but basic relations that get assigned to the diagonal of the
+ * grid have domains and ranges of the same dimension and so
+ * the standard algorithm can be used because the nested transitive
+ * closures are only applied to diagonal elements and because all
+ * compositions are peformed on relations with compatible domains and ranges.
+ */
+static __isl_give isl_union_map *union_floyd_warshall_on_list(isl_ctx *ctx,
+       __isl_keep isl_basic_map **list, int n, int *exact)
+{
+       int i, j, k;
+       int n_group;
+       int *group = NULL;
+       isl_set **set = NULL;
+       isl_map ***grid = NULL;
+       isl_union_map *app;
+
+       group = setup_groups(ctx, list, n, &set, &n_group);
+       if (!group)
+               goto error;
+
+       grid = isl_calloc_array(ctx, isl_map **, n_group);
+       if (!grid)
+               goto error;
+       for (i = 0; i < n_group; ++i) {
+               grid[i] = isl_calloc_array(map->ctx, isl_map *, n_group);
+               if (!grid[i])
+                       goto error;
+               for (j = 0; j < n_group; ++j) {
+                       isl_dim *dim1, *dim2, *dim;
+                       dim1 = isl_dim_reverse(isl_set_get_dim(set[i]));
+                       dim2 = isl_set_get_dim(set[j]);
+                       dim = isl_dim_join(dim1, dim2);
+                       grid[i][j] = isl_map_empty(dim);
+               }
+       }
+
+       for (k = 0; k < n; ++k) {
+               i = group[2 * k];
+               j = group[2 * k + 1];
+               grid[i][j] = isl_map_union(grid[i][j],
+                               isl_map_from_basic_map(
+                                       isl_basic_map_copy(list[k])));
+       }
+       
+       floyd_warshall_iterate(grid, n_group, exact);
+
+       app = isl_union_map_empty(isl_map_get_dim(grid[0][0]));
+
+       for (i = 0; i < n_group; ++i) {
+               for (j = 0; j < n_group; ++j)
+                       app = isl_union_map_add_map(app, grid[i][j]);
+               free(grid[i]);
+       }
+       free(grid);
+
+       for (i = 0; i < 2 * n; ++i)
+               isl_set_free(set[i]);
+       free(set);
+
+       free(group);
+       return app;
+error:
+       if (grid)
+               for (i = 0; i < n_group; ++i) {
+                       if (!grid[i])
+                               continue;
+                       for (j = 0; j < n_group; ++j)
+                               isl_map_free(grid[i][j]);
+                       free(grid[i]);
+               }
+       free(grid);
+       if (set) {
+               for (i = 0; i < 2 * n; ++i)
+                       isl_set_free(set[i]);
+               free(set);
+       }
+       free(group);
+       return NULL;
+}
+
+/* Perform Floyd-Warshall on the given union relation.
+ * The implementation is very similar to that for non-unions.
+ * The main difference is that it is applied unconditionally.
+ * We first extract a list of basic maps from the union map
+ * and then perform the algorithm on this list.
+ */
+static __isl_give isl_union_map *union_floyd_warshall(
+       __isl_take isl_union_map *umap, int *exact)
+{
+       int i, n;
+       isl_ctx *ctx;
+       isl_basic_map **list;
+       isl_basic_map **next;
+       isl_union_map *res;
+
+       n = 0;
+       if (isl_union_map_foreach_map(umap, inc_count, &n) < 0)
+               goto error;
+
+       ctx = isl_union_map_get_ctx(umap);
+       list = isl_calloc_array(ctx, isl_basic_map *, n);
+       if (!list)
+               goto error;
+
+       next = list;
+       if (isl_union_map_foreach_map(umap, collect_basic_map, &next) < 0)
+               goto error;
+
+       res = union_floyd_warshall_on_list(ctx, list, n, exact);
+
+       if (list) {
+               for (i = 0; i < n; ++i)
+                       isl_basic_map_free(list[i]);
+               free(list);
+       }
+
+       isl_union_map_free(umap);
+       return res;
+error:
+       if (list) {
+               for (i = 0; i < n; ++i)
+                       isl_basic_map_free(list[i]);
+               free(list);
+       }
+       isl_union_map_free(umap);
+       return NULL;
+}
+
+/* Decompose the give union relation into strongly connected components.
+ * The implementation is essentially the same as that of
+ * construct_power_components with the major difference that all
+ * operations are performed on union maps.
+ */
+static __isl_give isl_union_map *union_components(
+       __isl_take isl_union_map *umap, int *exact)
+{
+       int i;
+       int n;
+       isl_ctx *ctx;
+       isl_basic_map **list;
+       isl_basic_map **next;
+       isl_union_map *path = NULL;
+       struct basic_map_sort *s = NULL;
+       int c, l;
+       int recheck = 0;
+
+       n = 0;
+       if (isl_union_map_foreach_map(umap, inc_count, &n) < 0)
+               goto error;
+
+       if (n <= 1)
+               return union_floyd_warshall(umap, exact);
+
+       ctx = isl_union_map_get_ctx(umap);
+       list = isl_calloc_array(ctx, isl_basic_map *, n);
+       if (!list)
+               goto error;
+
+       next = list;
+       if (isl_union_map_foreach_map(umap, collect_basic_map, &next) < 0)
+               goto error;
+
+       s = basic_map_sort_init(ctx, n, list);
+       if (!s)
+               goto error;
+
+       c = 0;
+       i = 0;
+       l = n;
+       path = isl_union_map_empty(isl_union_map_get_dim(umap));
+       while (l) {
+               isl_union_map *comp;
+               isl_union_map *path_comp, *path_comb;
+               comp = isl_union_map_empty(isl_union_map_get_dim(umap));
+               while (s->order[i] != -1) {
+                       comp = isl_union_map_add_map(comp,
+                                   isl_map_from_basic_map(
+                                       isl_basic_map_copy(list[s->order[i]])));
+                       --l;
+                       ++i;
+               }
+               path_comp = union_floyd_warshall(comp, exact);
+               path_comb = isl_union_map_apply_range(isl_union_map_copy(path),
+                                               isl_union_map_copy(path_comp));
+               path = isl_union_map_union(path, path_comp);
+               path = isl_union_map_union(path, path_comb);
+               ++i;
+               ++c;
+       }
+
+       if (c > 1 && s->check_closed && !*exact) {
+               int closed;
+
+               closed = isl_union_map_is_transitively_closed(path);
+               if (closed < 0)
+                       goto error;
+               recheck = !closed;
+       }
+
+       basic_map_sort_free(s);
+
+       for (i = 0; i < n; ++i)
+               isl_basic_map_free(list[i]);
+       free(list);
+
+       if (recheck) {
+               isl_union_map_free(path);
+               return union_floyd_warshall(umap, exact);
+       }
+
+       isl_union_map_free(umap);
+
+       return path;
+error:
+       basic_map_sort_free(s);
+       if (list) {
+               for (i = 0; i < n; ++i)
+                       isl_basic_map_free(list[i]);
+               free(list);
+       }
+       isl_union_map_free(umap);
+       isl_union_map_free(path);
+       return NULL;
+}
+
+/* Compute the transitive closure  of "umap", or an overapproximation.
+ * If the result is exact, then *exact is set to 1.
+ */
+__isl_give isl_union_map *isl_union_map_transitive_closure(
+       __isl_take isl_union_map *umap, int *exact)
+{
+       int closed;
+
+       if (!umap)
+               return NULL;
+
+       if (exact)
+               *exact = 1;
+
+       umap = isl_union_map_compute_divs(umap);
+       umap = isl_union_map_coalesce(umap);
+       closed = isl_union_map_is_transitively_closed(umap);
+       if (closed < 0)
+               goto error;
+       if (closed)
+               return umap;
+       umap = union_components(umap, exact);
+       return umap;
+error:
+       isl_union_map_free(umap);
+       return NULL;
+}