X-Git-Url: http://review.tizen.org/git/?a=blobdiff_plain;f=isl_transitive_closure.c;h=ba7153be1ec236e6da6d235cd8aad98ab3d16b90;hb=e06f25819944b4c591fdfdce6b346ac2a7936396;hp=d989f4542ba4e8b8e15711b76c2ba92097832559;hpb=7d49acbe46e23d2a27a4a87195c23a944edb01b2;p=platform%2Fupstream%2Fisl.git

diff --git a/isl_transitive_closure.c b/isl_transitive_closure.c
index d989f45..ba7153b 100644
--- a/isl_transitive_closure.c
+++ b/isl_transitive_closure.c
@@ -12,6 +12,157 @@
 #include "isl_map_private.h"
 #include "isl_seq.h"
  
+/* Given a map that represents a path with the length of the path
+ * encoded as the difference between the last output coordindate
+ * and the last input coordinate, set this length to either
+ * exactly "length" (if "exactly" is set) or at least "length"
+ * (if "exactly" is not set).
+ */
+static __isl_give isl_map *set_path_length(__isl_take isl_map *map,
+	int exactly, int length)
+{
+	struct isl_dim *dim;
+	struct isl_basic_map *bmap;
+	unsigned d;
+	unsigned nparam;
+	int k;
+	isl_int *c;
+
+	if (!map)
+		return NULL;
+
+	dim = isl_map_get_dim(map);
+	d = isl_dim_size(dim, isl_dim_in);
+	nparam = isl_dim_size(dim, isl_dim_param);
+	bmap = isl_basic_map_alloc_dim(dim, 0, 1, 1);
+	if (exactly) {
+		k = isl_basic_map_alloc_equality(bmap);
+		c = bmap->eq[k];
+	} else {
+		k = isl_basic_map_alloc_inequality(bmap);
+		c = bmap->ineq[k];
+	}
+	if (k < 0)
+		goto error;
+	isl_seq_clr(c, 1 + isl_basic_map_total_dim(bmap));
+	isl_int_set_si(c[0], -length);
+	isl_int_set_si(c[1 + nparam + d - 1], -1);
+	isl_int_set_si(c[1 + nparam + d + d - 1], 1);
+
+	bmap = isl_basic_map_finalize(bmap);
+	map = isl_map_intersect(map, isl_map_from_basic_map(bmap));
+
+	return map;
+error:
+	isl_basic_map_free(bmap);
+	isl_map_free(map);
+	return NULL;
+}
+
+/* Check whether the overapproximation of the power of "map" is exactly
+ * the power of "map".  Let R be "map" and A_k the overapproximation.
+ * The approximation is exact if
+ *
+ *	A_1 = R
+ *	A_k = A_{k-1} \circ R			k >= 2
+ *
+ * Since A_k is known to be an overapproximation, we only need to check
+ *
+ *	A_1 \subset R
+ *	A_k \subset A_{k-1} \circ R		k >= 2
+ *
+ * In practice, "app" has an extra input and output coordinate
+ * to encode the length of the path.  So, we first need to add
+ * this coordinate to "map" and set the length of the path to
+ * one.
+ */
+static int check_power_exactness(__isl_take isl_map *map,
+	__isl_take isl_map *app)
+{
+	int exact;
+	isl_map *app_1;
+	isl_map *app_2;
+
+	map = isl_map_add(map, isl_dim_in, 1);
+	map = isl_map_add(map, isl_dim_out, 1);
+	map = set_path_length(map, 1, 1);
+
+	app_1 = set_path_length(isl_map_copy(app), 1, 1);
+
+	exact = isl_map_is_subset(app_1, map);
+	isl_map_free(app_1);
+
+	if (!exact || exact < 0) {
+		isl_map_free(app);
+		isl_map_free(map);
+		return exact;
+	}
+
+	app_1 = set_path_length(isl_map_copy(app), 0, 1);
+	app_2 = set_path_length(app, 0, 2);
+	app_1 = isl_map_apply_range(map, app_1);
+
+	exact = isl_map_is_subset(app_2, app_1);
+
+	isl_map_free(app_1);
+	isl_map_free(app_2);
+
+	return exact;
+}
+
+/* Check whether the overapproximation of the power of "map" is exactly
+ * the power of "map", possibly after projecting out the power (if "project"
+ * is set).
+ *
+ * If "project" is set and if "steps" can only result in acyclic paths,
+ * then we check
+ *
+ *	A = R \cup (A \circ R)
+ *
+ * where A is the overapproximation with the power projected out, i.e.,
+ * an overapproximation of the transitive closure.
+ * More specifically, since A is known to be an overapproximation, we check
+ *
+ *	A \subset R \cup (A \circ R)
+ *
+ * Otherwise, we check if the power is exact.
+ *
+ * Note that "app" has an extra input and output coordinate to encode
+ * the length of the part.  If we are only interested in the transitive
+ * closure, then we can simply project out these coordinates first.
+ */
+static int check_exactness(__isl_take isl_map *map, __isl_take isl_map *app,
+	int project)
+{
+	isl_map *test;
+	int exact;
+	unsigned d;
+
+	if (!project)
+		return check_power_exactness(map, app);
+
+	d = isl_map_dim(map, isl_dim_in);
+	app = set_path_length(app, 0, 1);
+	app = isl_map_project_out(app, isl_dim_in, d, 1);
+	app = isl_map_project_out(app, isl_dim_out, d, 1);
+
+	test = isl_map_apply_range(isl_map_copy(map), isl_map_copy(app));
+	test = isl_map_union(test, isl_map_copy(map));
+
+	exact = isl_map_is_subset(app, test);
+
+	isl_map_free(app);
+	isl_map_free(test);
+
+	isl_map_free(map);
+
+	return exact;
+error:
+	isl_map_free(app);
+	isl_map_free(map);
+	return -1;
+}
+
 /*
  * The transitive closure implementation is based on the paper
  * "Computing the Transitive Closure of a Union of Affine Integer
@@ -95,16 +246,22 @@ error:
 #define IMPURE		0
 #define PURE_PARAM	1
 #define PURE_VAR	2
+#define MIXED		3
 
 /* Return PURE_PARAM if only the coefficients of the parameters are non-zero.
  * Return PURE_VAR if only the coefficients of the set variables are non-zero.
+ * Return MIXED if only the coefficients of the parameters and the set
+ * 	variables are non-zero and if moreover the parametric constant
+ * 	can never attain positive values.
  * Return IMPURE otherwise.
  */
-static int purity(__isl_keep isl_basic_set *bset, isl_int *c)
+static int purity(__isl_keep isl_basic_set *bset, isl_int *c, int eq)
 {
 	unsigned d;
 	unsigned n_div;
 	unsigned nparam;
+	int k;
+	int empty;
 
 	n_div = isl_basic_set_dim(bset, isl_dim_div);
 	d = isl_basic_set_dim(bset, isl_dim_set);
@@ -116,7 +273,25 @@ static int purity(__isl_keep isl_basic_set *bset, isl_int *c)
 		return PURE_VAR;
 	if (isl_seq_first_non_zero(c + 1 + nparam, d) == -1)
 		return PURE_PARAM;
-	return IMPURE;
+	if (eq)
+		return IMPURE;
+
+	bset = isl_basic_set_copy(bset);
+	bset = isl_basic_set_cow(bset);
+	bset = isl_basic_set_extend_constraints(bset, 0, 1);
+	k = isl_basic_set_alloc_inequality(bset);
+	if (k < 0)
+		goto error;
+	isl_seq_clr(bset->ineq[k], 1 + isl_basic_set_total_dim(bset));
+	isl_seq_cpy(bset->ineq[k], c, 1 + nparam);
+	isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
+	empty = isl_basic_set_is_empty(bset);
+	isl_basic_set_free(bset);
+
+	return empty < 0 ? -1 : empty ? MIXED : IMPURE;
+error:
+	isl_basic_set_free(bset);
+	return -1;
 }
 
 /* Given a set of offsets "delta", construct a relation of the
@@ -136,12 +311,16 @@ static int purity(__isl_keep isl_basic_set *bset, isl_int *c)
  * In particular, let delta be defined as
  *
  *	\delta = [p] -> { [x] : A x + a >= and B p + b >= 0 and
- *				C x + C'p + c >= 0 }
+ *				C x + C'p + c >= 0 and
+ *				D x + D'p + d >= 0 }
  *
- * then the relation is constructed as
+ * where the constraints C x + C'p + c >= 0 are such that the parametric
+ * constant term of each constraint j, "C_j x + C'_j p + c_j",
+ * can never attain positive values, then the relation is constructed as
  *
  *	{ [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and
- *			A f + k a >= 0 and B p + b >= 0 and k >= 1 }
+ *			A f + k a >= 0 and B p + b >= 0 and
+ *			C f + C'p + c >= 0 and k >= 1 }
  *	union { [x] -> [x] }
  *
  * Existentially quantified variables in \delta are currently ignored.
@@ -184,7 +363,9 @@ static __isl_give isl_map *path_along_delta(__isl_take isl_dim *dim,
 	}
 
 	for (i = 0; i < delta->n_eq; ++i) {
-		int p = purity(delta, delta->eq[i]);
+		int p = purity(delta, delta->eq[i], 1);
+		if (p < 0)
+			goto error;
 		if (p == IMPURE)
 			continue;
 		k = isl_basic_map_alloc_equality(path);
@@ -200,7 +381,9 @@ static __isl_give isl_map *path_along_delta(__isl_take isl_dim *dim,
 	}
 
 	for (i = 0; i < delta->n_ineq; ++i) {
-		int p = purity(delta, delta->ineq[i]);
+		int p = purity(delta, delta->ineq[i], 0);
+		if (p < 0)
+			goto error;
 		if (p == IMPURE)
 			continue;
 		k = isl_basic_map_alloc_inequality(path);
@@ -211,8 +394,13 @@ static __isl_give isl_map *path_along_delta(__isl_take isl_dim *dim,
 			isl_seq_cpy(path->ineq[k] + off,
 				    delta->ineq[i] + 1 + nparam, d);
 			isl_int_set(path->ineq[k][off + d], delta->ineq[i][0]);
-		} else
+		} else if (p == PURE_PARAM) {
 			isl_seq_cpy(path->ineq[k], delta->ineq[i], 1 + nparam);
+		} else {
+			isl_seq_cpy(path->ineq[k] + off,
+				    delta->ineq[i] + 1 + nparam, d);
+			isl_seq_cpy(path->ineq[k], delta->ineq[i], 1 + nparam);
+		}
 	}
 
 	k = isl_basic_map_alloc_inequality(path);
@@ -368,6 +556,7 @@ static __isl_give isl_map *construct_extended_path(__isl_take isl_dim *dim,
 		if (j < d) {
 			path = isl_map_apply_range(path,
 				path_along_delta(isl_dim_copy(dim), delta));
+			path = isl_map_coalesce(path);
 		} else {
 			isl_basic_set_free(delta);
 			++n;
@@ -398,7 +587,8 @@ error:
 
 /* 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
+ * construct a map that is the union of the identity map and
+ * an overapproximation of the map
  * that takes an element from the dom R \times Z to an
  * element from ran R \times Z, such that the first n coordinates of the
  * difference between them is a sum of differences between images
@@ -412,13 +602,15 @@ error:
  *
  *	{ (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
  *				d = (\sum_i k_i \delta_i, \sum_i k_i) and
- *				x in dom R and x + d in ran R}
+ *				x in dom R and x + d in ran R } union
+ *	{ (x) -> (x) }
  */
-static __isl_give isl_map *construct_extended_power(__isl_take isl_dim *dim,
-	__isl_keep isl_map *map, int *project)
+static __isl_give isl_map *construct_component(__isl_take isl_dim *dim,
+	__isl_keep isl_map *map, int *exact, int project)
 {
 	struct isl_set *domain = NULL;
 	struct isl_set *range = NULL;
+	struct isl_set *overlap;
 	struct isl_map *app = NULL;
 	struct isl_map *path = NULL;
 
@@ -426,207 +618,345 @@ static __isl_give isl_map *construct_extended_power(__isl_take isl_dim *dim,
 	domain = isl_set_coalesce(domain);
 	range = isl_map_range(isl_map_copy(map));
 	range = isl_set_coalesce(range);
+	overlap = isl_set_intersect(isl_set_copy(domain), isl_set_copy(range));
+	if (isl_set_is_empty(overlap) == 1) {
+		isl_set_free(domain);
+		isl_set_free(range);
+		isl_set_free(overlap);
+		isl_dim_free(dim);
+
+		map = isl_map_copy(map);
+		map = isl_map_add(map, isl_dim_in, 1);
+		map = isl_map_add(map, isl_dim_out, 1);
+		map = set_path_length(map, 1, 1);
+		return map;
+	}
+	isl_set_free(overlap);
 	app = isl_map_from_domain_and_range(domain, range);
 	app = isl_map_add(app, isl_dim_in, 1);
 	app = isl_map_add(app, isl_dim_out, 1);
 
-	path = construct_extended_path(dim, map, project);
+	path = construct_extended_path(isl_dim_copy(dim), map,
+					exact && *exact ? &project : NULL);
 	app = isl_map_intersect(app, path);
 
-	return app;
+	if (exact && *exact &&
+	    (*exact = check_exactness(isl_map_copy(map), isl_map_copy(app),
+				      project)) < 0)
+		goto error;
+
+	return isl_map_union(app, isl_map_identity(isl_dim_domain(dim)));
+error:
+	isl_dim_free(dim);
+	isl_map_free(app);
+	return NULL;
 }
 
-/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D,
- * construct a map that is an overapproximation of the map
- * that takes an element from the space D to another
- * element from the same space, such that the difference between
- * them is a strictly positive sum of differences between images
- * and pre-images in one of the R_i.
- * The number of differences in the sum is equated to parameter "param".
- * That is, let
+/* 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
+ */
+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)
+{
+	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;
+
+	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).
  *
- *	\Delta_i = { y - x | (x, y) in R_i }
+ * That is check whether
  *
- * then the constructed map is an overapproximation of
+ *	R_1 \circ R_2
  *
- *	{ (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
- *				d = \sum_i k_i \delta_i and k = \sum_i k_i > 0 }
+ * is a subset of
  *
- * We first 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.
+ *	R_2 \circ R_1
+ *
+ * If so, then there is no reason for R_1 to immediately follow R_2
+ * in any path.
  */
-static __isl_give isl_map *construct_power(__isl_keep isl_map *map,
-	unsigned param, int *project)
+static int basic_map_follows(__isl_keep isl_basic_map *bmap1,
+	__isl_keep isl_basic_map *bmap2)
 {
-	struct isl_map *app = NULL;
-	struct isl_map *diff;
-	struct isl_dim *dim = NULL;
-	unsigned d;
+	struct isl_map *map12 = NULL;
+	struct isl_map *map21 = NULL;
+	int subset;
+
+	map21 = isl_map_from_basic_map(
+			isl_basic_map_apply_range(
+				isl_basic_map_copy(bmap2),
+				isl_basic_map_copy(bmap1)));
+	subset = isl_map_is_empty(map21);
+	if (subset < 0)
+		goto error;
+	if (subset) {
+		isl_map_free(map21);
+		return 0;
+	}
 
-	if (!map)
-		return NULL;
+	map12 = isl_map_from_basic_map(
+			isl_basic_map_apply_range(
+				isl_basic_map_copy(bmap1),
+				isl_basic_map_copy(bmap2)));
 
-	dim = isl_map_get_dim(map);
+	subset = isl_map_is_subset(map21, map12);
 
-	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);
-
-	app = construct_extended_power(isl_dim_copy(dim), map, project);
+	isl_map_free(map12);
+	isl_map_free(map21);
 
-	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);
-
-	return app;
+	return subset < 0 ? -1 : !subset;
+error:
+	isl_map_free(map21);
+	return -1;
 }
 
-/* Shift variable at position "pos" up by one.
- * That is, replace the corresponding variable v by v - 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 __isl_give isl_basic_map *basic_map_shift_pos(
-	__isl_take isl_basic_map *bmap, unsigned pos)
+static int power_components_tarjan(struct basic_map_sort *s,
+	__isl_keep isl_map *map, int i)
 {
-	int i;
+	int j;
 
-	bmap = isl_basic_map_cow(bmap);
-	if (!bmap)
-		return NULL;
+	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 (i = 0; i < bmap->n_eq; ++i)
-		isl_int_sub(bmap->eq[i][0], bmap->eq[i][0], bmap->eq[i][pos]);
+	for (j = s->len - 1; j >= 0; --j) {
+		int f;
 
-	for (i = 0; i < bmap->n_ineq; ++i)
-		isl_int_sub(bmap->ineq[i][0],
-				bmap->ineq[i][0], bmap->ineq[i][pos]);
+		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;
 
-	for (i = 0; i < bmap->n_div; ++i) {
-		if (isl_int_is_zero(bmap->div[i][0]))
+		f = basic_map_follows(map->p[i], map->p[j]);
+		if (f < 0)
+			return -1;
+		if (!f)
 			continue;
-		isl_int_sub(bmap->div[i][1],
-				bmap->div[i][1], bmap->div[i][1 + pos]);
+
+		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;
 	}
 
-	return bmap;
-}
+	if (s->node[i].index != s->node[i].min_index)
+		return 0;
 
-/* Shift variable at position "pos" up by one.
- * That is, replace the corresponding variable v by v - 1.
- */
-static __isl_give isl_map *map_shift_pos(__isl_take isl_map *map, unsigned pos)
-{
-	int i;
+	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;
 
-	map = isl_map_cow(map);
-	if (!map)
-		return NULL;
-
-	for (i = 0; i < map->n; ++i) {
-		map->p[i] = basic_map_shift_pos(map->p[i], pos);
-		if (!map->p[i])
-			goto error;
-	}
-	ISL_F_CLR(map, ISL_MAP_NORMALIZED);
-	return map;
-error:
-	isl_map_free(map);
-	return NULL;
+	return 0;
 }
 
-/* Check whether the overapproximation of the power of "map" is exactly
- * the power of "map".  Let R be "map" and A_k the overapproximation.
- * The approximation is exact if
+/* 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 the union of the identity map and
+ * an overapproximation of the map
+ * that takes an element from the dom R \times Z to an
+ * element from ran R \times Z, such that the first n coordinates of the
+ * difference between them is a sum of differences between images
+ * and pre-images in one of the R_i and such that the last coordinate
+ * is equal to the number of steps taken.
+ * That is, let
  *
- *	A_1 = R
- *	A_k = A_{k-1} \circ R			k >= 2
+ *	\Delta_i = { y - x | (x, y) in R_i }
  *
- * Since A_k is known to be an overapproximation, we only need to check
+ * then the constructed map is an overapproximation of
  *
- *	A_1 \subset R
- *	A_k \subset A_{k-1} \circ R		k >= 2
+ *	{ (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
+ *				d = (\sum_i k_i \delta_i, \sum_i k_i) and
+ *				x in dom R and x + d in ran R } union
+ *	{ (x) -> (x) }
  *
+ * We first split the map into strongly connected components, perform
+ * the above on each component and the join the results in the correct
+ * order.  The power of each of the components needs to be extended
+ * with the identity map because a path in the global result need
+ * not go through every component.
+ * The final result will then also contain the identity map, but
+ * this part will be removed when the length of the path is forced
+ * to be strictly positive.
  */
-static int check_power_exactness(__isl_take isl_map *map,
-	__isl_take isl_map *app, unsigned param)
+static __isl_give isl_map *construct_power_components(__isl_take isl_dim *dim,
+	__isl_keep isl_map *map, int *exact, int project)
 {
-	int exact;
-	isl_map *app_1;
-	isl_map *app_2;
-
-	app_1 = isl_map_fix_si(isl_map_copy(app), isl_dim_param, param, 1);
+	int i, n;
+	struct isl_map *path = NULL;
+	struct basic_map_sort *s = NULL;
 
-	exact = isl_map_is_subset(app_1, map);
-	isl_map_free(app_1);
+	if (!map)
+		goto error;
+	if (map->n <= 1)
+		return construct_component(dim, map, exact, project);
 
-	if (!exact || exact < 0) {
-		isl_map_free(app);
-		isl_map_free(map);
-		return exact;
+	s = basic_map_sort_alloc(map->ctx, map->n);
+	if (!s)
+		goto error;
+	for (i = map->n - 1; i >= 0; --i) {
+		if (s->node[i].index >= 0)
+			continue;
+		if (power_components_tarjan(s, map, i) < 0)
+			goto error;
 	}
 
-	app_2 = isl_map_lower_bound_si(isl_map_copy(app),
-					isl_dim_param, param, 2);
-	app_1 = map_shift_pos(app, 1 + param);
-	app_1 = isl_map_apply_range(map, app_1);
-
-	exact = isl_map_is_subset(app_2, app_1);
+	i = 0;
+	n = map->n;
+	path = isl_map_identity(isl_dim_domain(isl_dim_copy(dim)));
+	while (n) {
+		struct isl_map *comp;
+		comp = isl_map_alloc_dim(isl_map_get_dim(map), n, 0);
+		while (s->order[i] != -1) {
+			comp = isl_map_add_basic_map(comp,
+				    isl_basic_map_copy(map->p[s->order[i]]));
+			--n;
+			++i;
+		}
+		path = isl_map_apply_range(path,
+			    construct_component(isl_dim_copy(dim), comp,
+						exact, project));
+		isl_map_free(comp);
+		++i;
+	}
 
-	isl_map_free(app_1);
-	isl_map_free(app_2);
+	basic_map_sort_free(s);
+	isl_dim_free(dim);
 
-	return exact;
+	return path;
+error:
+	basic_map_sort_free(s);
+	isl_dim_free(dim);
+	return NULL;
 }
 
-/* Check whether the overapproximation of the power of "map" is exactly
- * the power of "map", possibly after projecting out the power (if "project"
- * is set).
- *
- * If "project" is set and if "steps" can only result in acyclic paths,
- * then we check
+/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D,
+ * construct a map that is an overapproximation of the map
+ * that takes an element from the space D to another
+ * element from the same space, such that the difference between
+ * them is a strictly positive sum of differences between images
+ * and pre-images in one of the R_i.
+ * The number of differences in the sum is equated to parameter "param".
+ * That is, let
  *
- *	A = R \cup (A \circ R)
+ *	\Delta_i = { y - x | (x, y) in R_i }
  *
- * where A is the overapproximation with the power projected out, i.e.,
- * an overapproximation of the transitive closure.
- * More specifically, since A is known to be an overapproximation, we check
+ * then the constructed map is an overapproximation of
  *
- *	A \subset R \cup (A \circ R)
+ *	{ (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
+ *				d = \sum_i k_i \delta_i and k = \sum_i k_i > 0 }
  *
- * Otherwise, we check if the power is exact.
+ * We first 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 int check_exactness(__isl_take isl_map *map, __isl_take isl_map *app,
-	unsigned param, int project)
+static __isl_give isl_map *construct_power(__isl_keep isl_map *map,
+	unsigned param, int *exact, int project)
 {
-	isl_map *test;
-	int exact;
-
-	if (!project)
-		return check_power_exactness(map, app, param);
+	struct isl_map *app = NULL;
+	struct isl_map *diff;
+	struct isl_dim *dim = NULL;
+	unsigned d;
 
-	map = isl_map_project_out(map, isl_dim_param, param, 1);
-	app = isl_map_project_out(app, isl_dim_param, param, 1);
+	if (!map)
+		return NULL;
 
-	test = isl_map_apply_range(isl_map_copy(map), isl_map_copy(app));
-	test = isl_map_union(test, isl_map_copy(map));
+	dim = isl_map_get_dim(map);
 
-	exact = isl_map_is_subset(app, test);
+	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);
 
-	isl_map_free(app);
-	isl_map_free(test);
+	app = construct_power_components(isl_dim_copy(dim), map,
+					exact, project);
 
-	isl_map_free(map);
+	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);
 
-	return exact;
-error:
-	isl_map_free(app);
-	isl_map_free(map);
-	return -1;
+	return app;
 }
 
 /* Compute the positive powers of "map", or an overapproximation.
@@ -655,12 +985,7 @@ static __isl_give isl_map *map_power(__isl_take isl_map *map, unsigned param,
 		isl_map_dim(map, isl_dim_in) == isl_map_dim(map, isl_dim_out),
 		goto error);
 
-	app = construct_power(map, param, exact ? &project : NULL);
-
-	if (exact &&
-	    (*exact = check_exactness(isl_map_copy(map), isl_map_copy(app),
-				      param, project)) < 0)
-		goto error;
+	app = construct_power(map, param, exact, project);
 
 	isl_map_free(map);
 	return app;