2 * Copyright 2010 INRIA Saclay
4 * Use of this software is governed by the GNU LGPLv2.1 license
6 * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
7 * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
12 #include "isl_map_private.h"
15 /* Given a map that represents a path with the length of the path
16 * encoded as the difference between the last output coordindate
17 * and the last input coordinate, set this length to either
18 * exactly "length" (if "exactly" is set) or at least "length"
19 * (if "exactly" is not set).
21 static __isl_give isl_map *set_path_length(__isl_take isl_map *map,
22 int exactly, int length)
25 struct isl_basic_map *bmap;
34 dim = isl_map_get_dim(map);
35 d = isl_dim_size(dim, isl_dim_in);
36 nparam = isl_dim_size(dim, isl_dim_param);
37 bmap = isl_basic_map_alloc_dim(dim, 0, 1, 1);
39 k = isl_basic_map_alloc_equality(bmap);
42 k = isl_basic_map_alloc_inequality(bmap);
47 isl_seq_clr(c, 1 + isl_basic_map_total_dim(bmap));
48 isl_int_set_si(c[0], -length);
49 isl_int_set_si(c[1 + nparam + d - 1], -1);
50 isl_int_set_si(c[1 + nparam + d + d - 1], 1);
52 bmap = isl_basic_map_finalize(bmap);
53 map = isl_map_intersect(map, isl_map_from_basic_map(bmap));
57 isl_basic_map_free(bmap);
62 /* Check whether the overapproximation of the power of "map" is exactly
63 * the power of "map". Let R be "map" and A_k the overapproximation.
64 * The approximation is exact if
67 * A_k = A_{k-1} \circ R k >= 2
69 * Since A_k is known to be an overapproximation, we only need to check
72 * A_k \subset A_{k-1} \circ R k >= 2
74 * In practice, "app" has an extra input and output coordinate
75 * to encode the length of the path. So, we first need to add
76 * this coordinate to "map" and set the length of the path to
79 static int check_power_exactness(__isl_take isl_map *map,
80 __isl_take isl_map *app)
86 map = isl_map_add(map, isl_dim_in, 1);
87 map = isl_map_add(map, isl_dim_out, 1);
88 map = set_path_length(map, 1, 1);
90 app_1 = set_path_length(isl_map_copy(app), 1, 1);
92 exact = isl_map_is_subset(app_1, map);
95 if (!exact || exact < 0) {
101 app_1 = set_path_length(isl_map_copy(app), 0, 1);
102 app_2 = set_path_length(app, 0, 2);
103 app_1 = isl_map_apply_range(map, app_1);
105 exact = isl_map_is_subset(app_2, app_1);
113 /* Check whether the overapproximation of the power of "map" is exactly
114 * the power of "map", possibly after projecting out the power (if "project"
117 * If "project" is set and if "steps" can only result in acyclic paths,
120 * A = R \cup (A \circ R)
122 * where A is the overapproximation with the power projected out, i.e.,
123 * an overapproximation of the transitive closure.
124 * More specifically, since A is known to be an overapproximation, we check
126 * A \subset R \cup (A \circ R)
128 * Otherwise, we check if the power is exact.
130 * Note that "app" has an extra input and output coordinate to encode
131 * the length of the part. If we are only interested in the transitive
132 * closure, then we can simply project out these coordinates first.
134 static int check_exactness(__isl_take isl_map *map, __isl_take isl_map *app,
142 return check_power_exactness(map, app);
144 d = isl_map_dim(map, isl_dim_in);
145 app = set_path_length(app, 0, 1);
146 app = isl_map_project_out(app, isl_dim_in, d, 1);
147 app = isl_map_project_out(app, isl_dim_out, d, 1);
149 test = isl_map_apply_range(isl_map_copy(map), isl_map_copy(app));
150 test = isl_map_union(test, isl_map_copy(map));
152 exact = isl_map_is_subset(app, test);
167 * The transitive closure implementation is based on the paper
168 * "Computing the Transitive Closure of a Union of Affine Integer
169 * Tuple Relations" by Anna Beletska, Denis Barthou, Wlodzimierz Bielecki and
173 /* Given a set of n offsets v_i (the rows of "steps"), construct a relation
174 * of the given dimension specification (Z^{n+1} -> Z^{n+1})
175 * that maps an element x to any element that can be reached
176 * by taking a non-negative number of steps along any of
177 * the extended offsets v'_i = [v_i 1].
180 * { [x] -> [y] : exists k_i >= 0, y = x + \sum_i k_i v'_i }
182 * For any element in this relation, the number of steps taken
183 * is equal to the difference in the final coordinates.
185 static __isl_give isl_map *path_along_steps(__isl_take isl_dim *dim,
186 __isl_keep isl_mat *steps)
189 struct isl_basic_map *path = NULL;
197 d = isl_dim_size(dim, isl_dim_in);
199 nparam = isl_dim_size(dim, isl_dim_param);
201 path = isl_basic_map_alloc_dim(isl_dim_copy(dim), n, d, n);
203 for (i = 0; i < n; ++i) {
204 k = isl_basic_map_alloc_div(path);
207 isl_assert(steps->ctx, i == k, goto error);
208 isl_int_set_si(path->div[k][0], 0);
211 for (i = 0; i < d; ++i) {
212 k = isl_basic_map_alloc_equality(path);
215 isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
216 isl_int_set_si(path->eq[k][1 + nparam + i], 1);
217 isl_int_set_si(path->eq[k][1 + nparam + d + i], -1);
219 for (j = 0; j < n; ++j)
220 isl_int_set_si(path->eq[k][1 + nparam + 2 * d + j], 1);
222 for (j = 0; j < n; ++j)
223 isl_int_set(path->eq[k][1 + nparam + 2 * d + j],
227 for (i = 0; i < n; ++i) {
228 k = isl_basic_map_alloc_inequality(path);
231 isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
232 isl_int_set_si(path->ineq[k][1 + nparam + 2 * d + i], 1);
237 path = isl_basic_map_simplify(path);
238 path = isl_basic_map_finalize(path);
239 return isl_map_from_basic_map(path);
242 isl_basic_map_free(path);
251 /* Return PURE_PARAM if only the coefficients of the parameters are non-zero.
252 * Return PURE_VAR if only the coefficients of the set variables are non-zero.
253 * Return MIXED if only the coefficients of the parameters and the set
254 * variables are non-zero and if moreover the parametric constant
255 * can never attain positive values.
256 * Return IMPURE otherwise.
258 static int purity(__isl_keep isl_basic_set *bset, isl_int *c, int eq)
266 n_div = isl_basic_set_dim(bset, isl_dim_div);
267 d = isl_basic_set_dim(bset, isl_dim_set);
268 nparam = isl_basic_set_dim(bset, isl_dim_param);
270 if (isl_seq_first_non_zero(c + 1 + nparam + d, n_div) != -1)
272 if (isl_seq_first_non_zero(c + 1, nparam) == -1)
274 if (isl_seq_first_non_zero(c + 1 + nparam, d) == -1)
279 bset = isl_basic_set_copy(bset);
280 bset = isl_basic_set_cow(bset);
281 bset = isl_basic_set_extend_constraints(bset, 0, 1);
282 k = isl_basic_set_alloc_inequality(bset);
285 isl_seq_clr(bset->ineq[k], 1 + isl_basic_set_total_dim(bset));
286 isl_seq_cpy(bset->ineq[k], c, 1 + nparam);
287 isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
288 empty = isl_basic_set_is_empty(bset);
289 isl_basic_set_free(bset);
291 return empty < 0 ? -1 : empty ? MIXED : IMPURE;
293 isl_basic_set_free(bset);
297 /* Given a set of offsets "delta", construct a relation of the
298 * given dimension specification (Z^{n+1} -> Z^{n+1}) that
299 * is an overapproximation of the relations that
300 * maps an element x to any element that can be reached
301 * by taking a non-negative number of steps along any of
302 * the elements in "delta".
303 * That is, construct an approximation of
305 * { [x] -> [y] : exists f \in \delta, k \in Z :
306 * y = x + k [f, 1] and k >= 0 }
308 * For any element in this relation, the number of steps taken
309 * is equal to the difference in the final coordinates.
311 * In particular, let delta be defined as
313 * \delta = [p] -> { [x] : A x + a >= and B p + b >= 0 and
314 * C x + C'p + c >= 0 and
315 * D x + D'p + d >= 0 }
317 * where the constraints C x + C'p + c >= 0 are such that the parametric
318 * constant term of each constraint j, "C_j x + C'_j p + c_j",
319 * can never attain positive values, then the relation is constructed as
321 * { [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and
322 * A f + k a >= 0 and B p + b >= 0 and
323 * C f + C'p + c >= 0 and k >= 1 }
324 * union { [x] -> [x] }
326 * Existentially quantified variables in \delta are currently ignored.
327 * This is safe, but leads to an additional overapproximation.
329 static __isl_give isl_map *path_along_delta(__isl_take isl_dim *dim,
330 __isl_take isl_basic_set *delta)
332 isl_basic_map *path = NULL;
341 n_div = isl_basic_set_dim(delta, isl_dim_div);
342 d = isl_basic_set_dim(delta, isl_dim_set);
343 nparam = isl_basic_set_dim(delta, isl_dim_param);
344 path = isl_basic_map_alloc_dim(isl_dim_copy(dim), n_div + d + 1,
345 d + 1 + delta->n_eq, delta->n_ineq + 1);
346 off = 1 + nparam + 2 * (d + 1) + n_div;
348 for (i = 0; i < n_div + d + 1; ++i) {
349 k = isl_basic_map_alloc_div(path);
352 isl_int_set_si(path->div[k][0], 0);
355 for (i = 0; i < d + 1; ++i) {
356 k = isl_basic_map_alloc_equality(path);
359 isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
360 isl_int_set_si(path->eq[k][1 + nparam + i], 1);
361 isl_int_set_si(path->eq[k][1 + nparam + d + 1 + i], -1);
362 isl_int_set_si(path->eq[k][off + i], 1);
365 for (i = 0; i < delta->n_eq; ++i) {
366 int p = purity(delta, delta->eq[i], 1);
371 k = isl_basic_map_alloc_equality(path);
374 isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
376 isl_seq_cpy(path->eq[k] + off,
377 delta->eq[i] + 1 + nparam, d);
378 isl_int_set(path->eq[k][off + d], delta->eq[i][0]);
380 isl_seq_cpy(path->eq[k], delta->eq[i], 1 + nparam);
383 for (i = 0; i < delta->n_ineq; ++i) {
384 int p = purity(delta, delta->ineq[i], 0);
389 k = isl_basic_map_alloc_inequality(path);
392 isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
394 isl_seq_cpy(path->ineq[k] + off,
395 delta->ineq[i] + 1 + nparam, d);
396 isl_int_set(path->ineq[k][off + d], delta->ineq[i][0]);
397 } else if (p == PURE_PARAM) {
398 isl_seq_cpy(path->ineq[k], delta->ineq[i], 1 + nparam);
400 isl_seq_cpy(path->ineq[k] + off,
401 delta->ineq[i] + 1 + nparam, d);
402 isl_seq_cpy(path->ineq[k], delta->ineq[i], 1 + nparam);
406 k = isl_basic_map_alloc_inequality(path);
409 isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
410 isl_int_set_si(path->ineq[k][0], -1);
411 isl_int_set_si(path->ineq[k][off + d], 1);
413 isl_basic_set_free(delta);
414 path = isl_basic_map_finalize(path);
415 return isl_basic_map_union(path,
416 isl_basic_map_identity(isl_dim_domain(dim)));
419 isl_basic_set_free(delta);
420 isl_basic_map_free(path);
424 /* Given a dimenion specification Z^{n+1} -> Z^{n+1} and a parameter "param",
425 * construct a map that equates the parameter to the difference
426 * in the final coordinates and imposes that this difference is positive.
429 * { [x,x_s] -> [y,y_s] : k = y_s - x_s > 0 }
431 static __isl_give isl_map *equate_parameter_to_length(__isl_take isl_dim *dim,
434 struct isl_basic_map *bmap;
439 d = isl_dim_size(dim, isl_dim_in);
440 nparam = isl_dim_size(dim, isl_dim_param);
441 bmap = isl_basic_map_alloc_dim(dim, 0, 1, 1);
442 k = isl_basic_map_alloc_equality(bmap);
445 isl_seq_clr(bmap->eq[k], 1 + isl_basic_map_total_dim(bmap));
446 isl_int_set_si(bmap->eq[k][1 + param], -1);
447 isl_int_set_si(bmap->eq[k][1 + nparam + d - 1], -1);
448 isl_int_set_si(bmap->eq[k][1 + nparam + d + d - 1], 1);
450 k = isl_basic_map_alloc_inequality(bmap);
453 isl_seq_clr(bmap->ineq[k], 1 + isl_basic_map_total_dim(bmap));
454 isl_int_set_si(bmap->ineq[k][1 + param], 1);
455 isl_int_set_si(bmap->ineq[k][0], -1);
457 bmap = isl_basic_map_finalize(bmap);
458 return isl_map_from_basic_map(bmap);
460 isl_basic_map_free(bmap);
464 /* Check whether "path" is acyclic, where the last coordinates of domain
465 * and range of path encode the number of steps taken.
466 * That is, check whether
468 * { d | d = y - x and (x,y) in path }
470 * does not contain any element with positive last coordinate (positive length)
471 * and zero remaining coordinates (cycle).
473 static int is_acyclic(__isl_take isl_map *path)
478 struct isl_set *delta;
480 delta = isl_map_deltas(path);
481 dim = isl_set_dim(delta, isl_dim_set);
482 for (i = 0; i < dim; ++i) {
484 delta = isl_set_lower_bound_si(delta, isl_dim_set, i, 1);
486 delta = isl_set_fix_si(delta, isl_dim_set, i, 0);
489 acyclic = isl_set_is_empty(delta);
495 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
496 * and a dimension specification (Z^{n+1} -> Z^{n+1}),
497 * construct a map that is an overapproximation of the map
498 * that takes an element from the space D \times Z to another
499 * element from the same space, such that the first n coordinates of the
500 * difference between them is a sum of differences between images
501 * and pre-images in one of the R_i and such that the last coordinate
502 * is equal to the number of steps taken.
505 * \Delta_i = { y - x | (x, y) in R_i }
507 * then the constructed map is an overapproximation of
509 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
510 * d = (\sum_i k_i \delta_i, \sum_i k_i) }
512 * The elements of the singleton \Delta_i's are collected as the
513 * rows of the steps matrix. For all these \Delta_i's together,
514 * a single path is constructed.
515 * For each of the other \Delta_i's, we compute an overapproximation
516 * of the paths along elements of \Delta_i.
517 * Since each of these paths performs an addition, composition is
518 * symmetric and we can simply compose all resulting paths in any order.
520 static __isl_give isl_map *construct_extended_path(__isl_take isl_dim *dim,
521 __isl_keep isl_map *map, int *project)
523 struct isl_mat *steps = NULL;
524 struct isl_map *path = NULL;
528 d = isl_map_dim(map, isl_dim_in);
530 path = isl_map_identity(isl_dim_domain(isl_dim_copy(dim)));
532 steps = isl_mat_alloc(map->ctx, map->n, d);
537 for (i = 0; i < map->n; ++i) {
538 struct isl_basic_set *delta;
540 delta = isl_basic_map_deltas(isl_basic_map_copy(map->p[i]));
542 for (j = 0; j < d; ++j) {
545 fixed = isl_basic_set_fast_dim_is_fixed(delta, j,
548 isl_basic_set_free(delta);
557 path = isl_map_apply_range(path,
558 path_along_delta(isl_dim_copy(dim), delta));
560 isl_basic_set_free(delta);
567 path = isl_map_apply_range(path,
568 path_along_steps(isl_dim_copy(dim), steps));
571 if (project && *project) {
572 *project = is_acyclic(isl_map_copy(path));
587 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
588 * and a dimension specification (Z^{n+1} -> Z^{n+1}),
589 * construct a map that is the union of the identity map and
590 * an overapproximation of the map
591 * that takes an element from the dom R \times Z to an
592 * element from ran R \times Z, such that the first n coordinates of the
593 * difference between them is a sum of differences between images
594 * and pre-images in one of the R_i and such that the last coordinate
595 * is equal to the number of steps taken.
598 * \Delta_i = { y - x | (x, y) in R_i }
600 * then the constructed map is an overapproximation of
602 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
603 * d = (\sum_i k_i \delta_i, \sum_i k_i) and
604 * x in dom R and x + d in ran R } union
607 static __isl_give isl_map *construct_component(__isl_take isl_dim *dim,
608 __isl_keep isl_map *map, int *exact, int project)
610 struct isl_set *domain = NULL;
611 struct isl_set *range = NULL;
612 struct isl_map *app = NULL;
613 struct isl_map *path = NULL;
615 domain = isl_map_domain(isl_map_copy(map));
616 domain = isl_set_coalesce(domain);
617 range = isl_map_range(isl_map_copy(map));
618 range = isl_set_coalesce(range);
619 app = isl_map_from_domain_and_range(domain, range);
620 app = isl_map_add(app, isl_dim_in, 1);
621 app = isl_map_add(app, isl_dim_out, 1);
623 path = construct_extended_path(isl_dim_copy(dim), map,
624 exact && *exact ? &project : NULL);
625 app = isl_map_intersect(app, path);
627 if (exact && *exact &&
628 (*exact = check_exactness(isl_map_copy(map), isl_map_copy(app),
632 return isl_map_union(app, isl_map_identity(isl_dim_domain(dim)));
639 /* Structure for representing the nodes in the graph being traversed
640 * using Tarjan's algorithm.
641 * index represents the order in which nodes are visited.
642 * min_index is the index of the root of a (sub)component.
643 * on_stack indicates whether the node is currently on the stack.
645 struct basic_map_sort_node {
650 /* Structure for representing the graph being traversed
651 * using Tarjan's algorithm.
652 * len is the number of nodes
653 * node is an array of nodes
654 * stack contains the nodes on the path from the root to the current node
655 * sp is the stack pointer
656 * index is the index of the last node visited
657 * order contains the elements of the components separated by -1
658 * op represents the current position in order
660 struct basic_map_sort {
662 struct basic_map_sort_node *node;
670 static void basic_map_sort_free(struct basic_map_sort *s)
680 static struct basic_map_sort *basic_map_sort_alloc(struct isl_ctx *ctx, int len)
682 struct basic_map_sort *s;
685 s = isl_calloc_type(ctx, struct basic_map_sort);
689 s->node = isl_alloc_array(ctx, struct basic_map_sort_node, len);
692 for (i = 0; i < len; ++i)
693 s->node[i].index = -1;
694 s->stack = isl_alloc_array(ctx, int, len);
697 s->order = isl_alloc_array(ctx, int, 2 * len);
707 basic_map_sort_free(s);
711 /* Check whether in the computation of the transitive closure
712 * "bmap1" (R_1) should follow (or be part of the same component as)
715 * That is check whether
719 * is non-empty and that moreover, it is non-empty on the set
720 * of elements that do not get mapped to the same set of elements
721 * by both "R_1 \circ R_2" and "R_2 \circ R_1".
722 * For elements that do get mapped to the same elements by these
723 * two compositions, R_1 and R_2 are commutative, so if these
724 * elements are the only ones for which R_1 \circ R_2 is non-empty,
725 * then you may just as well apply R_1 first.
727 static int basic_map_follows(__isl_keep isl_basic_map *bmap1,
728 __isl_keep isl_basic_map *bmap2)
730 struct isl_map *map12 = NULL;
731 struct isl_map *map21 = NULL;
732 struct isl_map *d = NULL;
733 struct isl_set *dom = NULL;
736 map21 = isl_map_from_basic_map(
737 isl_basic_map_apply_range(
738 isl_basic_map_copy(bmap2),
739 isl_basic_map_copy(bmap1)));
740 empty = isl_map_is_empty(map21);
748 map12 = isl_map_from_basic_map(
749 isl_basic_map_apply_range(
750 isl_basic_map_copy(bmap1),
751 isl_basic_map_copy(bmap2)));
752 d = isl_map_subtract(isl_map_copy(map12), isl_map_copy(map21));
754 isl_map_subtract(isl_map_copy(map21), isl_map_copy(map12)));
755 dom = isl_map_domain(d);
757 map21 = isl_map_intersect_domain(map21, dom);
758 empty = isl_map_is_empty(map21);
763 return empty < 0 ? -1 : !empty;
769 /* Perform Tarjan's algorithm for computing the strongly connected components
770 * in the graph with the disjuncts of "map" as vertices and with an
771 * edge between any pair of disjuncts such that the first has
772 * to be applied after the second.
774 static int power_components_tarjan(struct basic_map_sort *s,
775 __isl_keep isl_map *map, int i)
779 s->node[i].index = s->index;
780 s->node[i].min_index = s->index;
781 s->node[i].on_stack = 1;
783 s->stack[s->sp++] = i;
785 for (j = s->len - 1; j >= 0; --j) {
790 if (s->node[j].index >= 0 &&
791 (!s->node[j].on_stack ||
792 s->node[j].index > s->node[i].min_index))
795 f = basic_map_follows(map->p[i], map->p[j]);
801 if (s->node[j].index < 0) {
802 power_components_tarjan(s, map, j);
803 if (s->node[j].min_index < s->node[i].min_index)
804 s->node[i].min_index = s->node[j].min_index;
805 } else if (s->node[j].index < s->node[i].min_index)
806 s->node[i].min_index = s->node[j].index;
809 if (s->node[i].index != s->node[i].min_index)
813 j = s->stack[--s->sp];
814 s->node[j].on_stack = 0;
815 s->order[s->op++] = j;
817 s->order[s->op++] = -1;
822 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
823 * and a dimension specification (Z^{n+1} -> Z^{n+1}),
824 * construct a map that is the union of the identity map and
825 * an overapproximation of the map
826 * that takes an element from the dom R \times Z to an
827 * element from ran R \times Z, such that the first n coordinates of the
828 * difference between them is a sum of differences between images
829 * and pre-images in one of the R_i and such that the last coordinate
830 * is equal to the number of steps taken.
833 * \Delta_i = { y - x | (x, y) in R_i }
835 * then the constructed map is an overapproximation of
837 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
838 * d = (\sum_i k_i \delta_i, \sum_i k_i) and
839 * x in dom R and x + d in ran R } union
842 * We first split the map into strongly connected components, perform
843 * the above on each component and the join the results in the correct
844 * order. The power of each of the components needs to be extended
845 * with the identity map because a path in the global result need
846 * not go through every component.
847 * The final result will then also contain the identity map, but
848 * this part will be removed when the length of the path is forced
849 * to be strictly positive.
851 static __isl_give isl_map *construct_power_components(__isl_take isl_dim *dim,
852 __isl_keep isl_map *map, int *exact, int project)
855 struct isl_map *path = NULL;
856 struct basic_map_sort *s = NULL;
861 return construct_component(dim, map, exact, project);
863 s = basic_map_sort_alloc(map->ctx, map->n);
866 for (i = map->n - 1; i >= 0; --i) {
867 if (s->node[i].index >= 0)
869 if (power_components_tarjan(s, map, i) < 0)
875 path = isl_map_identity(isl_dim_domain(isl_dim_copy(dim)));
877 struct isl_map *comp;
878 comp = isl_map_alloc_dim(isl_map_get_dim(map), n, 0);
879 while (s->order[i] != -1) {
880 comp = isl_map_add_basic_map(comp,
881 isl_basic_map_copy(map->p[s->order[i]]));
885 path = isl_map_apply_range(path,
886 construct_component(isl_dim_copy(dim), comp,
892 basic_map_sort_free(s);
897 basic_map_sort_free(s);
902 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D,
903 * construct a map that is an overapproximation of the map
904 * that takes an element from the space D to another
905 * element from the same space, such that the difference between
906 * them is a strictly positive sum of differences between images
907 * and pre-images in one of the R_i.
908 * The number of differences in the sum is equated to parameter "param".
911 * \Delta_i = { y - x | (x, y) in R_i }
913 * then the constructed map is an overapproximation of
915 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
916 * d = \sum_i k_i \delta_i and k = \sum_i k_i > 0 }
918 * We first construct an extended mapping with an extra coordinate
919 * that indicates the number of steps taken. In particular,
920 * the difference in the last coordinate is equal to the number
921 * of steps taken to move from a domain element to the corresponding
923 * In the final step, this difference is equated to the parameter "param"
924 * and made positive. The extra coordinates are subsequently projected out.
926 static __isl_give isl_map *construct_power(__isl_keep isl_map *map,
927 unsigned param, int *exact, int project)
929 struct isl_map *app = NULL;
930 struct isl_map *diff;
931 struct isl_dim *dim = NULL;
937 dim = isl_map_get_dim(map);
939 d = isl_dim_size(dim, isl_dim_in);
940 dim = isl_dim_add(dim, isl_dim_in, 1);
941 dim = isl_dim_add(dim, isl_dim_out, 1);
943 app = construct_power_components(isl_dim_copy(dim), map,
946 diff = equate_parameter_to_length(dim, param);
947 app = isl_map_intersect(app, diff);
948 app = isl_map_project_out(app, isl_dim_in, d, 1);
949 app = isl_map_project_out(app, isl_dim_out, d, 1);
954 /* Compute the positive powers of "map", or an overapproximation.
955 * The power is given by parameter "param". If the result is exact,
956 * then *exact is set to 1.
957 * If project is set, then we are actually interested in the transitive
958 * closure, so we can use a more relaxed exactness check.
960 static __isl_give isl_map *map_power(__isl_take isl_map *map, unsigned param,
961 int *exact, int project)
963 struct isl_map *app = NULL;
968 map = isl_map_remove_empty_parts(map);
972 if (isl_map_fast_is_empty(map))
975 isl_assert(map->ctx, param < isl_map_dim(map, isl_dim_param), goto error);
977 isl_map_dim(map, isl_dim_in) == isl_map_dim(map, isl_dim_out),
980 app = construct_power(map, param, exact, project);
990 /* Compute the positive powers of "map", or an overapproximation.
991 * The power is given by parameter "param". If the result is exact,
992 * then *exact is set to 1.
994 __isl_give isl_map *isl_map_power(__isl_take isl_map *map, unsigned param,
997 return map_power(map, param, exact, 0);
1000 /* Compute the transitive closure of "map", or an overapproximation.
1001 * If the result is exact, then *exact is set to 1.
1002 * Simply compute the powers of map and then project out the parameter
1003 * describing the power.
1005 __isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
1013 param = isl_map_dim(map, isl_dim_param);
1014 map = isl_map_add(map, isl_dim_param, 1);
1015 map = map_power(map, param, exact, 1);
1016 map = isl_map_project_out(map, isl_dim_param, param, 1);