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 path with the as yet unconstrained length at position "pos",
298 * check if setting the length to zero results in only the identity
301 int empty_path_is_identity(__isl_keep isl_basic_map *path, unsigned pos)
303 isl_basic_map *test = NULL;
304 isl_basic_map *id = NULL;
308 test = isl_basic_map_copy(path);
309 test = isl_basic_map_extend_constraints(test, 1, 0);
310 k = isl_basic_map_alloc_equality(test);
313 isl_seq_clr(test->eq[k], 1 + isl_basic_map_total_dim(test));
314 isl_int_set_si(test->eq[k][pos], 1);
315 id = isl_basic_map_identity(isl_dim_domain(isl_basic_map_get_dim(path)));
316 is_id = isl_basic_map_is_subset(test, id);
317 isl_basic_map_free(test);
318 isl_basic_map_free(id);
321 isl_basic_map_free(test);
325 __isl_give isl_basic_map *add_delta_constraints(__isl_take isl_basic_map *path,
326 __isl_keep isl_basic_set *delta, unsigned off, unsigned nparam,
330 int n = eq ? delta->n_eq : delta->n_ineq;
331 isl_int **delta_c = eq ? delta->eq : delta->ineq;
332 isl_int **path_c = eq ? path->eq : path->ineq;
334 for (i = 0; i < n; ++i) {
335 int p = purity(delta, delta_c[i], eq);
341 k = isl_basic_map_alloc_equality(path);
343 k = isl_basic_map_alloc_inequality(path);
346 isl_seq_clr(path_c[k], 1 + isl_basic_map_total_dim(path));
348 isl_seq_cpy(path_c[k] + off,
349 delta_c[i] + 1 + nparam, d);
350 isl_int_set(path_c[k][off + d], delta_c[i][0]);
351 } else if (p == PURE_PARAM) {
352 isl_seq_cpy(path_c[k], delta_c[i], 1 + nparam);
354 isl_seq_cpy(path_c[k] + off,
355 delta_c[i] + 1 + nparam, d);
356 isl_seq_cpy(path_c[k], delta_c[i], 1 + nparam);
362 isl_basic_map_free(path);
366 /* Given a set of offsets "delta", construct a relation of the
367 * given dimension specification (Z^{n+1} -> Z^{n+1}) that
368 * is an overapproximation of the relations that
369 * maps an element x to any element that can be reached
370 * by taking a non-negative number of steps along any of
371 * the elements in "delta".
372 * That is, construct an approximation of
374 * { [x] -> [y] : exists f \in \delta, k \in Z :
375 * y = x + k [f, 1] and k >= 0 }
377 * For any element in this relation, the number of steps taken
378 * is equal to the difference in the final coordinates.
380 * In particular, let delta be defined as
382 * \delta = [p] -> { [x] : A x + a >= and B p + b >= 0 and
383 * C x + C'p + c >= 0 and
384 * D x + D'p + d >= 0 }
386 * where the constraints C x + C'p + c >= 0 are such that the parametric
387 * constant term of each constraint j, "C_j x + C'_j p + c_j",
388 * can never attain positive values, then the relation is constructed as
390 * { [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and
391 * A f + k a >= 0 and B p + b >= 0 and
392 * C f + C'p + c >= 0 and k >= 1 }
393 * union { [x] -> [x] }
395 * If the zero-length paths happen to correspond exactly to the identity
396 * mapping, then we return
398 * { [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and
399 * A f + k a >= 0 and B p + b >= 0 and
400 * C f + C'p + c >= 0 and k >= 0 }
404 * Existentially quantified variables in \delta are currently ignored.
405 * This is safe, but leads to an additional overapproximation.
407 static __isl_give isl_map *path_along_delta(__isl_take isl_dim *dim,
408 __isl_take isl_basic_set *delta)
410 isl_basic_map *path = NULL;
420 n_div = isl_basic_set_dim(delta, isl_dim_div);
421 d = isl_basic_set_dim(delta, isl_dim_set);
422 nparam = isl_basic_set_dim(delta, isl_dim_param);
423 path = isl_basic_map_alloc_dim(isl_dim_copy(dim), n_div + d + 1,
424 d + 1 + delta->n_eq, delta->n_ineq + 1);
425 off = 1 + nparam + 2 * (d + 1) + n_div;
427 for (i = 0; i < n_div + d + 1; ++i) {
428 k = isl_basic_map_alloc_div(path);
431 isl_int_set_si(path->div[k][0], 0);
434 for (i = 0; i < d + 1; ++i) {
435 k = isl_basic_map_alloc_equality(path);
438 isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
439 isl_int_set_si(path->eq[k][1 + nparam + i], 1);
440 isl_int_set_si(path->eq[k][1 + nparam + d + 1 + i], -1);
441 isl_int_set_si(path->eq[k][off + i], 1);
444 path = add_delta_constraints(path, delta, off, nparam, d, 1);
445 path = add_delta_constraints(path, delta, off, nparam, d, 0);
447 is_id = empty_path_is_identity(path, off + d);
451 k = isl_basic_map_alloc_inequality(path);
454 isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
456 isl_int_set_si(path->ineq[k][0], -1);
457 isl_int_set_si(path->ineq[k][off + d], 1);
459 isl_basic_set_free(delta);
460 path = isl_basic_map_finalize(path);
463 return isl_map_from_basic_map(path);
465 return isl_basic_map_union(path,
466 isl_basic_map_identity(isl_dim_domain(dim)));
469 isl_basic_set_free(delta);
470 isl_basic_map_free(path);
474 /* Given a dimenion specification Z^{n+1} -> Z^{n+1} and a parameter "param",
475 * construct a map that equates the parameter to the difference
476 * in the final coordinates and imposes that this difference is positive.
479 * { [x,x_s] -> [y,y_s] : k = y_s - x_s > 0 }
481 static __isl_give isl_map *equate_parameter_to_length(__isl_take isl_dim *dim,
484 struct isl_basic_map *bmap;
489 d = isl_dim_size(dim, isl_dim_in);
490 nparam = isl_dim_size(dim, isl_dim_param);
491 bmap = isl_basic_map_alloc_dim(dim, 0, 1, 1);
492 k = isl_basic_map_alloc_equality(bmap);
495 isl_seq_clr(bmap->eq[k], 1 + isl_basic_map_total_dim(bmap));
496 isl_int_set_si(bmap->eq[k][1 + param], -1);
497 isl_int_set_si(bmap->eq[k][1 + nparam + d - 1], -1);
498 isl_int_set_si(bmap->eq[k][1 + nparam + d + d - 1], 1);
500 k = isl_basic_map_alloc_inequality(bmap);
503 isl_seq_clr(bmap->ineq[k], 1 + isl_basic_map_total_dim(bmap));
504 isl_int_set_si(bmap->ineq[k][1 + param], 1);
505 isl_int_set_si(bmap->ineq[k][0], -1);
507 bmap = isl_basic_map_finalize(bmap);
508 return isl_map_from_basic_map(bmap);
510 isl_basic_map_free(bmap);
514 /* Check whether "path" is acyclic, where the last coordinates of domain
515 * and range of path encode the number of steps taken.
516 * That is, check whether
518 * { d | d = y - x and (x,y) in path }
520 * does not contain any element with positive last coordinate (positive length)
521 * and zero remaining coordinates (cycle).
523 static int is_acyclic(__isl_take isl_map *path)
528 struct isl_set *delta;
530 delta = isl_map_deltas(path);
531 dim = isl_set_dim(delta, isl_dim_set);
532 for (i = 0; i < dim; ++i) {
534 delta = isl_set_lower_bound_si(delta, isl_dim_set, i, 1);
536 delta = isl_set_fix_si(delta, isl_dim_set, i, 0);
539 acyclic = isl_set_is_empty(delta);
545 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
546 * and a dimension specification (Z^{n+1} -> Z^{n+1}),
547 * construct a map that is an overapproximation of the map
548 * that takes an element from the space D \times Z to another
549 * element from the same space, such that the first n coordinates of the
550 * difference between them is a sum of differences between images
551 * and pre-images in one of the R_i and such that the last coordinate
552 * is equal to the number of steps taken.
555 * \Delta_i = { y - x | (x, y) in R_i }
557 * then the constructed map is an overapproximation of
559 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
560 * d = (\sum_i k_i \delta_i, \sum_i k_i) }
562 * The elements of the singleton \Delta_i's are collected as the
563 * rows of the steps matrix. For all these \Delta_i's together,
564 * a single path is constructed.
565 * For each of the other \Delta_i's, we compute an overapproximation
566 * of the paths along elements of \Delta_i.
567 * Since each of these paths performs an addition, composition is
568 * symmetric and we can simply compose all resulting paths in any order.
570 static __isl_give isl_map *construct_extended_path(__isl_take isl_dim *dim,
571 __isl_keep isl_map *map, int *project)
573 struct isl_mat *steps = NULL;
574 struct isl_map *path = NULL;
578 d = isl_map_dim(map, isl_dim_in);
580 path = isl_map_identity(isl_dim_domain(isl_dim_copy(dim)));
582 steps = isl_mat_alloc(map->ctx, map->n, d);
587 for (i = 0; i < map->n; ++i) {
588 struct isl_basic_set *delta;
590 delta = isl_basic_map_deltas(isl_basic_map_copy(map->p[i]));
592 for (j = 0; j < d; ++j) {
595 fixed = isl_basic_set_fast_dim_is_fixed(delta, j,
598 isl_basic_set_free(delta);
607 path = isl_map_apply_range(path,
608 path_along_delta(isl_dim_copy(dim), delta));
609 path = isl_map_coalesce(path);
611 isl_basic_set_free(delta);
618 path = isl_map_apply_range(path,
619 path_along_steps(isl_dim_copy(dim), steps));
622 if (project && *project) {
623 *project = is_acyclic(isl_map_copy(path));
638 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
639 * and a dimension specification (Z^{n+1} -> Z^{n+1}),
640 * construct a map that is the union of the identity map and
641 * an overapproximation of the map
642 * that takes an element from the dom R \times Z to an
643 * element from ran R \times Z, such that the first n coordinates of the
644 * difference between them is a sum of differences between images
645 * and pre-images in one of the R_i and such that the last coordinate
646 * is equal to the number of steps taken.
649 * \Delta_i = { y - x | (x, y) in R_i }
651 * then the constructed map is an overapproximation of
653 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
654 * d = (\sum_i k_i \delta_i, \sum_i k_i) and
655 * x in dom R and x + d in ran R } union
658 static __isl_give isl_map *construct_component(__isl_take isl_dim *dim,
659 __isl_keep isl_map *map, int *exact, int project)
661 struct isl_set *domain = NULL;
662 struct isl_set *range = NULL;
663 struct isl_set *overlap;
664 struct isl_map *app = NULL;
665 struct isl_map *path = NULL;
667 domain = isl_map_domain(isl_map_copy(map));
668 domain = isl_set_coalesce(domain);
669 range = isl_map_range(isl_map_copy(map));
670 range = isl_set_coalesce(range);
671 overlap = isl_set_intersect(isl_set_copy(domain), isl_set_copy(range));
672 if (isl_set_is_empty(overlap) == 1) {
673 isl_set_free(domain);
675 isl_set_free(overlap);
678 map = isl_map_copy(map);
679 map = isl_map_add(map, isl_dim_in, 1);
680 map = isl_map_add(map, isl_dim_out, 1);
681 map = set_path_length(map, 1, 1);
684 isl_set_free(overlap);
685 app = isl_map_from_domain_and_range(domain, range);
686 app = isl_map_add(app, isl_dim_in, 1);
687 app = isl_map_add(app, isl_dim_out, 1);
689 path = construct_extended_path(isl_dim_copy(dim), map,
690 exact && *exact ? &project : NULL);
691 app = isl_map_intersect(app, path);
693 if (exact && *exact &&
694 (*exact = check_exactness(isl_map_copy(map), isl_map_copy(app),
698 return isl_map_union(app, isl_map_identity(isl_dim_domain(dim)));
705 /* Structure for representing the nodes in the graph being traversed
706 * using Tarjan's algorithm.
707 * index represents the order in which nodes are visited.
708 * min_index is the index of the root of a (sub)component.
709 * on_stack indicates whether the node is currently on the stack.
711 struct basic_map_sort_node {
716 /* Structure for representing the graph being traversed
717 * using Tarjan's algorithm.
718 * len is the number of nodes
719 * node is an array of nodes
720 * stack contains the nodes on the path from the root to the current node
721 * sp is the stack pointer
722 * index is the index of the last node visited
723 * order contains the elements of the components separated by -1
724 * op represents the current position in order
726 struct basic_map_sort {
728 struct basic_map_sort_node *node;
736 static void basic_map_sort_free(struct basic_map_sort *s)
746 static struct basic_map_sort *basic_map_sort_alloc(struct isl_ctx *ctx, int len)
748 struct basic_map_sort *s;
751 s = isl_calloc_type(ctx, struct basic_map_sort);
755 s->node = isl_alloc_array(ctx, struct basic_map_sort_node, len);
758 for (i = 0; i < len; ++i)
759 s->node[i].index = -1;
760 s->stack = isl_alloc_array(ctx, int, len);
763 s->order = isl_alloc_array(ctx, int, 2 * len);
773 basic_map_sort_free(s);
777 /* Check whether in the computation of the transitive closure
778 * "bmap1" (R_1) should follow (or be part of the same component as)
781 * That is check whether
789 * If so, then there is no reason for R_1 to immediately follow R_2
792 static int basic_map_follows(__isl_keep isl_basic_map *bmap1,
793 __isl_keep isl_basic_map *bmap2)
795 struct isl_map *map12 = NULL;
796 struct isl_map *map21 = NULL;
799 map21 = isl_map_from_basic_map(
800 isl_basic_map_apply_range(
801 isl_basic_map_copy(bmap2),
802 isl_basic_map_copy(bmap1)));
803 subset = isl_map_is_empty(map21);
811 map12 = isl_map_from_basic_map(
812 isl_basic_map_apply_range(
813 isl_basic_map_copy(bmap1),
814 isl_basic_map_copy(bmap2)));
816 subset = isl_map_is_subset(map21, map12);
821 return subset < 0 ? -1 : !subset;
827 /* Perform Tarjan's algorithm for computing the strongly connected components
828 * in the graph with the disjuncts of "map" as vertices and with an
829 * edge between any pair of disjuncts such that the first has
830 * to be applied after the second.
832 static int power_components_tarjan(struct basic_map_sort *s,
833 __isl_keep isl_map *map, int i)
837 s->node[i].index = s->index;
838 s->node[i].min_index = s->index;
839 s->node[i].on_stack = 1;
841 s->stack[s->sp++] = i;
843 for (j = s->len - 1; j >= 0; --j) {
848 if (s->node[j].index >= 0 &&
849 (!s->node[j].on_stack ||
850 s->node[j].index > s->node[i].min_index))
853 f = basic_map_follows(map->p[i], map->p[j]);
859 if (s->node[j].index < 0) {
860 power_components_tarjan(s, map, j);
861 if (s->node[j].min_index < s->node[i].min_index)
862 s->node[i].min_index = s->node[j].min_index;
863 } else if (s->node[j].index < s->node[i].min_index)
864 s->node[i].min_index = s->node[j].index;
867 if (s->node[i].index != s->node[i].min_index)
871 j = s->stack[--s->sp];
872 s->node[j].on_stack = 0;
873 s->order[s->op++] = j;
875 s->order[s->op++] = -1;
880 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
881 * and a dimension specification (Z^{n+1} -> Z^{n+1}),
882 * construct a map that is the union of the identity map and
883 * an overapproximation of the map
884 * that takes an element from the dom R \times Z to an
885 * element from ran R \times Z, such that the first n coordinates of the
886 * difference between them is a sum of differences between images
887 * and pre-images in one of the R_i and such that the last coordinate
888 * is equal to the number of steps taken.
891 * \Delta_i = { y - x | (x, y) in R_i }
893 * then the constructed map is an overapproximation of
895 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
896 * d = (\sum_i k_i \delta_i, \sum_i k_i) and
897 * x in dom R and x + d in ran R } union
900 * We first split the map into strongly connected components, perform
901 * the above on each component and the join the results in the correct
902 * order. The power of each of the components needs to be extended
903 * with the identity map because a path in the global result need
904 * not go through every component.
905 * The final result will then also contain the identity map, but
906 * this part will be removed when the length of the path is forced
907 * to be strictly positive.
909 static __isl_give isl_map *construct_power_components(__isl_take isl_dim *dim,
910 __isl_keep isl_map *map, int *exact, int project)
913 struct isl_map *path = NULL;
914 struct basic_map_sort *s = NULL;
919 return construct_component(dim, map, exact, project);
921 s = basic_map_sort_alloc(map->ctx, map->n);
924 for (i = map->n - 1; i >= 0; --i) {
925 if (s->node[i].index >= 0)
927 if (power_components_tarjan(s, map, i) < 0)
933 path = isl_map_identity(isl_dim_domain(isl_dim_copy(dim)));
935 struct isl_map *comp;
936 comp = isl_map_alloc_dim(isl_map_get_dim(map), n, 0);
937 while (s->order[i] != -1) {
938 comp = isl_map_add_basic_map(comp,
939 isl_basic_map_copy(map->p[s->order[i]]));
943 path = isl_map_apply_range(path,
944 construct_component(isl_dim_copy(dim), comp,
950 basic_map_sort_free(s);
955 basic_map_sort_free(s);
960 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D,
961 * construct a map that is an overapproximation of the map
962 * that takes an element from the space D to another
963 * element from the same space, such that the difference between
964 * them is a strictly positive sum of differences between images
965 * and pre-images in one of the R_i.
966 * The number of differences in the sum is equated to parameter "param".
969 * \Delta_i = { y - x | (x, y) in R_i }
971 * then the constructed map is an overapproximation of
973 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
974 * d = \sum_i k_i \delta_i and k = \sum_i k_i > 0 }
976 * We first construct an extended mapping with an extra coordinate
977 * that indicates the number of steps taken. In particular,
978 * the difference in the last coordinate is equal to the number
979 * of steps taken to move from a domain element to the corresponding
981 * In the final step, this difference is equated to the parameter "param"
982 * and made positive. The extra coordinates are subsequently projected out.
984 static __isl_give isl_map *construct_power(__isl_keep isl_map *map,
985 unsigned param, int *exact, int project)
987 struct isl_map *app = NULL;
988 struct isl_map *diff;
989 struct isl_dim *dim = NULL;
995 dim = isl_map_get_dim(map);
997 d = isl_dim_size(dim, isl_dim_in);
998 dim = isl_dim_add(dim, isl_dim_in, 1);
999 dim = isl_dim_add(dim, isl_dim_out, 1);
1001 app = construct_power_components(isl_dim_copy(dim), map,
1004 diff = equate_parameter_to_length(dim, param);
1005 app = isl_map_intersect(app, diff);
1006 app = isl_map_project_out(app, isl_dim_in, d, 1);
1007 app = isl_map_project_out(app, isl_dim_out, d, 1);
1012 /* Compute the positive powers of "map", or an overapproximation.
1013 * The power is given by parameter "param". If the result is exact,
1014 * then *exact is set to 1.
1015 * If project is set, then we are actually interested in the transitive
1016 * closure, so we can use a more relaxed exactness check.
1018 static __isl_give isl_map *map_power(__isl_take isl_map *map, unsigned param,
1019 int *exact, int project)
1021 struct isl_map *app = NULL;
1026 map = isl_map_remove_empty_parts(map);
1030 if (isl_map_fast_is_empty(map))
1033 isl_assert(map->ctx, param < isl_map_dim(map, isl_dim_param), goto error);
1034 isl_assert(map->ctx,
1035 isl_map_dim(map, isl_dim_in) == isl_map_dim(map, isl_dim_out),
1038 app = construct_power(map, param, exact, project);
1048 /* Compute the positive powers of "map", or an overapproximation.
1049 * The power is given by parameter "param". If the result is exact,
1050 * then *exact is set to 1.
1052 __isl_give isl_map *isl_map_power(__isl_take isl_map *map, unsigned param,
1055 return map_power(map, param, exact, 0);
1058 /* Compute the transitive closure of "map", or an overapproximation.
1059 * If the result is exact, then *exact is set to 1.
1060 * Simply compute the powers of map and then project out the parameter
1061 * describing the power.
1063 __isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
1071 param = isl_map_dim(map, isl_dim_param);
1072 map = isl_map_add(map, isl_dim_param, 1);
1073 map = map_power(map, param, exact, 1);
1074 map = isl_map_project_out(map, isl_dim_param, param, 1);