X-Git-Url: http://review.tizen.org/git/?a=blobdiff_plain;f=isl_transitive_closure.c;h=70d15b91ee3f4fa76c9f6a8711b5a1468e2902dd;hb=de51a9bc4da5dd3f1f9f57c2362da6f9752c44e0;hp=4767040e6a4f7dceaf2c887a70672b462d038781;hpb=9146e946bd4e6d0b2515775e297a66115a320523;p=platform%2Fupstream%2Fisl.git

diff --git a/isl_transitive_closure.c b/isl_transitive_closure.c
index 4767040..70d15b9 100644
--- a/isl_transitive_closure.c
+++ b/isl_transitive_closure.c
@@ -1,17 +1,48 @@
 /*
  * 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)
+{
+	isl_map *map2;
+	int closed;
+
+	map2 = isl_map_apply_range(isl_map_copy(map), isl_map_copy(map));
+	closed = isl_map_is_subset(map2, map);
+	isl_map_free(map2);
+
+	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
@@ -22,7 +53,7 @@
 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;
@@ -32,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];
@@ -84,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);
@@ -147,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));
 
@@ -158,10 +191,6 @@ static int check_exactness(__isl_take isl_map *map, __isl_take isl_map *app,
 	isl_map_free(map);
 
 	return exact;
-error:
-	isl_map_free(app);
-	isl_map_free(map);
-	return -1;
 }
 
 /*
@@ -183,7 +212,7 @@ error:
  * 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;
@@ -195,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);
@@ -233,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;
 }
@@ -296,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)
@@ -307,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);
@@ -389,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;
@@ -403,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);
@@ -413,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;
@@ -429,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) {
@@ -478,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 >= 0 and B p + b >= 0 and
  *				C x + C'p + c >= 0 and
  *				D x + D'p + d >= 0 }
  *
@@ -505,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;
@@ -517,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;
 
@@ -548,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)
@@ -567,27 +638,26 @@ 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;
 }
 
-/* Given a dimenion specification Z^{n+1} -> Z^{n+1} and a parameter "param",
+/* Given a dimension specification Z^{n+1} -> Z^{n+1} and a parameter "param",
  * construct a map that equates the parameter to the difference
  * in the final coordinates and imposes that this difference is positive.
  * That is, construct
  *
  *	{ [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;
@@ -595,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;
@@ -676,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;
@@ -686,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)
@@ -701,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);
@@ -714,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);
@@ -725,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) {
@@ -734,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;
@@ -749,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);
@@ -775,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;
@@ -790,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);
 
@@ -811,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;
 }
@@ -824,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;
@@ -832,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) {
@@ -847,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;
@@ -855,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);
@@ -1005,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".
@@ -1015,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;
 
@@ -1023,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])
@@ -1054,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;
@@ -1066,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);
 
@@ -1074,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;
 }
@@ -1113,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)
@@ -1144,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;
@@ -1171,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;
@@ -1197,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;
@@ -1260,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) {
@@ -1285,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)
@@ -1306,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;
 	}
 
@@ -1322,7 +1403,7 @@ error:
 	free(ran);
 	free(left);
 	free(right);
-	isl_dim_free(dim);
+	isl_space_free(dim);
 	return NULL;
 }
 
@@ -1355,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;
 	}
@@ -1369,6 +1450,91 @@ error:
 	return -1;
 }
 
+/* Replace each entry in the n by n grid of maps by the cross product
+ * with the relation { [i] -> [i + 1] }.
+ */
+static int add_length(__isl_keep isl_map *map, isl_map ***grid, int n)
+{
+	int i, j, k;
+	isl_space *dim;
+	isl_basic_map *bstep;
+	isl_map *step;
+	unsigned nparam;
+
+	if (!map)
+		return -1;
+
+	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);
+		return -1;
+	}
+	isl_seq_clr(bstep->eq[k], 1 + isl_basic_map_total_dim(bstep));
+	isl_int_set_si(bstep->eq[k][0], 1);
+	isl_int_set_si(bstep->eq[k][1 + nparam], 1);
+	isl_int_set_si(bstep->eq[k][1 + nparam + 1], -1);
+	bstep = isl_basic_map_finalize(bstep);
+	step = isl_map_from_basic_map(bstep);
+
+	for (i = 0; i < n; ++i)
+		for (j = 0; j < n; ++j)
+			grid[i][j] = isl_map_product(grid[i][j],
+						     isl_map_copy(step));
+
+	isl_map_free(step);
+
+	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.
@@ -1380,19 +1546,18 @@ error:
  * 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
+ * coordinates.  We therefore apply the nested transitive closures
+ * to relations that include these lengths.  In particular, we replace
+ * 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;
 
@@ -1412,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) {
@@ -1423,33 +1588,12 @@ static __isl_give isl_map *floyd_warshall_with_groups(__isl_take isl_dim *dim,
 					isl_basic_map_copy(map->p[k])));
 	}
 
-	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;
+	if (!project && add_length(map, grid, n) < 0)
+		goto error;
 
-		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)
@@ -1459,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:
@@ -1473,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 can only apply this algorithm
- * if we want to compute the transitive closure, i.e., when "project"
- * is set.  If we want to compute the power, we need to keep track
- * of the lengths and the recursive calls inside the Floyd-Warshall
- * would result in non-linear lengths.
+/* 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.
@@ -1496,134 +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;
+	int g;
 
-	if (!map)
-		goto error;
-	if (!project || map->n <= 1)
-		return incremental_closure(dim, map, exact, project);
+	*set = isl_calloc_array(ctx, isl_set *, 2 * n);
+	group = isl_alloc_array(ctx, int, 2 * n);
 
-	set = isl_calloc_array(map->ctx, isl_set *, 2 * map->n);
-	group = isl_alloc_array(map->ctx, int, 2 * map->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]);
+	*n_group = g;
 
-	free(set);
-
-	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 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;
-};
-
-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]);
 
-	return s;
+	free(set);
+
+	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
  *
@@ -1635,18 +1747,25 @@ error:
  *
  * If so, then there is no reason for R_1 to immediately follow R_2
  * in any path.
+ *
+ * *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)
+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;
@@ -1655,75 +1774,33 @@ 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);
 
 	isl_map_free(map12);
 	isl_map_free(map21);
 
+	if (subset)
+		data->check_closed = 1;
+
 	return subset < 0 ? -1 : !subset;
 error:
 	isl_map_free(map21);
 	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]);
-		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
@@ -1757,61 +1834,82 @@ 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;
+	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;
 
 	if (!map)
 		goto error;
 	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 (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);
 		path = isl_map_union(path, path_comb);
 		isl_map_free(comp);
 		++i;
+		++c;
+	}
+
+	if (c > 1 && data.check_closed && !*exact) {
+		int closed;
+
+		closed = isl_map_is_transitively_closed(path);
+		if (closed < 0)
+			goto error;
+		if (!closed) {
+			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;
 }
 
@@ -1838,56 +1936,45 @@ error:
  * if "project" is set.
  *
  * If "project" is not set, then
- * we first construct an extended mapping with an extra coordinate
+ * we construct an extended mapping with an extra coordinate
  * that indicates the number of steps taken.  In particular,
  * the difference in the last coordinate is equal to the number
  * of steps taken to move from a domain element to the corresponding
  * image element(s).
- * In the final step, this difference is equated to the parameter "param"
- * and made positive.  The extra coordinates are subsequently projected out.
  */
 static __isl_give isl_map *construct_power(__isl_keep isl_map *map,
-	unsigned param, int *exact, int project)
+	int *exact, int project)
 {
 	struct isl_map *app = NULL;
-	struct isl_map *diff;
-	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);
 
-	if (project) {
-		isl_dim_free(dim);
-	} else {
-		diff = equate_parameter_to_length(dim, param);
-		app = isl_map_intersect(app, diff);
-		app = isl_map_project_out(app, isl_dim_in, d, 1);
-		app = isl_map_project_out(app, isl_dim_out, d, 1);
-	}
+	isl_space_free(dim);
 
 	return app;
 }
 
 /* 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.
+ * If the result is exact, then *exact is set to 1.
  *
  * If project is set, then we are actually interested in the transitive
  * closure, so we can use a more relaxed exactness check.
  * The lengths of the paths are also projected out instead of being
- * equated to "param" (which is then ignored in this case).
+ * encoded as the difference between an extra pair of final coordinates.
  */
-static __isl_give isl_map *map_power(__isl_take isl_map *map, unsigned param,
+static __isl_give isl_map *map_power(__isl_take isl_map *map,
 	int *exact, int project)
 {
 	struct isl_map *app = NULL;
@@ -1895,21 +1982,14 @@ static __isl_give isl_map *map_power(__isl_take isl_map *map, unsigned param,
 	if (exact)
 		*exact = 1;
 
-	map = isl_map_compute_divs(map);
-	map = isl_map_coalesce(map);
 	if (!map)
 		return NULL;
 
-	if (isl_map_fast_is_empty(map))
-		return map;
-
-	isl_assert(map->ctx, project || param < isl_map_dim(map, isl_dim_param),
-		goto error);
 	isl_assert(map->ctx,
 		isl_map_dim(map, isl_dim_in) == isl_map_dim(map, isl_dim_out),
 		goto error);
 
-	app = construct_power(map, param, exact, project);
+	app = construct_power(map, exact, project);
 
 	isl_map_free(map);
 	return app;
@@ -1920,13 +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 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,
+__isl_give isl_map *isl_map_power(__isl_take isl_map *map, int *exact)
+{
+	isl_space *target_dim;
+	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)) {
+		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);
+
+	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)
 {
-	return map_power(map, param, exact, 0);
+	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;
 }
 
 /* Check whether equality i of bset is a pure stride constraint
@@ -1938,7 +2110,6 @@ __isl_give isl_map *isl_map_power(__isl_take isl_map *map, unsigned param,
  */
 static int is_eq_stride(__isl_keep isl_basic_set *bset, int i)
 {
-	int k;
 	unsigned nparam;
 	unsigned d;
 	unsigned n_div;
@@ -2014,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;
@@ -2029,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);
@@ -2356,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)
@@ -2396,19 +2567,412 @@ error:
 __isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
 	int *exact)
 {
-	unsigned param;
+	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);
 
-	param = isl_map_dim(map, isl_dim_param);
-	map = map_power(map, param, exact, 1);
+	map = isl_map_compute_divs(map);
+	map = isl_map_coalesce(map);
+	closed = isl_map_is_transitively_closed(map);
+	if (closed < 0)
+		goto error;
+	if (closed) {
+		if (exact)
+			*exact = 1;
+		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"