isl_basic_set_opt: avoid invalid access on error path
[platform/upstream/isl.git] / isl_transitive_closure.c
index 7e65915..70d15b9 100644 (file)
@@ -1,17 +1,23 @@
 /*
  * Copyright 2010      INRIA Saclay
  *
- * Use of this software is governed by the GNU LGPLv2.1 license
+ * Use of this software is governed by the MIT license
  *
  * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
  * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
  * 91893 Orsay, France 
  */
 
-#include "isl_map.h"
-#include "isl_map_private.h"
-#include "isl_seq.h"
-#include <isl_lp.h>
+#include <isl_ctx_private.h>
+#include <isl_map_private.h>
+#include <isl/map.h>
+#include <isl/seq.h>
+#include <isl_space_private.h>
+#include <isl/lp.h>
+#include <isl/union_map.h>
+#include <isl_mat_private.h>
+#include <isl_options_private.h>
+#include <isl_tarjan.h>
 
 int isl_map_is_transitively_closed(__isl_keep isl_map *map)
 {
@@ -24,6 +30,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
@@ -34,7 +53,7 @@ int isl_map_is_transitively_closed(__isl_keep isl_map *map)
 static __isl_give isl_map *set_path_length(__isl_take isl_map *map,
        int exactly, int length)
 {
-       struct isl_dim *dim;
+       isl_space *dim;
        struct isl_basic_map *bmap;
        unsigned d;
        unsigned nparam;
@@ -44,10 +63,10 @@ static __isl_give isl_map *set_path_length(__isl_take isl_map *map,
        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);
+       dim = isl_map_get_space(map);
+       d = isl_space_dim(dim, isl_dim_in);
+       nparam = isl_space_dim(dim, isl_dim_param);
+       bmap = isl_basic_map_alloc_space(dim, 0, 1, 1);
        if (exactly) {
                k = isl_basic_map_alloc_equality(bmap);
                c = bmap->eq[k];
@@ -96,8 +115,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 +178,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_space(app, isl_map_get_space(map));
+
        test = isl_map_apply_range(isl_map_copy(map), isl_map_copy(app));
        test = isl_map_union(test, isl_map_copy(map));
 
@@ -191,7 +212,7 @@ static int check_exactness(__isl_take isl_map *map, __isl_take isl_map *app,
  * For any element in this relation, the number of steps taken
  * is equal to the difference in the final coordinates.
  */
-static __isl_give isl_map *path_along_steps(__isl_take isl_dim *dim,
+static __isl_give isl_map *path_along_steps(__isl_take isl_space *dim,
        __isl_keep isl_mat *steps)
 {
        int i, j, k;
@@ -203,11 +224,11 @@ static __isl_give isl_map *path_along_steps(__isl_take isl_dim *dim,
        if (!dim || !steps)
                goto error;
 
-       d = isl_dim_size(dim, isl_dim_in);
+       d = isl_space_dim(dim, isl_dim_in);
        n = steps->n_row;
-       nparam = isl_dim_size(dim, isl_dim_param);
+       nparam = isl_space_dim(dim, isl_dim_param);
 
-       path = isl_basic_map_alloc_dim(isl_dim_copy(dim), n, d, n);
+       path = isl_basic_map_alloc_space(isl_space_copy(dim), n, d, n);
 
        for (i = 0; i < n; ++i) {
                k = isl_basic_map_alloc_div(path);
@@ -241,13 +262,13 @@ static __isl_give isl_map *path_along_steps(__isl_take isl_dim *dim,
                isl_int_set_si(path->ineq[k][1 + nparam + 2 * d + i], 1);
        }
 
-       isl_dim_free(dim);
+       isl_space_free(dim);
 
        path = isl_basic_map_simplify(path);
        path = isl_basic_map_finalize(path);
        return isl_map_from_basic_map(path);
 error:
-       isl_dim_free(dim);
+       isl_space_free(dim);
        isl_basic_map_free(path);
        return NULL;
 }
@@ -304,6 +325,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 +339,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 +424,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 +438,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_space(path));
        is_id = isl_basic_map_is_equal(test, id);
        isl_basic_map_free(test);
        isl_basic_map_free(id);
@@ -421,9 +448,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 +468,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 +519,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,8 +546,19 @@ 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,
+static __isl_give isl_map *path_along_delta(__isl_take isl_space *dim,
        __isl_take isl_basic_set *delta)
 {
        isl_basic_map *path = NULL;
@@ -525,13 +569,14 @@ 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;
        n_div = isl_basic_set_dim(delta, isl_dim_div);
        d = isl_basic_set_dim(delta, isl_dim_set);
        nparam = isl_basic_set_dim(delta, isl_dim_param);
-       path = isl_basic_map_alloc_dim(isl_dim_copy(dim), n_div + d + 1,
+       path = isl_basic_map_alloc_space(isl_space_copy(dim), n_div + d + 1,
                        d + 1 + delta->n_eq, delta->n_eq + delta->n_ineq + 1);
        off = 1 + nparam + 2 * (d + 1) + n_div;
 
@@ -556,8 +601,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_space *dim = isl_basic_set_get_space(delta);
+               delta = isl_basic_set_project_out(delta,
+                                                 isl_dim_param, 0, nparam);
+               delta = isl_basic_set_add_dims(delta, isl_dim_param, nparam);
+               delta = isl_basic_set_reset_space(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)
@@ -575,14 +638,13 @@ static __isl_give isl_map *path_along_delta(__isl_take isl_dim *dim,
        isl_basic_set_free(delta);
        path = isl_basic_map_finalize(path);
        if (is_id) {
-               isl_dim_free(dim);
+               isl_space_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);
+       isl_space_free(dim);
        isl_basic_set_free(delta);
        isl_basic_map_free(path);
        return NULL;
@@ -595,7 +657,7 @@ error:
  *
  *     { [x,x_s] -> [y,y_s] : k = y_s - x_s > 0 }
  */
-static __isl_give isl_map *equate_parameter_to_length(__isl_take isl_dim *dim,
+static __isl_give isl_map *equate_parameter_to_length(__isl_take isl_space *dim,
        unsigned param)
 {
        struct isl_basic_map *bmap;
@@ -603,9 +665,9 @@ static __isl_give isl_map *equate_parameter_to_length(__isl_take isl_dim *dim,
        unsigned nparam;
        int k;
 
-       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);
+       d = isl_space_dim(dim, isl_dim_in);
+       nparam = isl_space_dim(dim, isl_dim_param);
+       bmap = isl_basic_map_alloc_space(dim, 0, 1, 1);
        k = isl_basic_map_alloc_equality(bmap);
        if (k < 0)
                goto error;
@@ -684,7 +746,7 @@ static int is_acyclic(__isl_take isl_map *path)
  * Since each of these paths performs an addition, composition is
  * symmetric and we can simply compose all resulting paths in any order.
  */
-static __isl_give isl_map *construct_extended_path(__isl_take isl_dim *dim,
+static __isl_give isl_map *construct_extended_path(__isl_take isl_space *dim,
        __isl_keep isl_map *map, int *project)
 {
        struct isl_mat *steps = NULL;
@@ -694,7 +756,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_space_copy(dim));
 
        steps = isl_mat_alloc(map->ctx, map->n, d);
        if (!steps)
@@ -709,7 +771,7 @@ static __isl_give isl_map *construct_extended_path(__isl_take isl_dim *dim,
                for (j = 0; j < d; ++j) {
                        int fixed;
 
-                       fixed = isl_basic_set_fast_dim_is_fixed(delta, j,
+                       fixed = isl_basic_set_plain_dim_is_fixed(delta, j,
                                                            &steps->row[n][j]);
                        if (fixed < 0) {
                                isl_basic_set_free(delta);
@@ -722,7 +784,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_along_delta(isl_space_copy(dim), delta));
                        path = isl_map_coalesce(path);
                } else {
                        isl_basic_set_free(delta);
@@ -733,7 +795,7 @@ static __isl_give isl_map *construct_extended_path(__isl_take isl_dim *dim,
        if (n > 0) {
                steps->n_row = n;
                path = isl_map_apply_range(path,
-                               path_along_steps(isl_dim_copy(dim), steps));
+                               path_along_steps(isl_space_copy(dim), steps));
        }
 
        if (project && *project) {
@@ -742,11 +804,11 @@ static __isl_give isl_map *construct_extended_path(__isl_take isl_dim *dim,
                        goto error;
        }
 
-       isl_dim_free(dim);
+       isl_space_free(dim);
        isl_mat_free(steps);
        return path;
 error:
-       isl_dim_free(dim);
+       isl_space_free(dim);
        isl_mat_free(steps);
        isl_map_free(path);
        return NULL;
@@ -757,6 +819,9 @@ static int isl_set_overlaps(__isl_keep isl_set *set1, __isl_keep isl_set *set2)
        isl_set *i;
        int no_overlap;
 
+       if (!isl_space_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);
@@ -783,7 +848,7 @@ static int isl_set_overlaps(__isl_keep isl_set *set1, __isl_keep isl_set *set2)
  *                             x in dom R and x + d in ran R and
  *                             \sum_i k_i >= 1 }
  */
-static __isl_give isl_map *construct_component(__isl_take isl_dim *dim,
+static __isl_give isl_map *construct_component(__isl_take isl_space *dim,
        __isl_keep isl_map *map, int *exact, int project)
 {
        struct isl_set *domain = NULL;
@@ -798,19 +863,19 @@ static __isl_give isl_map *construct_component(__isl_take isl_dim *dim,
        if (!isl_set_overlaps(domain, range)) {
                isl_set_free(domain);
                isl_set_free(range);
-               isl_dim_free(dim);
+               isl_space_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,
+       path = construct_extended_path(isl_space_copy(dim), map,
                                        exact && *exact ? &project : NULL);
        app = isl_map_intersect(app, path);
 
@@ -819,11 +884,11 @@ static __isl_give isl_map *construct_component(__isl_take isl_dim *dim,
                                      project)) < 0)
                goto error;
 
-       isl_dim_free(dim);
+       isl_space_free(dim);
        app = set_path_length(app, 0, 1);
        return app;
 error:
-       isl_dim_free(dim);
+       isl_space_free(dim);
        isl_map_free(app);
        return NULL;
 }
@@ -832,7 +897,7 @@ error:
  * the final coordinates.
  */
 static __isl_give isl_map *construct_projected_component(
-       __isl_take isl_dim *dim,
+       __isl_take isl_space *dim,
        __isl_keep isl_map *map, int *exact, int project)
 {
        isl_map *app;
@@ -840,7 +905,7 @@ static __isl_give isl_map *construct_projected_component(
 
        if (!dim)
                return NULL;
-       d = isl_dim_size(dim, isl_dim_in);
+       d = isl_space_dim(dim, isl_dim_in);
 
        app = construct_component(dim, map, exact, project);
        if (project) {
@@ -855,7 +920,7 @@ static __isl_give isl_map *construct_projected_component(
  * with path lengths greater than or equal to zero and with
  * domain and range equal to "dom".
  */
-static __isl_give isl_map *q_closure(__isl_take isl_dim *dim,
+static __isl_give isl_map *q_closure(__isl_take isl_space *dim,
        __isl_take isl_set *dom, __isl_keep isl_basic_map *bmap, int *exact)
 {
        int project = 1;
@@ -863,7 +928,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 +1078,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".
@@ -1023,7 +1095,7 @@ static __isl_give isl_map *compose(__isl_keep isl_map *map, int i,
        int j;
        isl_map *comp;
 
-       comp = isl_map_empty(isl_map_get_dim(map));
+       comp = isl_map_empty(isl_map_get_space(map));
        for (j = 0; j < map->n; ++j) {
                isl_map *map_j;
 
@@ -1031,6 +1103,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])
@@ -1062,7 +1135,7 @@ static __isl_give isl_map *compose(__isl_keep isl_map *map, int i,
  * depending on whether left or right are NULL.
  */
 static __isl_give isl_map *compute_incremental(
-       __isl_take isl_dim *dim, __isl_keep isl_map *map,
+       __isl_take isl_space *dim, __isl_keep isl_map *map,
        int i, __isl_take isl_map *qc, int *left, int *right, int *exact)
 {
        isl_map *map_i;
@@ -1074,7 +1147,7 @@ static __isl_give isl_map *compute_incremental(
        isl_assert(map->ctx, left || right, goto error);
 
        map_i = isl_map_from_basic_map(isl_basic_map_copy(map->p[i]));
-       tc = construct_projected_component(isl_dim_copy(dim), map_i,
+       tc = construct_projected_component(isl_space_copy(dim), map_i,
                                                exact, 1);
        isl_map_free(map_i);
 
@@ -1082,26 +1155,26 @@ static __isl_give isl_map *compute_incremental(
                qc = isl_map_transitive_closure(qc, exact);
 
        if (!*exact) {
-               isl_dim_free(dim);
+               isl_space_free(dim);
                isl_map_free(tc);
                isl_map_free(qc);
-               return isl_map_universe(isl_map_get_dim(map));
+               return isl_map_universe(isl_map_get_space(map));
        }
 
        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_space(tc)));
        if (!right)
                qc = isl_map_apply_range(rtc, qc);
        if (!left)
                qc = isl_map_apply_range(qc, rtc);
        qc = isl_map_union(tc, qc);
 
-       isl_dim_free(dim);
+       isl_space_free(dim);
 
        return qc;
 error:
-       isl_dim_free(dim);
+       isl_space_free(dim);
        isl_map_free(qc);
        return NULL;
 }
@@ -1121,7 +1194,7 @@ error:
  * after computing the integer divisions, is smaller than the number
  * of basic maps in the input map.
  */
-static int incemental_on_entire_domain(__isl_keep isl_dim *dim,
+static int incemental_on_entire_domain(__isl_keep isl_space *dim,
        __isl_keep isl_map *map,
        isl_set **dom, isl_set **ran, int *left, int *right,
        __isl_give isl_map **res)
@@ -1152,7 +1225,7 @@ static int incemental_on_entire_domain(__isl_keep isl_dim *dim,
                                        isl_basic_map_copy(map->p[i])));
                ran[i] = isl_set_from_basic_set(isl_basic_map_range(
                                        isl_basic_map_copy(map->p[i])));
-               qc = q_closure(isl_dim_copy(dim), isl_set_copy(C),
+               qc = q_closure(isl_space_copy(dim), isl_set_copy(C),
                                map->p[i], &exact_i);
                if (!qc)
                        goto error;
@@ -1179,7 +1252,7 @@ static int incemental_on_entire_domain(__isl_keep isl_dim *dim,
                        isl_map_free(qc);
                        continue;
                }
-               *res = compute_incremental(isl_dim_copy(dim), map, i, qc,
+               *res = compute_incremental(isl_space_copy(dim), map, i, qc,
                                left, right, &exact_i);
                if (!*res)
                        goto error;
@@ -1205,7 +1278,7 @@ error:
  * with C either the simple hull of the domain and range of the entire
  * map or the simple hull of domain and range of map_i.
  */
-static __isl_give isl_map *incremental_closure(__isl_take isl_dim *dim,
+static __isl_give isl_map *incremental_closure(__isl_take isl_space *dim,
        __isl_keep isl_map *map, int *exact, int project)
 {
        int i;
@@ -1268,7 +1341,7 @@ static __isl_give isl_map *incremental_closure(__isl_take isl_dim *dim,
                                goto error;
                        continue;
                }
-               qc = q_closure(isl_dim_copy(dim), C, map->p[i], &exact_i);
+               qc = q_closure(isl_space_copy(dim), C, map->p[i], &exact_i);
                if (!qc)
                        goto error;
                if (!exact_i) {
@@ -1293,7 +1366,7 @@ static __isl_give isl_map *incremental_closure(__isl_take isl_dim *dim,
                        isl_map_free(qc);
                        continue;
                }
-               res = compute_incremental(isl_dim_copy(dim), map, i, qc,
+               res = compute_incremental(isl_space_copy(dim), map, i, qc,
                                (comp & LEFT) ? left : NULL,
                                (comp & RIGHT) ? right : NULL, &exact_i);
                if (!res)
@@ -1314,7 +1387,7 @@ static __isl_give isl_map *incremental_closure(__isl_take isl_dim *dim,
        free(right);
 
        if (res) {
-               isl_dim_free(dim);
+               isl_space_free(dim);
                return res;
        }
 
@@ -1330,7 +1403,7 @@ error:
        free(ran);
        free(left);
        free(right);
-       isl_dim_free(dim);
+       isl_space_free(dim);
        return NULL;
 }
 
@@ -1363,9 +1436,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;
        }
@@ -1383,7 +1456,7 @@ error:
 static int add_length(__isl_keep isl_map *map, isl_map ***grid, int n)
 {
        int i, j, k;
-       isl_dim *dim;
+       isl_space *dim;
        isl_basic_map *bstep;
        isl_map *step;
        unsigned nparam;
@@ -1391,13 +1464,13 @@ static int add_length(__isl_keep isl_map *map, isl_map ***grid, int n)
        if (!map)
                return -1;
 
-       dim = isl_map_get_dim(map);
-       nparam = isl_dim_size(dim, isl_dim_param);
-       dim = isl_dim_drop(dim, isl_dim_in, 0, isl_dim_size(dim, isl_dim_in));
-       dim = isl_dim_drop(dim, isl_dim_out, 0, isl_dim_size(dim, isl_dim_out));
-       dim = isl_dim_add(dim, isl_dim_in, 1);
-       dim = isl_dim_add(dim, isl_dim_out, 1);
-       bstep = isl_basic_map_alloc_dim(dim, 0, 1, 0);
+       dim = isl_map_get_space(map);
+       nparam = isl_space_dim(dim, isl_dim_param);
+       dim = isl_space_drop_dims(dim, isl_dim_in, 0, isl_space_dim(dim, isl_dim_in));
+       dim = isl_space_drop_dims(dim, isl_dim_out, 0, isl_space_dim(dim, isl_dim_out));
+       dim = isl_space_add_dims(dim, isl_dim_in, 1);
+       dim = isl_space_add_dims(dim, isl_dim_out, 1);
+       bstep = isl_basic_map_alloc_space(dim, 0, 1, 0);
        k = isl_basic_map_alloc_equality(bstep);
        if (k < 0) {
                isl_basic_map_free(bstep);
@@ -1420,6 +1493,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 +1546,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
@@ -1447,11 +1554,10 @@ static int add_length(__isl_keep isl_map *map, isl_map ***grid, int n)
  * the input relation by the cross product with the unit length relation
  * { [i] -> [i + 1] }.
  */
-static __isl_give isl_map *floyd_warshall_with_groups(__isl_take isl_dim *dim,
+static __isl_give isl_map *floyd_warshall_with_groups(__isl_take isl_space *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;
 
@@ -1471,7 +1577,7 @@ static __isl_give isl_map *floyd_warshall_with_groups(__isl_take isl_dim *dim,
                if (!grid[i])
                        goto error;
                for (j = 0; j < n; ++j)
-                       grid[i][j] = isl_map_empty(isl_map_get_dim(map));
+                       grid[i][j] = isl_map_empty(isl_map_get_space(map));
        }
 
        for (k = 0; k < map->n; ++k) {
@@ -1485,33 +1591,9 @@ 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));
+       app = isl_map_empty(isl_map_get_space(map));
 
        for (i = 0; i < n; ++i) {
                for (j = 0; j < n; ++j)
@@ -1521,7 +1603,7 @@ static __isl_give isl_map *floyd_warshall_with_groups(__isl_take isl_dim *dim,
        free(grid);
 
        free(group);
-       isl_dim_free(dim);
+       isl_space_free(dim);
 
        return app;
 error:
@@ -1535,17 +1617,12 @@ error:
                }
        free(grid);
        free(group);
-       isl_dim_free(dim);
+       isl_space_free(dim);
        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,140 +1635,107 @@ 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;
-
-       if (!map)
-               goto error;
-       if (map->n <= 1)
-               return incremental_closure(dim, map, exact, project);
+       int g;
 
-       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]];
 
-       for (i = 0; i < 2 * map->n; ++i)
-               isl_set_free(set[i]);
-
-       free(set);
+       *n_group = g;
 
-       return floyd_warshall_with_groups(dim, map, exact, project, group, n);
+       return group;
 error:
-       for (i = 0; i < 2 * map->n; ++i)
-               isl_set_free(set[i]);
-       free(set);
+       if (*set) {
+               for (i = 0; i < 2 * n; ++i)
+                       isl_set_free((*set)[i]);
+               free(*set);
+               *set = NULL;
+       }
        free(group);
-       isl_dim_free(dim);
        return NULL;
 }
 
-/* Structure for representing the nodes in the graph being traversed
- * using Tarjan's algorithm.
- * index represents the order in which nodes are visited.
- * min_index is the index of the root of a (sub)component.
- * on_stack indicates whether the node is currently on the stack.
- */
-struct basic_map_sort_node {
-       int index;
-       int min_index;
-       int on_stack;
-};
-/* Structure for representing the graph being traversed
- * using Tarjan's algorithm.
- * len is the number of nodes
- * node is an array of nodes
- * stack contains the nodes on the path from the root to the current node
- * sp is the stack pointer
- * index is the index of the last node visited
- * order contains the elements of the components separated by -1
- * op represents the current position in order
- *
- * check_closed is set if we may have used the fact that
- * a pair of basic maps can be interchanged
+/* 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.
  */
-struct basic_map_sort {
-       int len;
-       struct basic_map_sort_node *node;
-       int *stack;
-       int sp;
-       int index;
-       int *order;
-       int op;
-       int check_closed;
-};
-
-static void basic_map_sort_free(struct basic_map_sort *s)
-{
-       if (!s)
-               return;
-       free(s->node);
-       free(s->stack);
-       free(s->order);
-       free(s);
-}
-
-static struct basic_map_sort *basic_map_sort_alloc(struct isl_ctx *ctx, int len)
+static __isl_give isl_map *floyd_warshall(__isl_take isl_space *dim,
+       __isl_keep isl_map *map, int *exact, int project)
 {
-       struct basic_map_sort *s;
        int i;
+       isl_set **set = NULL;
+       int *group = NULL;
+       int n;
 
-       s = isl_calloc_type(ctx, struct basic_map_sort);
-       if (!s)
-               return NULL;
-       s->len = len;
-       s->node = isl_alloc_array(ctx, struct basic_map_sort_node, len);
-       if (!s->node)
-               goto error;
-       for (i = 0; i < len; ++i)
-               s->node[i].index = -1;
-       s->stack = isl_alloc_array(ctx, int, len);
-       if (!s->stack)
+       if (!map)
                goto error;
-       s->order = isl_alloc_array(ctx, int, 2 * len);
-       if (!s->order)
+       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;
 
-       s->sp = 0;
-       s->index = 0;
-       s->op = 0;
+       for (i = 0; i < 2 * map->n; ++i)
+               isl_set_free(set[i]);
 
-       s->check_closed = 0;
+       free(set);
 
-       return s;
+       return floyd_warshall_with_groups(dim, map, exact, project, group, n);
 error:
-       basic_map_sort_free(s);
+       isl_space_free(dim);
        return NULL;
 }
 
+/* Structure for representing the nodes of the graph of which
+ * strongly connected components are being computed.
+ *
+ * list contains the actual nodes
+ * check_closed is set if we may have used the fact that
+ * a pair of basic maps can be interchanged
+ */
+struct isl_tc_follows_data {
+       isl_basic_map **list;
+       int check_closed;
+};
+
 /* Check whether in the computation of the transitive closure
- * "bmap1" (R_1) should follow (or be part of the same component as)
- * "bmap2" (R_2).
+ * "list[i]" (R_1) should follow (or be part of the same component as)
+ * "list[j]" (R_2).
  *
  * That is check whether
  *
@@ -1707,17 +1751,21 @@ error:
  * *check_closed is set if the subset relation holds while
  * R_1 \circ R_2 is not empty.
  */
-static int basic_map_follows(__isl_keep isl_basic_map *bmap1,
-       __isl_keep isl_basic_map *bmap2, int *check_closed)
+static int basic_map_follows(int i, int j, void *user)
 {
+       struct isl_tc_follows_data *data = user;
        struct isl_map *map12 = NULL;
        struct isl_map *map21 = NULL;
        int subset;
 
+       if (!isl_space_tuple_match(data->list[i]->dim, isl_dim_in,
+                                   data->list[j]->dim, isl_dim_out))
+               return 0;
+
        map21 = isl_map_from_basic_map(
                        isl_basic_map_apply_range(
-                               isl_basic_map_copy(bmap2),
-                               isl_basic_map_copy(bmap1)));
+                               isl_basic_map_copy(data->list[j]),
+                               isl_basic_map_copy(data->list[i])));
        subset = isl_map_is_empty(map21);
        if (subset < 0)
                goto error;
@@ -1726,10 +1774,18 @@ static int basic_map_follows(__isl_keep isl_basic_map *bmap1,
                return 0;
        }
 
+       if (!isl_space_tuple_match(data->list[i]->dim, isl_dim_in,
+                                   data->list[i]->dim, isl_dim_out) ||
+           !isl_space_tuple_match(data->list[j]->dim, isl_dim_in,
+                                   data->list[j]->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),
-                               isl_basic_map_copy(bmap2)));
+                               isl_basic_map_copy(data->list[i]),
+                               isl_basic_map_copy(data->list[j])));
 
        subset = isl_map_is_subset(map21, map12);
 
@@ -1737,7 +1793,7 @@ static int basic_map_follows(__isl_keep isl_basic_map *bmap1,
        isl_map_free(map21);
 
        if (subset)
-               *check_closed = 1;
+               data->check_closed = 1;
 
        return subset < 0 ? -1 : !subset;
 error:
@@ -1745,59 +1801,6 @@ error:
        return -1;
 }
 
-/* Perform Tarjan's algorithm for computing the strongly connected components
- * in the graph with the disjuncts of "map" as vertices and with an
- * edge between any pair of disjuncts such that the first has
- * to be applied after the second.
- */
-static int power_components_tarjan(struct basic_map_sort *s,
-       __isl_keep isl_map *map, int i)
-{
-       int j;
-
-       s->node[i].index = s->index;
-       s->node[i].min_index = s->index;
-       s->node[i].on_stack = 1;
-       s->index++;
-       s->stack[s->sp++] = i;
-
-       for (j = s->len - 1; j >= 0; --j) {
-               int f;
-
-               if (j == i)
-                       continue;
-               if (s->node[j].index >= 0 &&
-                       (!s->node[j].on_stack ||
-                        s->node[j].index > s->node[i].min_index))
-                       continue;
-
-               f = basic_map_follows(map->p[i], map->p[j], &s->check_closed);
-               if (f < 0)
-                       return -1;
-               if (!f)
-                       continue;
-
-               if (s->node[j].index < 0) {
-                       power_components_tarjan(s, map, 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)
-                       s->node[i].min_index = s->node[j].index;
-       }
-
-       if (s->node[i].index != s->node[i].min_index)
-               return 0;
-
-       do {
-               j = s->stack[--s->sp];
-               s->node[j].on_stack = 0;
-               s->order[s->op++] = j;
-       } while (j != i);
-       s->order[s->op++] = -1;
-
-       return 0;
-}
-
 /* 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
@@ -1831,12 +1834,13 @@ static int power_components_tarjan(struct basic_map_sort *s,
  * order, at each join also taking in the union of both arguments
  * to allow for paths that do not go through one of the two arguments.
  */
-static __isl_give isl_map *construct_power_components(__isl_take isl_dim *dim,
+static __isl_give isl_map *construct_power_components(__isl_take isl_space *dim,
        __isl_keep isl_map *map, int *exact, int project)
 {
        int i, n, c;
        struct isl_map *path = NULL;
-       struct basic_map_sort *s = NULL;
+       struct isl_tc_follows_data data;
+       struct isl_tarjan_graph *g = NULL;
        int *orig_exact;
        int local_exact;
 
@@ -1845,39 +1849,37 @@ 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);
-       if (!s)
+       data.list = map->p;
+       data.check_closed = 0;
+       g = isl_tarjan_graph_init(map->ctx, map->n, &basic_map_follows, &data);
+       if (!g)
                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)
+       if (data.check_closed && !exact)
                exact = &local_exact;
 
        c = 0;
        i = 0;
        n = map->n;
        if (project)
-               path = isl_map_empty(isl_map_get_dim(map));
+               path = isl_map_empty(isl_map_get_space(map));
        else
-               path = isl_map_empty(isl_dim_copy(dim));
+               path = isl_map_empty(isl_space_copy(dim));
+       path = anonymize(path);
        while (n) {
                struct isl_map *comp;
                isl_map *path_comp, *path_comb;
-               comp = isl_map_alloc_dim(isl_map_get_dim(map), n, 0);
-               while (s->order[i] != -1) {
+               comp = isl_map_alloc_space(isl_map_get_space(map), n, 0);
+               while (g->order[i] != -1) {
                        comp = isl_map_add_basic_map(comp,
-                                   isl_basic_map_copy(map->p[s->order[i]]));
+                                   isl_basic_map_copy(map->p[g->order[i]]));
                        --n;
                        ++i;
                }
-               path_comp = floyd_warshall(isl_dim_copy(dim),
+               path_comp = floyd_warshall(isl_space_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);
@@ -1887,26 +1889,26 @@ static __isl_give isl_map *construct_power_components(__isl_take isl_dim *dim,
                ++c;
        }
 
-       if (c > 1 && s->check_closed && !*exact) {
+       if (c > 1 && data.check_closed && !*exact) {
                int closed;
 
                closed = isl_map_is_transitively_closed(path);
                if (closed < 0)
                        goto error;
                if (!closed) {
-                       basic_map_sort_free(s);
+                       isl_tarjan_graph_free(g);
                        isl_map_free(path);
                        return floyd_warshall(dim, map, orig_exact, project);
                }
        }
 
-       basic_map_sort_free(s);
-       isl_dim_free(dim);
+       isl_tarjan_graph_free(g);
+       isl_space_free(dim);
 
        return path;
 error:
-       basic_map_sort_free(s);
-       isl_dim_free(dim);
+       isl_tarjan_graph_free(g);
+       isl_space_free(dim);
        isl_map_free(path);
        return NULL;
 }
@@ -1944,22 +1946,22 @@ static __isl_give isl_map *construct_power(__isl_keep isl_map *map,
        int *exact, int project)
 {
        struct isl_map *app = NULL;
-       struct isl_dim *dim = NULL;
+       isl_space *dim = NULL;
        unsigned d;
 
        if (!map)
                return NULL;
 
-       dim = isl_map_get_dim(map);
+       dim = isl_map_get_space(map);
 
-       d = isl_dim_size(dim, isl_dim_in);
-       dim = isl_dim_add(dim, isl_dim_in, 1);
-       dim = isl_dim_add(dim, isl_dim_out, 1);
+       d = isl_space_dim(dim, isl_dim_in);
+       dim = isl_space_add_dims(dim, isl_dim_in, 1);
+       dim = isl_space_add_dims(dim, isl_dim_out, 1);
 
-       app = construct_power_components(isl_dim_copy(dim), map,
+       app = construct_power_components(isl_space_copy(dim), map,
                                        exact, project);
 
-       isl_dim_free(dim);
+       isl_space_free(dim);
 
        return app;
 }
@@ -1998,46 +2000,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 *dim;
+       isl_space *target_dim;
+       isl_space *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_plain_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_space(map);
+       target_dim = isl_space_from_range(isl_space_wrap(target_dim));
+       target_dim = isl_space_add_dims(target_dim, isl_dim_in, 1);
+       target_dim = isl_space_set_dim_name(target_dim, isl_dim_in, 0, "k");
 
        map = map_power(map, exact, 0);
 
-       dim = isl_map_get_dim(map);
+       map = isl_map_add_dims(map, isl_dim_param, 1);
+       dim = isl_map_get_space(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_space(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_space *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_plain_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_space(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
@@ -2049,7 +2110,6 @@ error:
  */
 static int is_eq_stride(__isl_keep isl_basic_set *bset, int i)
 {
-       int k;
        unsigned nparam;
        unsigned d;
        unsigned n_div;
@@ -2125,7 +2185,7 @@ static __isl_give isl_map *box_closure_on_domain(__isl_take isl_map *map,
        unsigned d;
        unsigned nparam;
        unsigned total;
-       isl_dim *dim;
+       isl_space *dim;
        isl_set *delta;
        isl_map *app = NULL;
        isl_basic_set *aff = NULL;
@@ -2140,11 +2200,11 @@ static __isl_give isl_map *box_closure_on_domain(__isl_take isl_map *map,
        aff = isl_set_affine_hull(isl_set_copy(delta));
        if (!aff)
                goto error;
-       dim = isl_map_get_dim(map);
-       d = isl_dim_size(dim, isl_dim_in);
-       nparam = isl_dim_size(dim, isl_dim_param);
-       total = isl_dim_total(dim);
-       bmap = isl_basic_map_alloc_dim(dim,
+       dim = isl_map_get_space(map);
+       d = isl_space_dim(dim, isl_dim_in);
+       nparam = isl_space_dim(dim, isl_dim_param);
+       total = isl_space_dim(dim, isl_dim_all);
+       bmap = isl_basic_map_alloc_space(dim,
                                        aff->n_div + 1, aff->n_div, 2 * d + 1);
        for (i = 0; i < aff->n_div + 1; ++i) {
                k = isl_basic_map_alloc_div(bmap);
@@ -2467,7 +2527,7 @@ static __isl_give isl_map *transitive_closure_omega(__isl_take isl_map *map,
                if (!ok)
                        continue;
 
-               app = isl_map_alloc_dim(isl_map_get_dim(map), map->n - 1, 0);
+               app = isl_map_alloc_space(isl_map_get_space(map), map->n - 1, 0);
 
                for (j = 0; j < map->n; ++j) {
                        if (j == i)
@@ -2507,12 +2567,13 @@ error:
 __isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
        int *exact)
 {
+       isl_space *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 +2587,392 @@ __isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
                return map;
        }
 
+       target_dim = isl_map_get_space(map);
        map = map_power(map, exact, 1);
+       map = isl_map_reset_space(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(ctx, isl_map *, n_group);
+               if (!grid[i])
+                       goto error;
+               for (j = 0; j < n_group; ++j) {
+                       isl_space *dim1, *dim2, *dim;
+                       dim1 = isl_space_reverse(isl_set_get_space(set[i]));
+                       dim2 = isl_set_get_space(set[j]);
+                       dim = isl_space_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_space(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 = NULL;
+       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 = NULL;
+       isl_basic_map **next;
+       isl_union_map *path = NULL;
+       struct isl_tc_follows_data data;
+       struct isl_tarjan_graph *g = 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;
+
+       data.list = list;
+       data.check_closed = 0;
+       g = isl_tarjan_graph_init(ctx, n, &basic_map_follows, &data);
+       if (!g)
+               goto error;
+
+       c = 0;
+       i = 0;
+       l = n;
+       path = isl_union_map_empty(isl_union_map_get_space(umap));
+       while (l) {
+               isl_union_map *comp;
+               isl_union_map *path_comp, *path_comb;
+               comp = isl_union_map_empty(isl_union_map_get_space(umap));
+               while (g->order[i] != -1) {
+                       comp = isl_union_map_add_map(comp,
+                                   isl_map_from_basic_map(
+                                       isl_basic_map_copy(list[g->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 && data.check_closed && !*exact) {
+               int closed;
+
+               closed = isl_union_map_is_transitively_closed(path);
+               if (closed < 0)
+                       goto error;
+               recheck = !closed;
+       }
+
+       isl_tarjan_graph_free(g);
+
+       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:
+       isl_tarjan_graph_free(g);
+       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;
+}
+
+struct isl_union_power {
+       isl_union_map *pow;
+       int *exact;
+};
+
+static int power(__isl_take isl_map *map, void *user)
+{
+       struct isl_union_power *up = user;
+
+       map = isl_map_power(map, up->exact);
+       up->pow = isl_union_map_from_map(map);
+
+       return -1;
+}
+
+/* Construct a map [x] -> [x+1], with parameters prescribed by "dim".
+ */
+static __isl_give isl_union_map *increment(__isl_take isl_space *dim)
+{
+       int k;
+       isl_basic_map *bmap;
+
+       dim = isl_space_add_dims(dim, isl_dim_in, 1);
+       dim = isl_space_add_dims(dim, isl_dim_out, 1);
+       bmap = isl_basic_map_alloc_space(dim, 0, 1, 0);
+       k = isl_basic_map_alloc_equality(bmap);
+       if (k < 0)
+               goto error;
+       isl_seq_clr(bmap->eq[k], isl_basic_map_total_dim(bmap));
+       isl_int_set_si(bmap->eq[k][0], 1);
+       isl_int_set_si(bmap->eq[k][isl_basic_map_offset(bmap, isl_dim_in)], 1);
+       isl_int_set_si(bmap->eq[k][isl_basic_map_offset(bmap, isl_dim_out)], -1);
+       return isl_union_map_from_map(isl_map_from_basic_map(bmap));
+error:
+       isl_basic_map_free(bmap);
+       return NULL;
+}
+
+/* Construct a map [[x]->[y]] -> [y-x], with parameters prescribed by "dim".
+ */
+static __isl_give isl_union_map *deltas_map(__isl_take isl_space *dim)
+{
+       isl_basic_map *bmap;
+
+       dim = isl_space_add_dims(dim, isl_dim_in, 1);
+       dim = isl_space_add_dims(dim, isl_dim_out, 1);
+       bmap = isl_basic_map_universe(dim);
+       bmap = isl_basic_map_deltas_map(bmap);
+
+       return isl_union_map_from_map(isl_map_from_basic_map(bmap));
+}
+
+/* Compute the positive powers of "map", or an overapproximation.
+ * The result maps the exponent to a nested copy of the corresponding power.
+ * If the result is exact, then *exact is set to 1.
+ */
+__isl_give isl_union_map *isl_union_map_power(__isl_take isl_union_map *umap,
+       int *exact)
+{
+       int n;
+       isl_union_map *inc;
+       isl_union_map *dm;
+
+       if (!umap)
+               return NULL;
+       n = isl_union_map_n_map(umap);
+       if (n == 0)
+               return umap;
+       if (n == 1) {
+               struct isl_union_power up = { NULL, exact };
+               isl_union_map_foreach_map(umap, &power, &up);
+               isl_union_map_free(umap);
+               return up.pow;
+       }
+       inc = increment(isl_union_map_get_space(umap));
+       umap = isl_union_map_product(inc, umap);
+       umap = isl_union_map_transitive_closure(umap, exact);
+       umap = isl_union_map_zip(umap);
+       dm = deltas_map(isl_union_map_get_space(umap));
+       umap = isl_union_map_apply_domain(umap, dm);
+       
+       return umap;
+}
+
+#undef TYPE
+#define TYPE isl_map
+#include "isl_power_templ.c"
+
+#undef TYPE
+#define TYPE isl_union_map
+#include "isl_power_templ.c"