add isl_printer
[isl.git] / isl_transitive_closure.c
blob879e0ad3facc0c0df429bbb36d8af815a3cc7a5a
1 /*
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,
8 * 91893 Orsay, France
9 */
11 #include "isl_map.h"
12 #include "isl_map_private.h"
13 #include "isl_seq.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)
24 struct isl_dim *dim;
25 struct isl_basic_map *bmap;
26 unsigned d;
27 unsigned nparam;
28 int k;
29 isl_int *c;
31 if (!map)
32 return NULL;
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);
38 if (exactly) {
39 k = isl_basic_map_alloc_equality(bmap);
40 c = bmap->eq[k];
41 } else {
42 k = isl_basic_map_alloc_inequality(bmap);
43 c = bmap->ineq[k];
45 if (k < 0)
46 goto error;
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));
55 return map;
56 error:
57 isl_basic_map_free(bmap);
58 isl_map_free(map);
59 return NULL;
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
66 * A_1 = R
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
71 * A_1 \subset R
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
77 * one.
79 static int check_power_exactness(__isl_take isl_map *map,
80 __isl_take isl_map *app)
82 int exact;
83 isl_map *app_1;
84 isl_map *app_2;
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);
93 isl_map_free(app_1);
95 if (!exact || exact < 0) {
96 isl_map_free(app);
97 isl_map_free(map);
98 return exact;
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);
107 isl_map_free(app_1);
108 isl_map_free(app_2);
110 return exact;
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"
115 * is set).
117 * If "project" is set and if "steps" can only result in acyclic paths,
118 * then we check
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,
135 int project)
137 isl_map *test;
138 int exact;
139 unsigned d;
141 if (!project)
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);
154 isl_map_free(app);
155 isl_map_free(test);
157 isl_map_free(map);
159 return exact;
160 error:
161 isl_map_free(app);
162 isl_map_free(map);
163 return -1;
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
170 * Albert Cohen.
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].
178 * That is, construct
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)
188 int i, j, k;
189 struct isl_basic_map *path = NULL;
190 unsigned d;
191 unsigned n;
192 unsigned nparam;
194 if (!dim || !steps)
195 goto error;
197 d = isl_dim_size(dim, isl_dim_in);
198 n = steps->n_row;
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);
205 if (k < 0)
206 goto error;
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);
213 if (k < 0)
214 goto error;
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);
218 if (i == d - 1)
219 for (j = 0; j < n; ++j)
220 isl_int_set_si(path->eq[k][1 + nparam + 2 * d + j], 1);
221 else
222 for (j = 0; j < n; ++j)
223 isl_int_set(path->eq[k][1 + nparam + 2 * d + j],
224 steps->row[j][i]);
227 for (i = 0; i < n; ++i) {
228 k = isl_basic_map_alloc_inequality(path);
229 if (k < 0)
230 goto error;
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);
235 isl_dim_free(dim);
237 path = isl_basic_map_simplify(path);
238 path = isl_basic_map_finalize(path);
239 return isl_map_from_basic_map(path);
240 error:
241 isl_dim_free(dim);
242 isl_basic_map_free(path);
243 return NULL;
246 #define IMPURE 0
247 #define PURE_PARAM 1
248 #define PURE_VAR 2
249 #define MIXED 3
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)
260 unsigned d;
261 unsigned n_div;
262 unsigned nparam;
263 int k;
264 int empty;
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)
271 return IMPURE;
272 if (isl_seq_first_non_zero(c + 1, nparam) == -1)
273 return PURE_VAR;
274 if (isl_seq_first_non_zero(c + 1 + nparam, d) == -1)
275 return PURE_PARAM;
276 if (eq)
277 return IMPURE;
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);
283 if (k < 0)
284 goto error;
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;
292 error:
293 isl_basic_set_free(bset);
294 return -1;
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;
333 unsigned d;
334 unsigned n_div;
335 unsigned nparam;
336 unsigned off;
337 int i, k;
339 if (!delta)
340 goto error;
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);
350 if (k < 0)
351 goto error;
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);
357 if (k < 0)
358 goto error;
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);
367 if (p < 0)
368 goto error;
369 if (p == IMPURE)
370 continue;
371 k = isl_basic_map_alloc_equality(path);
372 if (k < 0)
373 goto error;
374 isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
375 if (p == PURE_VAR) {
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]);
379 } else
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);
385 if (p < 0)
386 goto error;
387 if (p == IMPURE)
388 continue;
389 k = isl_basic_map_alloc_inequality(path);
390 if (k < 0)
391 goto error;
392 isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
393 if (p == PURE_VAR) {
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);
399 } else {
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);
407 if (k < 0)
408 goto error;
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)));
417 error:
418 isl_dim_free(dim);
419 isl_basic_set_free(delta);
420 isl_basic_map_free(path);
421 return NULL;
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.
427 * That is, construct
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,
432 unsigned param)
434 struct isl_basic_map *bmap;
435 unsigned d;
436 unsigned nparam;
437 int k;
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);
443 if (k < 0)
444 goto error;
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);
451 if (k < 0)
452 goto error;
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);
459 error:
460 isl_basic_map_free(bmap);
461 return NULL;
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)
475 int i;
476 int acyclic;
477 unsigned dim;
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) {
483 if (i == dim -1)
484 delta = isl_set_lower_bound_si(delta, isl_dim_set, i, 1);
485 else
486 delta = isl_set_fix_si(delta, isl_dim_set, i, 0);
489 acyclic = isl_set_is_empty(delta);
490 isl_set_free(delta);
492 return acyclic;
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.
503 * That is, let
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;
525 unsigned d;
526 int i, j, n;
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);
533 if (!steps)
534 goto error;
536 n = 0;
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) {
543 int fixed;
545 fixed = isl_basic_set_fast_dim_is_fixed(delta, j,
546 &steps->row[n][j]);
547 if (fixed < 0) {
548 isl_basic_set_free(delta);
549 goto error;
551 if (!fixed)
552 break;
556 if (j < d) {
557 path = isl_map_apply_range(path,
558 path_along_delta(isl_dim_copy(dim), delta));
559 } else {
560 isl_basic_set_free(delta);
561 ++n;
565 if (n > 0) {
566 steps->n_row = n;
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));
573 if (*project < 0)
574 goto error;
577 isl_dim_free(dim);
578 isl_mat_free(steps);
579 return path;
580 error:
581 isl_dim_free(dim);
582 isl_mat_free(steps);
583 isl_map_free(path);
584 return NULL;
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.
596 * That is, let
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
605 * { (x) -> (x) }
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),
629 project)) < 0)
630 goto error;
632 return isl_map_union(app, isl_map_identity(isl_dim_domain(dim)));
633 error:
634 isl_dim_free(dim);
635 isl_map_free(app);
636 return NULL;
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 {
646 int index;
647 int min_index;
648 int on_stack;
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 {
661 int len;
662 struct basic_map_sort_node *node;
663 int *stack;
664 int sp;
665 int index;
666 int *order;
667 int op;
670 static void basic_map_sort_free(struct basic_map_sort *s)
672 if (!s)
673 return;
674 free(s->node);
675 free(s->stack);
676 free(s->order);
677 free(s);
680 static struct basic_map_sort *basic_map_sort_alloc(struct isl_ctx *ctx, int len)
682 struct basic_map_sort *s;
683 int i;
685 s = isl_calloc_type(ctx, struct basic_map_sort);
686 if (!s)
687 return NULL;
688 s->len = len;
689 s->node = isl_alloc_array(ctx, struct basic_map_sort_node, len);
690 if (!s->node)
691 goto error;
692 for (i = 0; i < len; ++i)
693 s->node[i].index = -1;
694 s->stack = isl_alloc_array(ctx, int, len);
695 if (!s->stack)
696 goto error;
697 s->order = isl_alloc_array(ctx, int, 2 * len);
698 if (!s->order)
699 goto error;
701 s->sp = 0;
702 s->index = 0;
703 s->op = 0;
705 return s;
706 error:
707 basic_map_sort_free(s);
708 return NULL;
711 /* Check whether in the computation of the transitive closure
712 * "bmap1" (R_1) should follow (or be part of the same component as)
713 * "bmap2" (R_2).
715 * That is check whether
717 * R_1 \circ R_2
719 * is a subset of
721 * R_2 \circ R_1
723 * If so, then there is no reason for R_1 to immediately follow R_2
724 * in any path.
726 static int basic_map_follows(__isl_keep isl_basic_map *bmap1,
727 __isl_keep isl_basic_map *bmap2)
729 struct isl_map *map12 = NULL;
730 struct isl_map *map21 = NULL;
731 int subset;
733 map21 = isl_map_from_basic_map(
734 isl_basic_map_apply_range(
735 isl_basic_map_copy(bmap2),
736 isl_basic_map_copy(bmap1)));
737 subset = isl_map_is_empty(map21);
738 if (subset < 0)
739 goto error;
740 if (subset) {
741 isl_map_free(map21);
742 return 0;
745 map12 = isl_map_from_basic_map(
746 isl_basic_map_apply_range(
747 isl_basic_map_copy(bmap1),
748 isl_basic_map_copy(bmap2)));
750 subset = isl_map_is_subset(map21, map12);
752 isl_map_free(map12);
753 isl_map_free(map21);
755 return subset < 0 ? -1 : !subset;
756 error:
757 isl_map_free(map21);
758 return -1;
761 /* Perform Tarjan's algorithm for computing the strongly connected components
762 * in the graph with the disjuncts of "map" as vertices and with an
763 * edge between any pair of disjuncts such that the first has
764 * to be applied after the second.
766 static int power_components_tarjan(struct basic_map_sort *s,
767 __isl_keep isl_map *map, int i)
769 int j;
771 s->node[i].index = s->index;
772 s->node[i].min_index = s->index;
773 s->node[i].on_stack = 1;
774 s->index++;
775 s->stack[s->sp++] = i;
777 for (j = s->len - 1; j >= 0; --j) {
778 int f;
780 if (j == i)
781 continue;
782 if (s->node[j].index >= 0 &&
783 (!s->node[j].on_stack ||
784 s->node[j].index > s->node[i].min_index))
785 continue;
787 f = basic_map_follows(map->p[i], map->p[j]);
788 if (f < 0)
789 return -1;
790 if (!f)
791 continue;
793 if (s->node[j].index < 0) {
794 power_components_tarjan(s, map, j);
795 if (s->node[j].min_index < s->node[i].min_index)
796 s->node[i].min_index = s->node[j].min_index;
797 } else if (s->node[j].index < s->node[i].min_index)
798 s->node[i].min_index = s->node[j].index;
801 if (s->node[i].index != s->node[i].min_index)
802 return 0;
804 do {
805 j = s->stack[--s->sp];
806 s->node[j].on_stack = 0;
807 s->order[s->op++] = j;
808 } while (j != i);
809 s->order[s->op++] = -1;
811 return 0;
814 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
815 * and a dimension specification (Z^{n+1} -> Z^{n+1}),
816 * construct a map that is the union of the identity map and
817 * an overapproximation of the map
818 * that takes an element from the dom R \times Z to an
819 * element from ran R \times Z, such that the first n coordinates of the
820 * difference between them is a sum of differences between images
821 * and pre-images in one of the R_i and such that the last coordinate
822 * is equal to the number of steps taken.
823 * That is, let
825 * \Delta_i = { y - x | (x, y) in R_i }
827 * then the constructed map is an overapproximation of
829 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
830 * d = (\sum_i k_i \delta_i, \sum_i k_i) and
831 * x in dom R and x + d in ran R } union
832 * { (x) -> (x) }
834 * We first split the map into strongly connected components, perform
835 * the above on each component and the join the results in the correct
836 * order. The power of each of the components needs to be extended
837 * with the identity map because a path in the global result need
838 * not go through every component.
839 * The final result will then also contain the identity map, but
840 * this part will be removed when the length of the path is forced
841 * to be strictly positive.
843 static __isl_give isl_map *construct_power_components(__isl_take isl_dim *dim,
844 __isl_keep isl_map *map, int *exact, int project)
846 int i, n;
847 struct isl_map *path = NULL;
848 struct basic_map_sort *s = NULL;
850 if (!map)
851 goto error;
852 if (map->n <= 1)
853 return construct_component(dim, map, exact, project);
855 s = basic_map_sort_alloc(map->ctx, map->n);
856 if (!s)
857 goto error;
858 for (i = map->n - 1; i >= 0; --i) {
859 if (s->node[i].index >= 0)
860 continue;
861 if (power_components_tarjan(s, map, i) < 0)
862 goto error;
865 i = 0;
866 n = map->n;
867 path = isl_map_identity(isl_dim_domain(isl_dim_copy(dim)));
868 while (n) {
869 struct isl_map *comp;
870 comp = isl_map_alloc_dim(isl_map_get_dim(map), n, 0);
871 while (s->order[i] != -1) {
872 comp = isl_map_add_basic_map(comp,
873 isl_basic_map_copy(map->p[s->order[i]]));
874 --n;
875 ++i;
877 path = isl_map_apply_range(path,
878 construct_component(isl_dim_copy(dim), comp,
879 exact, project));
880 isl_map_free(comp);
881 ++i;
884 basic_map_sort_free(s);
885 isl_dim_free(dim);
887 return path;
888 error:
889 basic_map_sort_free(s);
890 isl_dim_free(dim);
891 return NULL;
894 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D,
895 * construct a map that is an overapproximation of the map
896 * that takes an element from the space D to another
897 * element from the same space, such that the difference between
898 * them is a strictly positive sum of differences between images
899 * and pre-images in one of the R_i.
900 * The number of differences in the sum is equated to parameter "param".
901 * That is, let
903 * \Delta_i = { y - x | (x, y) in R_i }
905 * then the constructed map is an overapproximation of
907 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
908 * d = \sum_i k_i \delta_i and k = \sum_i k_i > 0 }
910 * We first construct an extended mapping with an extra coordinate
911 * that indicates the number of steps taken. In particular,
912 * the difference in the last coordinate is equal to the number
913 * of steps taken to move from a domain element to the corresponding
914 * image element(s).
915 * In the final step, this difference is equated to the parameter "param"
916 * and made positive. The extra coordinates are subsequently projected out.
918 static __isl_give isl_map *construct_power(__isl_keep isl_map *map,
919 unsigned param, int *exact, int project)
921 struct isl_map *app = NULL;
922 struct isl_map *diff;
923 struct isl_dim *dim = NULL;
924 unsigned d;
926 if (!map)
927 return NULL;
929 dim = isl_map_get_dim(map);
931 d = isl_dim_size(dim, isl_dim_in);
932 dim = isl_dim_add(dim, isl_dim_in, 1);
933 dim = isl_dim_add(dim, isl_dim_out, 1);
935 app = construct_power_components(isl_dim_copy(dim), map,
936 exact, project);
938 diff = equate_parameter_to_length(dim, param);
939 app = isl_map_intersect(app, diff);
940 app = isl_map_project_out(app, isl_dim_in, d, 1);
941 app = isl_map_project_out(app, isl_dim_out, d, 1);
943 return app;
946 /* Compute the positive powers of "map", or an overapproximation.
947 * The power is given by parameter "param". If the result is exact,
948 * then *exact is set to 1.
949 * If project is set, then we are actually interested in the transitive
950 * closure, so we can use a more relaxed exactness check.
952 static __isl_give isl_map *map_power(__isl_take isl_map *map, unsigned param,
953 int *exact, int project)
955 struct isl_map *app = NULL;
957 if (exact)
958 *exact = 1;
960 map = isl_map_remove_empty_parts(map);
961 if (!map)
962 return NULL;
964 if (isl_map_fast_is_empty(map))
965 return map;
967 isl_assert(map->ctx, param < isl_map_dim(map, isl_dim_param), goto error);
968 isl_assert(map->ctx,
969 isl_map_dim(map, isl_dim_in) == isl_map_dim(map, isl_dim_out),
970 goto error);
972 app = construct_power(map, param, exact, project);
974 isl_map_free(map);
975 return app;
976 error:
977 isl_map_free(map);
978 isl_map_free(app);
979 return NULL;
982 /* Compute the positive powers of "map", or an overapproximation.
983 * The power is given by parameter "param". If the result is exact,
984 * then *exact is set to 1.
986 __isl_give isl_map *isl_map_power(__isl_take isl_map *map, unsigned param,
987 int *exact)
989 return map_power(map, param, exact, 0);
992 /* Compute the transitive closure of "map", or an overapproximation.
993 * If the result is exact, then *exact is set to 1.
994 * Simply compute the powers of map and then project out the parameter
995 * describing the power.
997 __isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
998 int *exact)
1000 unsigned param;
1002 if (!map)
1003 goto error;
1005 param = isl_map_dim(map, isl_dim_param);
1006 map = isl_map_add(map, isl_dim_param, 1);
1007 map = map_power(map, param, exact, 1);
1008 map = isl_map_project_out(map, isl_dim_param, param, 1);
1010 return map;
1011 error:
1012 isl_map_free(map);
1013 return NULL;