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"
16 * The transitive closure implementation is based on the paper
17 * "Computing the Transitive Closure of a Union of Affine Integer
18 * Tuple Relations" by Anna Beletska, Denis Barthou, Wlodzimierz Bielecki and
22 /* Given a set of n offsets v_i (the rows of "steps"), construct a relation
23 * of the given dimension specification (Z^{n+1} -> Z^{n+1})
24 * that maps an element x to any element that can be reached
25 * by taking a non-negative number of steps along any of
26 * the extended offsets v'_i = [v_i 1].
29 * { [x] -> [y] : exists k_i >= 0, y = x + \sum_i k_i v'_i }
31 * For any element in this relation, the number of steps taken
32 * is equal to the difference in the final coordinates.
34 static __isl_give isl_map *path_along_steps(__isl_take isl_dim *dim,
35 __isl_keep isl_mat *steps)
38 struct isl_basic_map *path = NULL;
46 d = isl_dim_size(dim, isl_dim_in);
48 nparam = isl_dim_size(dim, isl_dim_param);
50 path = isl_basic_map_alloc_dim(isl_dim_copy(dim), n, d, n);
52 for (i = 0; i < n; ++i) {
53 k = isl_basic_map_alloc_div(path);
56 isl_assert(steps->ctx, i == k, goto error);
57 isl_int_set_si(path->div[k][0], 0);
60 for (i = 0; i < d; ++i) {
61 k = isl_basic_map_alloc_equality(path);
64 isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
65 isl_int_set_si(path->eq[k][1 + nparam + i], 1);
66 isl_int_set_si(path->eq[k][1 + nparam + d + i], -1);
68 for (j = 0; j < n; ++j)
69 isl_int_set_si(path->eq[k][1 + nparam + 2 * d + j], 1);
71 for (j = 0; j < n; ++j)
72 isl_int_set(path->eq[k][1 + nparam + 2 * d + j],
76 for (i = 0; i < n; ++i) {
77 k = isl_basic_map_alloc_inequality(path);
80 isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
81 isl_int_set_si(path->ineq[k][1 + nparam + 2 * d + i], 1);
86 path = isl_basic_map_simplify(path);
87 path = isl_basic_map_finalize(path);
88 return isl_map_from_basic_map(path);
91 isl_basic_map_free(path);
99 /* Return PURE_PARAM if only the coefficients of the parameters are non-zero.
100 * Return PURE_VAR if only the coefficients of the set variables are non-zero.
101 * Return IMPURE otherwise.
103 static int purity(__isl_keep isl_basic_set *bset, isl_int *c)
109 n_div = isl_basic_set_dim(bset, isl_dim_div);
110 d = isl_basic_set_dim(bset, isl_dim_set);
111 nparam = isl_basic_set_dim(bset, isl_dim_param);
113 if (isl_seq_first_non_zero(c + 1 + nparam + d, n_div) != -1)
115 if (isl_seq_first_non_zero(c + 1, nparam) == -1)
117 if (isl_seq_first_non_zero(c + 1 + nparam, d) == -1)
122 /* Given a set of offsets "delta", construct a relation of the
123 * given dimension specification (Z^{n+1} -> Z^{n+1}) that
124 * is an overapproximation of the relations that
125 * maps an element x to any element that can be reached
126 * by taking a non-negative number of steps along any of
127 * the elements in "delta".
128 * That is, construct an approximation of
130 * { [x] -> [y] : exists f \in \delta, k \in Z :
131 * y = x + k [f, 1] and k >= 0 }
133 * For any element in this relation, the number of steps taken
134 * is equal to the difference in the final coordinates.
136 * In particular, let delta be defined as
138 * \delta = [p] -> { [x] : A x + a >= and B p + b >= 0 and
139 * C x + C'p + c >= 0 }
141 * then the relation is constructed as
143 * { [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and
144 * A f + k a >= 0 and B p + b >= 0 and k >= 1 }
145 * union { [x] -> [x] }
147 * Existentially quantified variables in \delta are currently ignored.
148 * This is safe, but leads to an additional overapproximation.
150 static __isl_give isl_map *path_along_delta(__isl_take isl_dim *dim,
151 __isl_take isl_basic_set *delta)
153 isl_basic_map *path = NULL;
162 n_div = isl_basic_set_dim(delta, isl_dim_div);
163 d = isl_basic_set_dim(delta, isl_dim_set);
164 nparam = isl_basic_set_dim(delta, isl_dim_param);
165 path = isl_basic_map_alloc_dim(isl_dim_copy(dim), n_div + d + 1,
166 d + 1 + delta->n_eq, delta->n_ineq + 1);
167 off = 1 + nparam + 2 * (d + 1) + n_div;
169 for (i = 0; i < n_div + d + 1; ++i) {
170 k = isl_basic_map_alloc_div(path);
173 isl_int_set_si(path->div[k][0], 0);
176 for (i = 0; i < d + 1; ++i) {
177 k = isl_basic_map_alloc_equality(path);
180 isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
181 isl_int_set_si(path->eq[k][1 + nparam + i], 1);
182 isl_int_set_si(path->eq[k][1 + nparam + d + 1 + i], -1);
183 isl_int_set_si(path->eq[k][off + i], 1);
186 for (i = 0; i < delta->n_eq; ++i) {
187 int p = purity(delta, delta->eq[i]);
190 k = isl_basic_map_alloc_equality(path);
193 isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
195 isl_seq_cpy(path->eq[k] + off,
196 delta->eq[i] + 1 + nparam, d);
197 isl_int_set(path->eq[k][off + d], delta->eq[i][0]);
199 isl_seq_cpy(path->eq[k], delta->eq[i], 1 + nparam);
202 for (i = 0; i < delta->n_ineq; ++i) {
203 int p = purity(delta, delta->ineq[i]);
206 k = isl_basic_map_alloc_inequality(path);
209 isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
211 isl_seq_cpy(path->ineq[k] + off,
212 delta->ineq[i] + 1 + nparam, d);
213 isl_int_set(path->ineq[k][off + d], delta->ineq[i][0]);
215 isl_seq_cpy(path->ineq[k], delta->ineq[i], 1 + nparam);
218 k = isl_basic_map_alloc_inequality(path);
221 isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
222 isl_int_set_si(path->ineq[k][0], -1);
223 isl_int_set_si(path->ineq[k][off + d], 1);
225 isl_basic_set_free(delta);
226 path = isl_basic_map_finalize(path);
227 return isl_basic_map_union(path,
228 isl_basic_map_identity(isl_dim_domain(dim)));
231 isl_basic_set_free(delta);
232 isl_basic_map_free(path);
236 /* Given a dimenion specification Z^{n+1} -> Z^{n+1} and a parameter "param",
237 * construct a map that equates the parameter to the difference
238 * in the final coordinates and imposes that this difference is positive.
241 * { [x,x_s] -> [y,y_s] : k = y_s - x_s > 0 }
243 static __isl_give isl_map *equate_parameter_to_length(__isl_take isl_dim *dim,
246 struct isl_basic_map *bmap;
251 d = isl_dim_size(dim, isl_dim_in);
252 nparam = isl_dim_size(dim, isl_dim_param);
253 bmap = isl_basic_map_alloc_dim(dim, 0, 1, 1);
254 k = isl_basic_map_alloc_equality(bmap);
257 isl_seq_clr(bmap->eq[k], 1 + isl_basic_map_total_dim(bmap));
258 isl_int_set_si(bmap->eq[k][1 + param], -1);
259 isl_int_set_si(bmap->eq[k][1 + nparam + d - 1], -1);
260 isl_int_set_si(bmap->eq[k][1 + nparam + d + d - 1], 1);
262 k = isl_basic_map_alloc_inequality(bmap);
265 isl_seq_clr(bmap->ineq[k], 1 + isl_basic_map_total_dim(bmap));
266 isl_int_set_si(bmap->ineq[k][1 + param], 1);
267 isl_int_set_si(bmap->ineq[k][0], -1);
269 bmap = isl_basic_map_finalize(bmap);
270 return isl_map_from_basic_map(bmap);
272 isl_basic_map_free(bmap);
276 /* Check whether "path" is acyclic, where the last coordinates of domain
277 * and range of path encode the number of steps taken.
278 * That is, check whether
280 * { d | d = y - x and (x,y) in path }
282 * does not contain any element with positive last coordinate (positive length)
283 * and zero remaining coordinates (cycle).
285 static int is_acyclic(__isl_take isl_map *path)
290 struct isl_set *delta;
292 delta = isl_map_deltas(path);
293 dim = isl_set_dim(delta, isl_dim_set);
294 for (i = 0; i < dim; ++i) {
296 delta = isl_set_lower_bound_si(delta, isl_dim_set, i, 1);
298 delta = isl_set_fix_si(delta, isl_dim_set, i, 0);
301 acyclic = isl_set_is_empty(delta);
307 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
308 * and a dimension specification (Z^{n+1} -> Z^{n+1}),
309 * construct a map that is an overapproximation of the map
310 * that takes an element from the space D \times Z to another
311 * element from the same space, such that the first n coordinates of the
312 * difference between them is a sum of differences between images
313 * and pre-images in one of the R_i and such that the last coordinate
314 * is equal to the number of steps taken.
317 * \Delta_i = { y - x | (x, y) in R_i }
319 * then the constructed map is an overapproximation of
321 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
322 * d = (\sum_i k_i \delta_i, \sum_i k_i) }
324 * The elements of the singleton \Delta_i's are collected as the
325 * rows of the steps matrix. For all these \Delta_i's together,
326 * a single path is constructed.
327 * For each of the other \Delta_i's, we compute an overapproximation
328 * of the paths along elements of \Delta_i.
329 * Since each of these paths performs an addition, composition is
330 * symmetric and we can simply compose all resulting paths in any order.
332 static __isl_give isl_map *construct_extended_path(__isl_take isl_dim *dim,
333 __isl_keep isl_map *map, int *project)
335 struct isl_mat *steps = NULL;
336 struct isl_map *path = NULL;
340 d = isl_map_dim(map, isl_dim_in);
342 path = isl_map_identity(isl_dim_domain(isl_dim_copy(dim)));
344 steps = isl_mat_alloc(map->ctx, map->n, d);
349 for (i = 0; i < map->n; ++i) {
350 struct isl_basic_set *delta;
352 delta = isl_basic_map_deltas(isl_basic_map_copy(map->p[i]));
354 for (j = 0; j < d; ++j) {
357 fixed = isl_basic_set_fast_dim_is_fixed(delta, j,
360 isl_basic_set_free(delta);
369 path = isl_map_apply_range(path,
370 path_along_delta(isl_dim_copy(dim), delta));
372 isl_basic_set_free(delta);
379 path = isl_map_apply_range(path,
380 path_along_steps(isl_dim_copy(dim), steps));
383 if (project && *project) {
384 *project = is_acyclic(isl_map_copy(path));
399 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D,
400 * construct a map that is an overapproximation of the map
401 * that takes an element from the space D to another
402 * element from the same space, such that the difference between
403 * them is a strictly positive sum of differences between images
404 * and pre-images in one of the R_i.
405 * The number of differences in the sum is equated to parameter "param".
408 * \Delta_i = { y - x | (x, y) in R_i }
410 * then the constructed map is an overapproximation of
412 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
413 * d = \sum_i k_i \delta_i and k = \sum_i k_i > 0 }
415 * We first construct an extended mapping with an extra coordinate
416 * that indicates the number of steps taken. In particular,
417 * the difference in the last coordinate is equal to the number
418 * of steps taken to move from a domain element to the corresponding
420 * In the final step, this difference is equated to the parameter "param"
421 * and made positive. The extra coordinates are subsequently projected out.
423 static __isl_give isl_map *construct_path(__isl_keep isl_map *map,
424 unsigned param, int *project)
426 struct isl_map *path = NULL;
427 struct isl_map *diff;
428 struct isl_dim *dim = NULL;
434 dim = isl_map_get_dim(map);
436 d = isl_dim_size(dim, isl_dim_in);
437 dim = isl_dim_add(dim, isl_dim_in, 1);
438 dim = isl_dim_add(dim, isl_dim_out, 1);
440 path = construct_extended_path(isl_dim_copy(dim), map, project);
442 diff = equate_parameter_to_length(dim, param);
443 path = isl_map_intersect(path, diff);
444 path = isl_map_project_out(path, isl_dim_in, d, 1);
445 path = isl_map_project_out(path, isl_dim_out, d, 1);
450 /* Shift variable at position "pos" up by one.
451 * That is, replace the corresponding variable v by v - 1.
453 static __isl_give isl_basic_map *basic_map_shift_pos(
454 __isl_take isl_basic_map *bmap, unsigned pos)
458 bmap = isl_basic_map_cow(bmap);
462 for (i = 0; i < bmap->n_eq; ++i)
463 isl_int_sub(bmap->eq[i][0], bmap->eq[i][0], bmap->eq[i][pos]);
465 for (i = 0; i < bmap->n_ineq; ++i)
466 isl_int_sub(bmap->ineq[i][0],
467 bmap->ineq[i][0], bmap->ineq[i][pos]);
469 for (i = 0; i < bmap->n_div; ++i) {
470 if (isl_int_is_zero(bmap->div[i][0]))
472 isl_int_sub(bmap->div[i][1],
473 bmap->div[i][1], bmap->div[i][1 + pos]);
479 /* Shift variable at position "pos" up by one.
480 * That is, replace the corresponding variable v by v - 1.
482 static __isl_give isl_map *map_shift_pos(__isl_take isl_map *map, unsigned pos)
486 map = isl_map_cow(map);
490 for (i = 0; i < map->n; ++i) {
491 map->p[i] = basic_map_shift_pos(map->p[i], pos);
495 ISL_F_CLR(map, ISL_MAP_NORMALIZED);
502 /* Check whether the overapproximation of the power of "map" is exactly
503 * the power of "map". Let R be "map" and A_k the overapproximation.
504 * The approximation is exact if
507 * A_k = A_{k-1} \circ R k >= 2
509 * Since A_k is known to be an overapproximation, we only need to check
512 * A_k \subset A_{k-1} \circ R k >= 2
515 static int check_power_exactness(__isl_take isl_map *map,
516 __isl_take isl_map *app, unsigned param)
522 app_1 = isl_map_fix_si(isl_map_copy(app), isl_dim_param, param, 1);
524 exact = isl_map_is_subset(app_1, map);
527 if (!exact || exact < 0) {
533 app_2 = isl_map_lower_bound_si(isl_map_copy(app),
534 isl_dim_param, param, 2);
535 app_1 = map_shift_pos(app, 1 + param);
536 app_1 = isl_map_apply_range(map, app_1);
538 exact = isl_map_is_subset(app_2, app_1);
546 /* Check whether the overapproximation of the power of "map" is exactly
547 * the power of "map", possibly after projecting out the power (if "project"
550 * If "project" is set and if "steps" can only result in acyclic paths,
553 * A = R \cup (A \circ R)
555 * where A is the overapproximation with the power projected out, i.e.,
556 * an overapproximation of the transitive closure.
557 * More specifically, since A is known to be an overapproximation, we check
559 * A \subset R \cup (A \circ R)
561 * Otherwise, we check if the power is exact.
563 static int check_exactness(__isl_take isl_map *map, __isl_take isl_map *app,
564 unsigned param, int project)
570 return check_power_exactness(map, app, param);
572 map = isl_map_project_out(map, isl_dim_param, param, 1);
573 app = isl_map_project_out(app, isl_dim_param, param, 1);
575 test = isl_map_apply_range(isl_map_copy(map), isl_map_copy(app));
576 test = isl_map_union(test, isl_map_copy(map));
578 exact = isl_map_is_subset(app, test);
592 /* Compute the positive powers of "map", or an overapproximation.
593 * The power is given by parameter "param". If the result is exact,
594 * then *exact is set to 1.
595 * If project is set, then we are actually interested in the transitive
596 * closure, so we can use a more relaxed exactness check.
598 static __isl_give isl_map *map_power(__isl_take isl_map *map, unsigned param,
599 int *exact, int project)
601 struct isl_set *domain = NULL;
602 struct isl_set *range = NULL;
603 struct isl_map *app = NULL;
604 struct isl_map *path = NULL;
609 map = isl_map_remove_empty_parts(map);
613 if (isl_map_fast_is_empty(map))
616 isl_assert(map->ctx, param < isl_map_dim(map, isl_dim_param), goto error);
618 isl_map_dim(map, isl_dim_in) == isl_map_dim(map, isl_dim_out),
621 domain = isl_map_domain(isl_map_copy(map));
622 domain = isl_set_coalesce(domain);
623 range = isl_map_range(isl_map_copy(map));
624 range = isl_set_coalesce(range);
625 app = isl_map_from_domain_and_range(isl_set_copy(domain),
626 isl_set_copy(range));
628 path = construct_path(map, param, exact ? &project : NULL);
629 app = isl_map_intersect(app, isl_map_copy(path));
632 (*exact = check_exactness(isl_map_copy(map), isl_map_copy(app),
633 param, project)) < 0)
636 isl_set_free(domain);
642 isl_set_free(domain);
650 /* Compute the positive powers of "map", or an overapproximation.
651 * The power is given by parameter "param". If the result is exact,
652 * then *exact is set to 1.
654 __isl_give isl_map *isl_map_power(__isl_take isl_map *map, unsigned param,
657 return map_power(map, param, exact, 0);
660 /* Compute the transitive closure of "map", or an overapproximation.
661 * If the result is exact, then *exact is set to 1.
662 * Simply compute the powers of map and then project out the parameter
663 * describing the power.
665 __isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
673 param = isl_map_dim(map, isl_dim_param);
674 map = isl_map_add(map, isl_dim_param, 1);
675 map = map_power(map, param, exact, 1);
676 map = isl_map_project_out(map, isl_dim_param, param, 1);