isl_basic_set_opt: avoid invalid access on error path
[platform/upstream/isl.git] / isl_transitive_closure.c
index 53930fb..70d15b9 100644 (file)
@@ -1,19 +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_dim_private.h>
-#include <isl_lp.h>
-#include <isl_union_map.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)
 {
@@ -49,7 +53,7 @@ int isl_union_map_is_transitively_closed(__isl_keep isl_union_map *umap)
 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;
@@ -59,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];
@@ -111,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);
@@ -174,7 +178,7 @@ 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));
+       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));
@@ -208,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;
@@ -220,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);
@@ -258,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;
 }
@@ -321,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)
@@ -332,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);
@@ -414,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;
@@ -428,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);
@@ -438,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;
@@ -454,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) {
@@ -503,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 }
  *
@@ -530,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;
@@ -542,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;
 
@@ -573,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)
@@ -592,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;
@@ -612,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;
@@ -620,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;
@@ -701,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;
@@ -711,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)
@@ -726,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);
@@ -739,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);
@@ -750,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) {
@@ -759,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;
@@ -774,7 +819,7 @@ 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))
+       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));
@@ -803,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;
@@ -818,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);
 
@@ -839,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;
 }
@@ -852,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;
@@ -860,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) {
@@ -875,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;
@@ -883,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);
@@ -1033,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".
@@ -1043,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;
 
@@ -1051,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])
@@ -1082,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;
@@ -1094,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);
 
@@ -1102,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;
 }
@@ -1141,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)
@@ -1172,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;
@@ -1199,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;
@@ -1225,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;
@@ -1288,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) {
@@ -1313,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)
@@ -1334,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;
        }
 
@@ -1350,7 +1403,7 @@ error:
        free(ran);
        free(left);
        free(right);
-       isl_dim_free(dim);
+       isl_space_free(dim);
        return NULL;
 }
 
@@ -1403,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;
@@ -1411,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);
@@ -1501,7 +1554,7 @@ static void floyd_warshall_iterate(isl_map ***grid, int n, int *exact)
  * 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;
@@ -1524,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) {
@@ -1540,7 +1593,7 @@ static __isl_give isl_map *floyd_warshall_with_groups(__isl_take isl_dim *dim,
 
        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)
@@ -1550,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:
@@ -1564,7 +1617,7 @@ error:
                }
        free(grid);
        free(group);
-       isl_dim_free(dim);
+       isl_space_free(dim);
        return NULL;
 }
 
@@ -1640,7 +1693,7 @@ error:
  * 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,
+static __isl_give isl_map *floyd_warshall(__isl_take isl_space *dim,
        __isl_keep isl_map *map, int *exact, int project)
 {
        int i;
@@ -1664,91 +1717,25 @@ static __isl_give isl_map *floyd_warshall(__isl_take isl_dim *dim,
 
        return floyd_warshall_with_groups(dim, map, exact, project, group, n);
 error:
-       isl_dim_free(dim);
+       isl_space_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
+/* 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 basic_map_sort {
-       int len;
-       struct basic_map_sort_node *node;
-       int *stack;
-       int sp;
-       int index;
-       int *order;
-       int op;
+struct isl_tc_follows_data {
+       isl_basic_map **list;
        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)
-{
-       struct basic_map_sort *s;
-       int i;
-
-       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)
-               goto error;
-       s->order = isl_alloc_array(ctx, int, 2 * len);
-       if (!s->order)
-               goto error;
-
-       s->sp = 0;
-       s->index = 0;
-       s->op = 0;
-
-       s->check_closed = 0;
-
-       return s;
-error:
-       basic_map_sort_free(s);
-       return NULL;
-}
-
 /* 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
  *
@@ -1764,20 +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_dim_tuple_match(bmap1->dim, isl_dim_in, bmap2->dim, isl_dim_out))
+       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;
@@ -1786,16 +1774,18 @@ 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)) {
+       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);
 
@@ -1803,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:
@@ -1811,84 +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_basic_map **list, 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(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, 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)
-                       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;
-}
-
-/* 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
@@ -1922,12 +1834,13 @@ error:
  * 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;
 
@@ -1936,33 +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_init(map->ctx, map->n, map->p);
-       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;
 
        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);
@@ -1972,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;
 }
@@ -2029,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;
 }
@@ -2083,50 +2000,57 @@ 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_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");
 
-       target_dim = isl_map_get_dim(map);
        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_dim(map, target_dim);
+       map = isl_map_reset_space(map, target_dim);
 
        return map;
-error:
-       isl_map_free(map);
-       return NULL;
 }
 
 /* Compute a relation that maps each element in the range of the input
@@ -2141,7 +2065,7 @@ error:
 __isl_give isl_map *isl_map_reaching_path_lengths(__isl_take isl_map *map,
        int *exact)
 {
-       isl_dim *dim;
+       isl_space *dim;
        isl_map *diff;
        unsigned d;
        unsigned param;
@@ -2155,18 +2079,18 @@ __isl_give isl_map *isl_map_reaching_path_lengths(__isl_take isl_map *map,
        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)) {
                if (exact)
                        *exact = 1;
                map = isl_map_project_out(map, isl_dim_out, 0, d);
-               map = isl_map_add(map, isl_dim_out, 1);
+               map = isl_map_add_dims(map, isl_dim_out, 1);
                return map;
        }
 
        map = map_power(map, exact, 0);
 
-       map = isl_map_add(map, isl_dim_param, 1);
-       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, 0, d + 1);
@@ -2186,7 +2110,6 @@ __isl_give isl_map *isl_map_reaching_path_lengths(__isl_take isl_map *map,
  */
 static int is_eq_stride(__isl_keep isl_basic_set *bset, int i)
 {
-       int k;
        unsigned nparam;
        unsigned d;
        unsigned n_div;
@@ -2262,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;
@@ -2277,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);
@@ -2604,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)
@@ -2644,13 +2567,13 @@ error:
 __isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
        int *exact)
 {
-       isl_dim *target_dim;
+       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);
@@ -2664,9 +2587,9 @@ __isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
                return map;
        }
 
-       target_dim = isl_map_get_dim(map);
+       target_dim = isl_map_get_space(map);
        map = map_power(map, exact, 1);
-       map = isl_map_reset_dim(map, target_dim);
+       map = isl_map_reset_space(map, target_dim);
 
        return map;
 error:
@@ -2730,14 +2653,14 @@ static __isl_give isl_union_map *union_floyd_warshall_on_list(isl_ctx *ctx,
        if (!grid)
                goto error;
        for (i = 0; i < n_group; ++i) {
-               grid[i] = isl_calloc_array(map->ctx, isl_map *, n_group);
+               grid[i] = isl_calloc_array(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);
+                       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);
                }
        }
@@ -2752,7 +2675,7 @@ static __isl_give isl_union_map *union_floyd_warshall_on_list(isl_ctx *ctx,
        
        floyd_warshall_iterate(grid, n_group, exact);
 
-       app = isl_union_map_empty(isl_map_get_dim(grid[0][0]));
+       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)
@@ -2797,7 +2720,7 @@ static __isl_give isl_union_map *union_floyd_warshall(
 {
        int i, n;
        isl_ctx *ctx;
-       isl_basic_map **list;
+       isl_basic_map **list = NULL;
        isl_basic_map **next;
        isl_union_map *res;
 
@@ -2845,10 +2768,11 @@ static __isl_give isl_union_map *union_components(
        int i;
        int n;
        isl_ctx *ctx;
-       isl_basic_map **list;
+       isl_basic_map **list = NULL;
        isl_basic_map **next;
        isl_union_map *path = NULL;
-       struct basic_map_sort *s = NULL;
+       struct isl_tc_follows_data data;
+       struct isl_tarjan_graph *g = NULL;
        int c, l;
        int recheck = 0;
 
@@ -2868,22 +2792,24 @@ static __isl_give isl_union_map *union_components(
        if (isl_union_map_foreach_map(umap, collect_basic_map, &next) < 0)
                goto error;
 
-       s = basic_map_sort_init(ctx, n, list);
-       if (!s)
+       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_dim(umap));
+       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_dim(umap));
-               while (s->order[i] != -1) {
+               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[s->order[i]])));
+                                       isl_basic_map_copy(list[g->order[i]])));
                        --l;
                        ++i;
                }
@@ -2896,7 +2822,7 @@ static __isl_give isl_union_map *union_components(
                ++c;
        }
 
-       if (c > 1 && s->check_closed && !*exact) {
+       if (c > 1 && data.check_closed && !*exact) {
                int closed;
 
                closed = isl_union_map_is_transitively_closed(path);
@@ -2905,7 +2831,7 @@ static __isl_give isl_union_map *union_components(
                recheck = !closed;
        }
 
-       basic_map_sort_free(s);
+       isl_tarjan_graph_free(g);
 
        for (i = 0; i < n; ++i)
                isl_basic_map_free(list[i]);
@@ -2920,7 +2846,7 @@ static __isl_give isl_union_map *union_components(
 
        return path;
 error:
-       basic_map_sort_free(s);
+       isl_tarjan_graph_free(g);
        if (list) {
                for (i = 0; i < n; ++i)
                        isl_basic_map_free(list[i]);
@@ -2958,3 +2884,95 @@ 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"