Line data Source code
1 : /*
2 : * Copyright 2010 INRIA Saclay
3 : *
4 : * Use of this software is governed by the MIT license
5 : *
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 : */
10 :
11 : #include <isl_ctx_private.h>
12 : #include <isl_map_private.h>
13 : #include <isl/map.h>
14 : #include <isl_seq.h>
15 : #include <isl_space_private.h>
16 : #include <isl_lp_private.h>
17 : #include <isl/union_map.h>
18 : #include <isl_mat_private.h>
19 : #include <isl_vec_private.h>
20 : #include <isl_options_private.h>
21 : #include <isl_tarjan.h>
22 :
23 0 : int isl_map_is_transitively_closed(__isl_keep isl_map *map)
24 : {
25 : isl_map *map2;
26 : int closed;
27 :
28 0 : map2 = isl_map_apply_range(isl_map_copy(map), isl_map_copy(map));
29 0 : closed = isl_map_is_subset(map2, map);
30 0 : isl_map_free(map2);
31 :
32 0 : return closed;
33 : }
34 :
35 0 : int isl_union_map_is_transitively_closed(__isl_keep isl_union_map *umap)
36 : {
37 : isl_union_map *umap2;
38 : int closed;
39 :
40 0 : umap2 = isl_union_map_apply_range(isl_union_map_copy(umap),
41 : isl_union_map_copy(umap));
42 0 : closed = isl_union_map_is_subset(umap2, umap);
43 0 : isl_union_map_free(umap2);
44 :
45 0 : return closed;
46 : }
47 :
48 : /* Given a map that represents a path with the length of the path
49 : * encoded as the difference between the last output coordindate
50 : * and the last input coordinate, set this length to either
51 : * exactly "length" (if "exactly" is set) or at least "length"
52 : * (if "exactly" is not set).
53 : */
54 0 : static __isl_give isl_map *set_path_length(__isl_take isl_map *map,
55 : int exactly, int length)
56 : {
57 : isl_space *dim;
58 : struct isl_basic_map *bmap;
59 : unsigned d;
60 : unsigned nparam;
61 : int k;
62 : isl_int *c;
63 :
64 0 : if (!map)
65 0 : return NULL;
66 :
67 0 : dim = isl_map_get_space(map);
68 0 : d = isl_space_dim(dim, isl_dim_in);
69 0 : nparam = isl_space_dim(dim, isl_dim_param);
70 0 : bmap = isl_basic_map_alloc_space(dim, 0, 1, 1);
71 0 : if (exactly) {
72 0 : k = isl_basic_map_alloc_equality(bmap);
73 0 : if (k < 0)
74 0 : goto error;
75 0 : c = bmap->eq[k];
76 : } else {
77 0 : k = isl_basic_map_alloc_inequality(bmap);
78 0 : if (k < 0)
79 0 : goto error;
80 0 : c = bmap->ineq[k];
81 : }
82 0 : isl_seq_clr(c, 1 + isl_basic_map_total_dim(bmap));
83 0 : isl_int_set_si(c[0], -length);
84 0 : isl_int_set_si(c[1 + nparam + d - 1], -1);
85 0 : isl_int_set_si(c[1 + nparam + d + d - 1], 1);
86 :
87 0 : bmap = isl_basic_map_finalize(bmap);
88 0 : map = isl_map_intersect(map, isl_map_from_basic_map(bmap));
89 :
90 0 : return map;
91 : error:
92 0 : isl_basic_map_free(bmap);
93 0 : isl_map_free(map);
94 0 : return NULL;
95 : }
96 :
97 : /* Check whether the overapproximation of the power of "map" is exactly
98 : * the power of "map". Let R be "map" and A_k the overapproximation.
99 : * The approximation is exact if
100 : *
101 : * A_1 = R
102 : * A_k = A_{k-1} \circ R k >= 2
103 : *
104 : * Since A_k is known to be an overapproximation, we only need to check
105 : *
106 : * A_1 \subset R
107 : * A_k \subset A_{k-1} \circ R k >= 2
108 : *
109 : * In practice, "app" has an extra input and output coordinate
110 : * to encode the length of the path. So, we first need to add
111 : * this coordinate to "map" and set the length of the path to
112 : * one.
113 : */
114 0 : static int check_power_exactness(__isl_take isl_map *map,
115 : __isl_take isl_map *app)
116 : {
117 : int exact;
118 : isl_map *app_1;
119 : isl_map *app_2;
120 :
121 0 : map = isl_map_add_dims(map, isl_dim_in, 1);
122 0 : map = isl_map_add_dims(map, isl_dim_out, 1);
123 0 : map = set_path_length(map, 1, 1);
124 :
125 0 : app_1 = set_path_length(isl_map_copy(app), 1, 1);
126 :
127 0 : exact = isl_map_is_subset(app_1, map);
128 0 : isl_map_free(app_1);
129 :
130 0 : if (!exact || exact < 0) {
131 0 : isl_map_free(app);
132 0 : isl_map_free(map);
133 0 : return exact;
134 : }
135 :
136 0 : app_1 = set_path_length(isl_map_copy(app), 0, 1);
137 0 : app_2 = set_path_length(app, 0, 2);
138 0 : app_1 = isl_map_apply_range(map, app_1);
139 :
140 0 : exact = isl_map_is_subset(app_2, app_1);
141 :
142 0 : isl_map_free(app_1);
143 0 : isl_map_free(app_2);
144 :
145 0 : return exact;
146 : }
147 :
148 : /* Check whether the overapproximation of the power of "map" is exactly
149 : * the power of "map", possibly after projecting out the power (if "project"
150 : * is set).
151 : *
152 : * If "project" is set and if "steps" can only result in acyclic paths,
153 : * then we check
154 : *
155 : * A = R \cup (A \circ R)
156 : *
157 : * where A is the overapproximation with the power projected out, i.e.,
158 : * an overapproximation of the transitive closure.
159 : * More specifically, since A is known to be an overapproximation, we check
160 : *
161 : * A \subset R \cup (A \circ R)
162 : *
163 : * Otherwise, we check if the power is exact.
164 : *
165 : * Note that "app" has an extra input and output coordinate to encode
166 : * the length of the part. If we are only interested in the transitive
167 : * closure, then we can simply project out these coordinates first.
168 : */
169 0 : static int check_exactness(__isl_take isl_map *map, __isl_take isl_map *app,
170 : int project)
171 : {
172 : isl_map *test;
173 : int exact;
174 : unsigned d;
175 :
176 0 : if (!project)
177 0 : return check_power_exactness(map, app);
178 :
179 0 : d = isl_map_dim(map, isl_dim_in);
180 0 : app = set_path_length(app, 0, 1);
181 0 : app = isl_map_project_out(app, isl_dim_in, d, 1);
182 0 : app = isl_map_project_out(app, isl_dim_out, d, 1);
183 :
184 0 : app = isl_map_reset_space(app, isl_map_get_space(map));
185 :
186 0 : test = isl_map_apply_range(isl_map_copy(map), isl_map_copy(app));
187 0 : test = isl_map_union(test, isl_map_copy(map));
188 :
189 0 : exact = isl_map_is_subset(app, test);
190 :
191 0 : isl_map_free(app);
192 0 : isl_map_free(test);
193 :
194 0 : isl_map_free(map);
195 :
196 0 : return exact;
197 : }
198 :
199 : /*
200 : * The transitive closure implementation is based on the paper
201 : * "Computing the Transitive Closure of a Union of Affine Integer
202 : * Tuple Relations" by Anna Beletska, Denis Barthou, Wlodzimierz Bielecki and
203 : * Albert Cohen.
204 : */
205 :
206 : /* Given a set of n offsets v_i (the rows of "steps"), construct a relation
207 : * of the given dimension specification (Z^{n+1} -> Z^{n+1})
208 : * that maps an element x to any element that can be reached
209 : * by taking a non-negative number of steps along any of
210 : * the extended offsets v'_i = [v_i 1].
211 : * That is, construct
212 : *
213 : * { [x] -> [y] : exists k_i >= 0, y = x + \sum_i k_i v'_i }
214 : *
215 : * For any element in this relation, the number of steps taken
216 : * is equal to the difference in the final coordinates.
217 : */
218 0 : static __isl_give isl_map *path_along_steps(__isl_take isl_space *dim,
219 : __isl_keep isl_mat *steps)
220 : {
221 : int i, j, k;
222 0 : struct isl_basic_map *path = NULL;
223 : unsigned d;
224 : unsigned n;
225 : unsigned nparam;
226 :
227 0 : if (!dim || !steps)
228 : goto error;
229 :
230 0 : d = isl_space_dim(dim, isl_dim_in);
231 0 : n = steps->n_row;
232 0 : nparam = isl_space_dim(dim, isl_dim_param);
233 :
234 0 : path = isl_basic_map_alloc_space(isl_space_copy(dim), n, d, n);
235 :
236 0 : for (i = 0; i < n; ++i) {
237 0 : k = isl_basic_map_alloc_div(path);
238 0 : if (k < 0)
239 0 : goto error;
240 0 : isl_assert(steps->ctx, i == k, goto error);
241 0 : isl_int_set_si(path->div[k][0], 0);
242 : }
243 :
244 0 : for (i = 0; i < d; ++i) {
245 0 : k = isl_basic_map_alloc_equality(path);
246 0 : if (k < 0)
247 0 : goto error;
248 0 : isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
249 0 : isl_int_set_si(path->eq[k][1 + nparam + i], 1);
250 0 : isl_int_set_si(path->eq[k][1 + nparam + d + i], -1);
251 0 : if (i == d - 1)
252 0 : for (j = 0; j < n; ++j)
253 0 : isl_int_set_si(path->eq[k][1 + nparam + 2 * d + j], 1);
254 : else
255 0 : for (j = 0; j < n; ++j)
256 0 : isl_int_set(path->eq[k][1 + nparam + 2 * d + j],
257 : steps->row[j][i]);
258 : }
259 :
260 0 : for (i = 0; i < n; ++i) {
261 0 : k = isl_basic_map_alloc_inequality(path);
262 0 : if (k < 0)
263 0 : goto error;
264 0 : isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
265 0 : isl_int_set_si(path->ineq[k][1 + nparam + 2 * d + i], 1);
266 : }
267 :
268 0 : isl_space_free(dim);
269 :
270 0 : path = isl_basic_map_simplify(path);
271 0 : path = isl_basic_map_finalize(path);
272 0 : return isl_map_from_basic_map(path);
273 : error:
274 0 : isl_space_free(dim);
275 0 : isl_basic_map_free(path);
276 0 : return NULL;
277 : }
278 :
279 : #define IMPURE 0
280 : #define PURE_PARAM 1
281 : #define PURE_VAR 2
282 : #define MIXED 3
283 :
284 : /* Check whether the parametric constant term of constraint c is never
285 : * positive in "bset".
286 : */
287 0 : static isl_bool parametric_constant_never_positive(
288 : __isl_keep isl_basic_set *bset, isl_int *c, int *div_purity)
289 : {
290 : unsigned d;
291 : unsigned n_div;
292 : unsigned nparam;
293 : int i;
294 : int k;
295 : isl_bool empty;
296 :
297 0 : n_div = isl_basic_set_dim(bset, isl_dim_div);
298 0 : d = isl_basic_set_dim(bset, isl_dim_set);
299 0 : nparam = isl_basic_set_dim(bset, isl_dim_param);
300 :
301 0 : bset = isl_basic_set_copy(bset);
302 0 : bset = isl_basic_set_cow(bset);
303 0 : bset = isl_basic_set_extend_constraints(bset, 0, 1);
304 0 : k = isl_basic_set_alloc_inequality(bset);
305 0 : if (k < 0)
306 0 : goto error;
307 0 : isl_seq_clr(bset->ineq[k], 1 + isl_basic_set_total_dim(bset));
308 0 : isl_seq_cpy(bset->ineq[k], c, 1 + nparam);
309 0 : for (i = 0; i < n_div; ++i) {
310 0 : if (div_purity[i] != PURE_PARAM)
311 0 : continue;
312 0 : isl_int_set(bset->ineq[k][1 + nparam + d + i],
313 : c[1 + nparam + d + i]);
314 : }
315 0 : isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
316 0 : empty = isl_basic_set_is_empty(bset);
317 0 : isl_basic_set_free(bset);
318 :
319 0 : return empty;
320 : error:
321 0 : isl_basic_set_free(bset);
322 0 : return isl_bool_error;
323 : }
324 :
325 : /* Return PURE_PARAM if only the coefficients of the parameters are non-zero.
326 : * Return PURE_VAR if only the coefficients of the set variables are non-zero.
327 : * Return MIXED if only the coefficients of the parameters and the set
328 : * variables are non-zero and if moreover the parametric constant
329 : * can never attain positive values.
330 : * Return IMPURE otherwise.
331 : */
332 0 : static int purity(__isl_keep isl_basic_set *bset, isl_int *c, int *div_purity,
333 : int eq)
334 : {
335 : unsigned d;
336 : unsigned n_div;
337 : unsigned nparam;
338 : isl_bool empty;
339 : int i;
340 0 : int p = 0, v = 0;
341 :
342 0 : n_div = isl_basic_set_dim(bset, isl_dim_div);
343 0 : d = isl_basic_set_dim(bset, isl_dim_set);
344 0 : nparam = isl_basic_set_dim(bset, isl_dim_param);
345 :
346 0 : for (i = 0; i < n_div; ++i) {
347 0 : if (isl_int_is_zero(c[1 + nparam + d + i]))
348 0 : continue;
349 0 : switch (div_purity[i]) {
350 0 : case PURE_PARAM: p = 1; break;
351 0 : case PURE_VAR: v = 1; break;
352 0 : default: return IMPURE;
353 : }
354 : }
355 0 : if (!p && isl_seq_first_non_zero(c + 1, nparam) == -1)
356 0 : return PURE_VAR;
357 0 : if (!v && isl_seq_first_non_zero(c + 1 + nparam, d) == -1)
358 0 : return PURE_PARAM;
359 :
360 0 : empty = parametric_constant_never_positive(bset, c, div_purity);
361 0 : if (eq && empty >= 0 && !empty) {
362 0 : isl_seq_neg(c, c, 1 + nparam + d + n_div);
363 0 : empty = parametric_constant_never_positive(bset, c, div_purity);
364 : }
365 :
366 0 : return empty < 0 ? -1 : empty ? MIXED : IMPURE;
367 : }
368 :
369 : /* Return an array of integers indicating the type of each div in bset.
370 : * If the div is (recursively) defined in terms of only the parameters,
371 : * then the type is PURE_PARAM.
372 : * If the div is (recursively) defined in terms of only the set variables,
373 : * then the type is PURE_VAR.
374 : * Otherwise, the type is IMPURE.
375 : */
376 0 : static __isl_give int *get_div_purity(__isl_keep isl_basic_set *bset)
377 : {
378 : int i, j;
379 : int *div_purity;
380 : unsigned d;
381 : unsigned n_div;
382 : unsigned nparam;
383 :
384 0 : if (!bset)
385 0 : return NULL;
386 :
387 0 : n_div = isl_basic_set_dim(bset, isl_dim_div);
388 0 : d = isl_basic_set_dim(bset, isl_dim_set);
389 0 : nparam = isl_basic_set_dim(bset, isl_dim_param);
390 :
391 0 : div_purity = isl_alloc_array(bset->ctx, int, n_div);
392 0 : if (n_div && !div_purity)
393 0 : return NULL;
394 :
395 0 : for (i = 0; i < bset->n_div; ++i) {
396 0 : int p = 0, v = 0;
397 0 : if (isl_int_is_zero(bset->div[i][0])) {
398 0 : div_purity[i] = IMPURE;
399 0 : continue;
400 : }
401 0 : if (isl_seq_first_non_zero(bset->div[i] + 2, nparam) != -1)
402 0 : p = 1;
403 0 : if (isl_seq_first_non_zero(bset->div[i] + 2 + nparam, d) != -1)
404 0 : v = 1;
405 0 : for (j = 0; j < i; ++j) {
406 0 : if (isl_int_is_zero(bset->div[i][2 + nparam + d + j]))
407 0 : continue;
408 0 : switch (div_purity[j]) {
409 0 : case PURE_PARAM: p = 1; break;
410 0 : case PURE_VAR: v = 1; break;
411 0 : default: p = v = 1; break;
412 : }
413 : }
414 0 : div_purity[i] = v ? p ? IMPURE : PURE_VAR : PURE_PARAM;
415 : }
416 :
417 0 : return div_purity;
418 : }
419 :
420 : /* Given a path with the as yet unconstrained length at position "pos",
421 : * check if setting the length to zero results in only the identity
422 : * mapping.
423 : */
424 0 : static int empty_path_is_identity(__isl_keep isl_basic_map *path, unsigned pos)
425 : {
426 0 : isl_basic_map *test = NULL;
427 0 : isl_basic_map *id = NULL;
428 : int k;
429 : int is_id;
430 :
431 0 : test = isl_basic_map_copy(path);
432 0 : test = isl_basic_map_extend_constraints(test, 1, 0);
433 0 : k = isl_basic_map_alloc_equality(test);
434 0 : if (k < 0)
435 0 : goto error;
436 0 : isl_seq_clr(test->eq[k], 1 + isl_basic_map_total_dim(test));
437 0 : isl_int_set_si(test->eq[k][pos], 1);
438 0 : test = isl_basic_map_gauss(test, NULL);
439 0 : id = isl_basic_map_identity(isl_basic_map_get_space(path));
440 0 : is_id = isl_basic_map_is_equal(test, id);
441 0 : isl_basic_map_free(test);
442 0 : isl_basic_map_free(id);
443 0 : return is_id;
444 : error:
445 0 : isl_basic_map_free(test);
446 0 : return -1;
447 : }
448 :
449 : /* If any of the constraints is found to be impure then this function
450 : * sets *impurity to 1.
451 : *
452 : * If impurity is NULL then we are dealing with a non-parametric set
453 : * and so the constraints are obviously PURE_VAR.
454 : */
455 0 : static __isl_give isl_basic_map *add_delta_constraints(
456 : __isl_take isl_basic_map *path,
457 : __isl_keep isl_basic_set *delta, unsigned off, unsigned nparam,
458 : unsigned d, int *div_purity, int eq, int *impurity)
459 : {
460 : int i, k;
461 0 : int n = eq ? delta->n_eq : delta->n_ineq;
462 0 : isl_int **delta_c = eq ? delta->eq : delta->ineq;
463 : unsigned n_div;
464 :
465 0 : n_div = isl_basic_set_dim(delta, isl_dim_div);
466 :
467 0 : for (i = 0; i < n; ++i) {
468 : isl_int *path_c;
469 0 : int p = PURE_VAR;
470 0 : if (impurity)
471 0 : p = purity(delta, delta_c[i], div_purity, eq);
472 0 : if (p < 0)
473 0 : goto error;
474 0 : if (p != PURE_VAR && p != PURE_PARAM && !*impurity)
475 0 : *impurity = 1;
476 0 : if (p == IMPURE)
477 0 : continue;
478 0 : if (eq && p != MIXED) {
479 0 : k = isl_basic_map_alloc_equality(path);
480 0 : if (k < 0)
481 0 : goto error;
482 0 : path_c = path->eq[k];
483 : } else {
484 0 : k = isl_basic_map_alloc_inequality(path);
485 0 : if (k < 0)
486 0 : goto error;
487 0 : path_c = path->ineq[k];
488 : }
489 0 : isl_seq_clr(path_c, 1 + isl_basic_map_total_dim(path));
490 0 : if (p == PURE_VAR) {
491 0 : isl_seq_cpy(path_c + off,
492 0 : delta_c[i] + 1 + nparam, d);
493 0 : isl_int_set(path_c[off + d], delta_c[i][0]);
494 0 : } else if (p == PURE_PARAM) {
495 0 : isl_seq_cpy(path_c, delta_c[i], 1 + nparam);
496 : } else {
497 0 : isl_seq_cpy(path_c + off,
498 0 : delta_c[i] + 1 + nparam, d);
499 0 : isl_seq_cpy(path_c, delta_c[i], 1 + nparam);
500 : }
501 0 : isl_seq_cpy(path_c + off - n_div,
502 0 : delta_c[i] + 1 + nparam + d, n_div);
503 : }
504 :
505 0 : return path;
506 : error:
507 0 : isl_basic_map_free(path);
508 0 : return NULL;
509 : }
510 :
511 : /* Given a set of offsets "delta", construct a relation of the
512 : * given dimension specification (Z^{n+1} -> Z^{n+1}) that
513 : * is an overapproximation of the relations that
514 : * maps an element x to any element that can be reached
515 : * by taking a non-negative number of steps along any of
516 : * the elements in "delta".
517 : * That is, construct an approximation of
518 : *
519 : * { [x] -> [y] : exists f \in \delta, k \in Z :
520 : * y = x + k [f, 1] and k >= 0 }
521 : *
522 : * For any element in this relation, the number of steps taken
523 : * is equal to the difference in the final coordinates.
524 : *
525 : * In particular, let delta be defined as
526 : *
527 : * \delta = [p] -> { [x] : A x + a >= 0 and B p + b >= 0 and
528 : * C x + C'p + c >= 0 and
529 : * D x + D'p + d >= 0 }
530 : *
531 : * where the constraints C x + C'p + c >= 0 are such that the parametric
532 : * constant term of each constraint j, "C_j x + C'_j p + c_j",
533 : * can never attain positive values, then the relation is constructed as
534 : *
535 : * { [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and
536 : * A f + k a >= 0 and B p + b >= 0 and
537 : * C f + C'p + c >= 0 and k >= 1 }
538 : * union { [x] -> [x] }
539 : *
540 : * If the zero-length paths happen to correspond exactly to the identity
541 : * mapping, then we return
542 : *
543 : * { [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and
544 : * A f + k a >= 0 and B p + b >= 0 and
545 : * C f + C'p + c >= 0 and k >= 0 }
546 : *
547 : * instead.
548 : *
549 : * Existentially quantified variables in \delta are handled by
550 : * classifying them as independent of the parameters, purely
551 : * parameter dependent and others. Constraints containing
552 : * any of the other existentially quantified variables are removed.
553 : * This is safe, but leads to an additional overapproximation.
554 : *
555 : * If there are any impure constraints, then we also eliminate
556 : * the parameters from \delta, resulting in a set
557 : *
558 : * \delta' = { [x] : E x + e >= 0 }
559 : *
560 : * and add the constraints
561 : *
562 : * E f + k e >= 0
563 : *
564 : * to the constructed relation.
565 : */
566 0 : static __isl_give isl_map *path_along_delta(__isl_take isl_space *dim,
567 : __isl_take isl_basic_set *delta)
568 : {
569 0 : isl_basic_map *path = NULL;
570 : unsigned d;
571 : unsigned n_div;
572 : unsigned nparam;
573 : unsigned off;
574 : int i, k;
575 : int is_id;
576 0 : int *div_purity = NULL;
577 0 : int impurity = 0;
578 :
579 0 : if (!delta)
580 0 : goto error;
581 0 : n_div = isl_basic_set_dim(delta, isl_dim_div);
582 0 : d = isl_basic_set_dim(delta, isl_dim_set);
583 0 : nparam = isl_basic_set_dim(delta, isl_dim_param);
584 0 : path = isl_basic_map_alloc_space(isl_space_copy(dim), n_div + d + 1,
585 0 : d + 1 + delta->n_eq, delta->n_eq + delta->n_ineq + 1);
586 0 : off = 1 + nparam + 2 * (d + 1) + n_div;
587 :
588 0 : for (i = 0; i < n_div + d + 1; ++i) {
589 0 : k = isl_basic_map_alloc_div(path);
590 0 : if (k < 0)
591 0 : goto error;
592 0 : isl_int_set_si(path->div[k][0], 0);
593 : }
594 :
595 0 : for (i = 0; i < d + 1; ++i) {
596 0 : k = isl_basic_map_alloc_equality(path);
597 0 : if (k < 0)
598 0 : goto error;
599 0 : isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
600 0 : isl_int_set_si(path->eq[k][1 + nparam + i], 1);
601 0 : isl_int_set_si(path->eq[k][1 + nparam + d + 1 + i], -1);
602 0 : isl_int_set_si(path->eq[k][off + i], 1);
603 : }
604 :
605 0 : div_purity = get_div_purity(delta);
606 0 : if (n_div && !div_purity)
607 0 : goto error;
608 :
609 0 : path = add_delta_constraints(path, delta, off, nparam, d,
610 : div_purity, 1, &impurity);
611 0 : path = add_delta_constraints(path, delta, off, nparam, d,
612 : div_purity, 0, &impurity);
613 0 : if (impurity) {
614 0 : isl_space *dim = isl_basic_set_get_space(delta);
615 0 : delta = isl_basic_set_project_out(delta,
616 : isl_dim_param, 0, nparam);
617 0 : delta = isl_basic_set_add_dims(delta, isl_dim_param, nparam);
618 0 : delta = isl_basic_set_reset_space(delta, dim);
619 0 : if (!delta)
620 0 : goto error;
621 0 : path = isl_basic_map_extend_constraints(path, delta->n_eq,
622 0 : delta->n_ineq + 1);
623 0 : path = add_delta_constraints(path, delta, off, nparam, d,
624 : NULL, 1, NULL);
625 0 : path = add_delta_constraints(path, delta, off, nparam, d,
626 : NULL, 0, NULL);
627 0 : path = isl_basic_map_gauss(path, NULL);
628 : }
629 :
630 0 : is_id = empty_path_is_identity(path, off + d);
631 0 : if (is_id < 0)
632 0 : goto error;
633 :
634 0 : k = isl_basic_map_alloc_inequality(path);
635 0 : if (k < 0)
636 0 : goto error;
637 0 : isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
638 0 : if (!is_id)
639 0 : isl_int_set_si(path->ineq[k][0], -1);
640 0 : isl_int_set_si(path->ineq[k][off + d], 1);
641 :
642 0 : free(div_purity);
643 0 : isl_basic_set_free(delta);
644 0 : path = isl_basic_map_finalize(path);
645 0 : if (is_id) {
646 0 : isl_space_free(dim);
647 0 : return isl_map_from_basic_map(path);
648 : }
649 0 : return isl_basic_map_union(path, isl_basic_map_identity(dim));
650 : error:
651 0 : free(div_purity);
652 0 : isl_space_free(dim);
653 0 : isl_basic_set_free(delta);
654 0 : isl_basic_map_free(path);
655 0 : return NULL;
656 : }
657 :
658 : /* Given a dimension specification Z^{n+1} -> Z^{n+1} and a parameter "param",
659 : * construct a map that equates the parameter to the difference
660 : * in the final coordinates and imposes that this difference is positive.
661 : * That is, construct
662 : *
663 : * { [x,x_s] -> [y,y_s] : k = y_s - x_s > 0 }
664 : */
665 0 : static __isl_give isl_map *equate_parameter_to_length(__isl_take isl_space *dim,
666 : unsigned param)
667 : {
668 : struct isl_basic_map *bmap;
669 : unsigned d;
670 : unsigned nparam;
671 : int k;
672 :
673 0 : d = isl_space_dim(dim, isl_dim_in);
674 0 : nparam = isl_space_dim(dim, isl_dim_param);
675 0 : bmap = isl_basic_map_alloc_space(dim, 0, 1, 1);
676 0 : k = isl_basic_map_alloc_equality(bmap);
677 0 : if (k < 0)
678 0 : goto error;
679 0 : isl_seq_clr(bmap->eq[k], 1 + isl_basic_map_total_dim(bmap));
680 0 : isl_int_set_si(bmap->eq[k][1 + param], -1);
681 0 : isl_int_set_si(bmap->eq[k][1 + nparam + d - 1], -1);
682 0 : isl_int_set_si(bmap->eq[k][1 + nparam + d + d - 1], 1);
683 :
684 0 : k = isl_basic_map_alloc_inequality(bmap);
685 0 : if (k < 0)
686 0 : goto error;
687 0 : isl_seq_clr(bmap->ineq[k], 1 + isl_basic_map_total_dim(bmap));
688 0 : isl_int_set_si(bmap->ineq[k][1 + param], 1);
689 0 : isl_int_set_si(bmap->ineq[k][0], -1);
690 :
691 0 : bmap = isl_basic_map_finalize(bmap);
692 0 : return isl_map_from_basic_map(bmap);
693 : error:
694 0 : isl_basic_map_free(bmap);
695 0 : return NULL;
696 : }
697 :
698 : /* Check whether "path" is acyclic, where the last coordinates of domain
699 : * and range of path encode the number of steps taken.
700 : * That is, check whether
701 : *
702 : * { d | d = y - x and (x,y) in path }
703 : *
704 : * does not contain any element with positive last coordinate (positive length)
705 : * and zero remaining coordinates (cycle).
706 : */
707 0 : static int is_acyclic(__isl_take isl_map *path)
708 : {
709 : int i;
710 : int acyclic;
711 : unsigned dim;
712 : struct isl_set *delta;
713 :
714 0 : delta = isl_map_deltas(path);
715 0 : dim = isl_set_dim(delta, isl_dim_set);
716 0 : for (i = 0; i < dim; ++i) {
717 0 : if (i == dim -1)
718 0 : delta = isl_set_lower_bound_si(delta, isl_dim_set, i, 1);
719 : else
720 0 : delta = isl_set_fix_si(delta, isl_dim_set, i, 0);
721 : }
722 :
723 0 : acyclic = isl_set_is_empty(delta);
724 0 : isl_set_free(delta);
725 :
726 0 : return acyclic;
727 : }
728 :
729 : /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
730 : * and a dimension specification (Z^{n+1} -> Z^{n+1}),
731 : * construct a map that is an overapproximation of the map
732 : * that takes an element from the space D \times Z to another
733 : * element from the same space, such that the first n coordinates of the
734 : * difference between them is a sum of differences between images
735 : * and pre-images in one of the R_i and such that the last coordinate
736 : * is equal to the number of steps taken.
737 : * That is, let
738 : *
739 : * \Delta_i = { y - x | (x, y) in R_i }
740 : *
741 : * then the constructed map is an overapproximation of
742 : *
743 : * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
744 : * d = (\sum_i k_i \delta_i, \sum_i k_i) }
745 : *
746 : * The elements of the singleton \Delta_i's are collected as the
747 : * rows of the steps matrix. For all these \Delta_i's together,
748 : * a single path is constructed.
749 : * For each of the other \Delta_i's, we compute an overapproximation
750 : * of the paths along elements of \Delta_i.
751 : * Since each of these paths performs an addition, composition is
752 : * symmetric and we can simply compose all resulting paths in any order.
753 : */
754 0 : static __isl_give isl_map *construct_extended_path(__isl_take isl_space *dim,
755 : __isl_keep isl_map *map, int *project)
756 : {
757 0 : struct isl_mat *steps = NULL;
758 0 : struct isl_map *path = NULL;
759 : unsigned d;
760 : int i, j, n;
761 :
762 0 : if (!map)
763 0 : goto error;
764 :
765 0 : d = isl_map_dim(map, isl_dim_in);
766 :
767 0 : path = isl_map_identity(isl_space_copy(dim));
768 :
769 0 : steps = isl_mat_alloc(map->ctx, map->n, d);
770 0 : if (!steps)
771 0 : goto error;
772 :
773 0 : n = 0;
774 0 : for (i = 0; i < map->n; ++i) {
775 : struct isl_basic_set *delta;
776 :
777 0 : delta = isl_basic_map_deltas(isl_basic_map_copy(map->p[i]));
778 :
779 0 : for (j = 0; j < d; ++j) {
780 : isl_bool fixed;
781 :
782 0 : fixed = isl_basic_set_plain_dim_is_fixed(delta, j,
783 0 : &steps->row[n][j]);
784 0 : if (fixed < 0) {
785 0 : isl_basic_set_free(delta);
786 0 : goto error;
787 : }
788 0 : if (!fixed)
789 0 : break;
790 : }
791 :
792 :
793 0 : if (j < d) {
794 0 : path = isl_map_apply_range(path,
795 : path_along_delta(isl_space_copy(dim), delta));
796 0 : path = isl_map_coalesce(path);
797 : } else {
798 0 : isl_basic_set_free(delta);
799 0 : ++n;
800 : }
801 : }
802 :
803 0 : if (n > 0) {
804 0 : steps->n_row = n;
805 0 : path = isl_map_apply_range(path,
806 : path_along_steps(isl_space_copy(dim), steps));
807 : }
808 :
809 0 : if (project && *project) {
810 0 : *project = is_acyclic(isl_map_copy(path));
811 0 : if (*project < 0)
812 0 : goto error;
813 : }
814 :
815 0 : isl_space_free(dim);
816 0 : isl_mat_free(steps);
817 0 : return path;
818 : error:
819 0 : isl_space_free(dim);
820 0 : isl_mat_free(steps);
821 0 : isl_map_free(path);
822 0 : return NULL;
823 : }
824 :
825 0 : static isl_bool isl_set_overlaps(__isl_keep isl_set *set1,
826 : __isl_keep isl_set *set2)
827 : {
828 : isl_set *i;
829 : isl_bool no_overlap;
830 :
831 0 : if (!set1 || !set2)
832 0 : return isl_bool_error;
833 :
834 0 : if (!isl_space_tuple_is_equal(set1->dim, isl_dim_set,
835 : set2->dim, isl_dim_set))
836 0 : return isl_bool_false;
837 :
838 0 : i = isl_set_intersect(isl_set_copy(set1), isl_set_copy(set2));
839 0 : no_overlap = isl_set_is_empty(i);
840 0 : isl_set_free(i);
841 :
842 0 : return isl_bool_not(no_overlap);
843 : }
844 :
845 : /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
846 : * and a dimension specification (Z^{n+1} -> Z^{n+1}),
847 : * construct a map that is an overapproximation of the map
848 : * that takes an element from the dom R \times Z to an
849 : * element from ran R \times Z, such that the first n coordinates of the
850 : * difference between them is a sum of differences between images
851 : * and pre-images in one of the R_i and such that the last coordinate
852 : * is equal to the number of steps taken.
853 : * That is, let
854 : *
855 : * \Delta_i = { y - x | (x, y) in R_i }
856 : *
857 : * then the constructed map is an overapproximation of
858 : *
859 : * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
860 : * d = (\sum_i k_i \delta_i, \sum_i k_i) and
861 : * x in dom R and x + d in ran R and
862 : * \sum_i k_i >= 1 }
863 : */
864 0 : static __isl_give isl_map *construct_component(__isl_take isl_space *dim,
865 : __isl_keep isl_map *map, int *exact, int project)
866 : {
867 0 : struct isl_set *domain = NULL;
868 0 : struct isl_set *range = NULL;
869 0 : struct isl_map *app = NULL;
870 0 : struct isl_map *path = NULL;
871 : isl_bool overlaps;
872 :
873 0 : domain = isl_map_domain(isl_map_copy(map));
874 0 : domain = isl_set_coalesce(domain);
875 0 : range = isl_map_range(isl_map_copy(map));
876 0 : range = isl_set_coalesce(range);
877 0 : overlaps = isl_set_overlaps(domain, range);
878 0 : if (overlaps < 0 || !overlaps) {
879 0 : isl_set_free(domain);
880 0 : isl_set_free(range);
881 0 : isl_space_free(dim);
882 :
883 0 : if (overlaps < 0)
884 0 : map = NULL;
885 0 : map = isl_map_copy(map);
886 0 : map = isl_map_add_dims(map, isl_dim_in, 1);
887 0 : map = isl_map_add_dims(map, isl_dim_out, 1);
888 0 : map = set_path_length(map, 1, 1);
889 0 : return map;
890 : }
891 0 : app = isl_map_from_domain_and_range(domain, range);
892 0 : app = isl_map_add_dims(app, isl_dim_in, 1);
893 0 : app = isl_map_add_dims(app, isl_dim_out, 1);
894 :
895 0 : path = construct_extended_path(isl_space_copy(dim), map,
896 0 : exact && *exact ? &project : NULL);
897 0 : app = isl_map_intersect(app, path);
898 :
899 0 : if (exact && *exact &&
900 0 : (*exact = check_exactness(isl_map_copy(map), isl_map_copy(app),
901 : project)) < 0)
902 0 : goto error;
903 :
904 0 : isl_space_free(dim);
905 0 : app = set_path_length(app, 0, 1);
906 0 : return app;
907 : error:
908 0 : isl_space_free(dim);
909 0 : isl_map_free(app);
910 0 : return NULL;
911 : }
912 :
913 : /* Call construct_component and, if "project" is set, project out
914 : * the final coordinates.
915 : */
916 0 : static __isl_give isl_map *construct_projected_component(
917 : __isl_take isl_space *dim,
918 : __isl_keep isl_map *map, int *exact, int project)
919 : {
920 : isl_map *app;
921 : unsigned d;
922 :
923 0 : if (!dim)
924 0 : return NULL;
925 0 : d = isl_space_dim(dim, isl_dim_in);
926 :
927 0 : app = construct_component(dim, map, exact, project);
928 0 : if (project) {
929 0 : app = isl_map_project_out(app, isl_dim_in, d - 1, 1);
930 0 : app = isl_map_project_out(app, isl_dim_out, d - 1, 1);
931 : }
932 0 : return app;
933 : }
934 :
935 : /* Compute an extended version, i.e., with path lengths, of
936 : * an overapproximation of the transitive closure of "bmap"
937 : * with path lengths greater than or equal to zero and with
938 : * domain and range equal to "dom".
939 : */
940 0 : static __isl_give isl_map *q_closure(__isl_take isl_space *dim,
941 : __isl_take isl_set *dom, __isl_keep isl_basic_map *bmap, int *exact)
942 : {
943 0 : int project = 1;
944 : isl_map *path;
945 : isl_map *map;
946 : isl_map *app;
947 :
948 0 : dom = isl_set_add_dims(dom, isl_dim_set, 1);
949 0 : app = isl_map_from_domain_and_range(dom, isl_set_copy(dom));
950 0 : map = isl_map_from_basic_map(isl_basic_map_copy(bmap));
951 0 : path = construct_extended_path(dim, map, &project);
952 0 : app = isl_map_intersect(app, path);
953 :
954 0 : if ((*exact = check_exactness(map, isl_map_copy(app), project)) < 0)
955 0 : goto error;
956 :
957 0 : return app;
958 : error:
959 0 : isl_map_free(app);
960 0 : return NULL;
961 : }
962 :
963 : /* Check whether qc has any elements of length at least one
964 : * with domain and/or range outside of dom and ran.
965 : */
966 0 : static int has_spurious_elements(__isl_keep isl_map *qc,
967 : __isl_keep isl_set *dom, __isl_keep isl_set *ran)
968 : {
969 : isl_set *s;
970 : int subset;
971 : unsigned d;
972 :
973 0 : if (!qc || !dom || !ran)
974 0 : return -1;
975 :
976 0 : d = isl_map_dim(qc, isl_dim_in);
977 :
978 0 : qc = isl_map_copy(qc);
979 0 : qc = set_path_length(qc, 0, 1);
980 0 : qc = isl_map_project_out(qc, isl_dim_in, d - 1, 1);
981 0 : qc = isl_map_project_out(qc, isl_dim_out, d - 1, 1);
982 :
983 0 : s = isl_map_domain(isl_map_copy(qc));
984 0 : subset = isl_set_is_subset(s, dom);
985 0 : isl_set_free(s);
986 0 : if (subset < 0)
987 0 : goto error;
988 0 : if (!subset) {
989 0 : isl_map_free(qc);
990 0 : return 1;
991 : }
992 :
993 0 : s = isl_map_range(qc);
994 0 : subset = isl_set_is_subset(s, ran);
995 0 : isl_set_free(s);
996 :
997 0 : return subset < 0 ? -1 : !subset;
998 : error:
999 0 : isl_map_free(qc);
1000 0 : return -1;
1001 : }
1002 :
1003 : #define LEFT 2
1004 : #define RIGHT 1
1005 :
1006 : /* For each basic map in "map", except i, check whether it combines
1007 : * with the transitive closure that is reflexive on C combines
1008 : * to the left and to the right.
1009 : *
1010 : * In particular, if
1011 : *
1012 : * dom map_j \subseteq C
1013 : *
1014 : * then right[j] is set to 1. Otherwise, if
1015 : *
1016 : * ran map_i \cap dom map_j = \emptyset
1017 : *
1018 : * then right[j] is set to 0. Otherwise, composing to the right
1019 : * is impossible.
1020 : *
1021 : * Similar, for composing to the left, we have if
1022 : *
1023 : * ran map_j \subseteq C
1024 : *
1025 : * then left[j] is set to 1. Otherwise, if
1026 : *
1027 : * dom map_i \cap ran map_j = \emptyset
1028 : *
1029 : * then left[j] is set to 0. Otherwise, composing to the left
1030 : * is impossible.
1031 : *
1032 : * The return value is or'd with LEFT if composing to the left
1033 : * is possible and with RIGHT if composing to the right is possible.
1034 : */
1035 0 : static int composability(__isl_keep isl_set *C, int i,
1036 : isl_set **dom, isl_set **ran, int *left, int *right,
1037 : __isl_keep isl_map *map)
1038 : {
1039 : int j;
1040 : int ok;
1041 :
1042 0 : ok = LEFT | RIGHT;
1043 0 : for (j = 0; j < map->n && ok; ++j) {
1044 : isl_bool overlaps, subset;
1045 0 : if (j == i)
1046 0 : continue;
1047 :
1048 0 : if (ok & RIGHT) {
1049 0 : if (!dom[j])
1050 0 : dom[j] = isl_set_from_basic_set(
1051 : isl_basic_map_domain(
1052 0 : isl_basic_map_copy(map->p[j])));
1053 0 : if (!dom[j])
1054 0 : return -1;
1055 0 : overlaps = isl_set_overlaps(ran[i], dom[j]);
1056 0 : if (overlaps < 0)
1057 0 : return -1;
1058 0 : if (!overlaps)
1059 0 : right[j] = 0;
1060 : else {
1061 0 : subset = isl_set_is_subset(dom[j], C);
1062 0 : if (subset < 0)
1063 0 : return -1;
1064 0 : if (subset)
1065 0 : right[j] = 1;
1066 : else
1067 0 : ok &= ~RIGHT;
1068 : }
1069 : }
1070 :
1071 0 : if (ok & LEFT) {
1072 0 : if (!ran[j])
1073 0 : ran[j] = isl_set_from_basic_set(
1074 : isl_basic_map_range(
1075 0 : isl_basic_map_copy(map->p[j])));
1076 0 : if (!ran[j])
1077 0 : return -1;
1078 0 : overlaps = isl_set_overlaps(dom[i], ran[j]);
1079 0 : if (overlaps < 0)
1080 0 : return -1;
1081 0 : if (!overlaps)
1082 0 : left[j] = 0;
1083 : else {
1084 0 : subset = isl_set_is_subset(ran[j], C);
1085 0 : if (subset < 0)
1086 0 : return -1;
1087 0 : if (subset)
1088 0 : left[j] = 1;
1089 : else
1090 0 : ok &= ~LEFT;
1091 : }
1092 : }
1093 : }
1094 :
1095 0 : return ok;
1096 : }
1097 :
1098 0 : static __isl_give isl_map *anonymize(__isl_take isl_map *map)
1099 : {
1100 0 : map = isl_map_reset(map, isl_dim_in);
1101 0 : map = isl_map_reset(map, isl_dim_out);
1102 0 : return map;
1103 : }
1104 :
1105 : /* Return a map that is a union of the basic maps in "map", except i,
1106 : * composed to left and right with qc based on the entries of "left"
1107 : * and "right".
1108 : */
1109 0 : static __isl_give isl_map *compose(__isl_keep isl_map *map, int i,
1110 : __isl_take isl_map *qc, int *left, int *right)
1111 : {
1112 : int j;
1113 : isl_map *comp;
1114 :
1115 0 : comp = isl_map_empty(isl_map_get_space(map));
1116 0 : for (j = 0; j < map->n; ++j) {
1117 : isl_map *map_j;
1118 :
1119 0 : if (j == i)
1120 0 : continue;
1121 :
1122 0 : map_j = isl_map_from_basic_map(isl_basic_map_copy(map->p[j]));
1123 0 : map_j = anonymize(map_j);
1124 0 : if (left && left[j])
1125 0 : map_j = isl_map_apply_range(map_j, isl_map_copy(qc));
1126 0 : if (right && right[j])
1127 0 : map_j = isl_map_apply_range(isl_map_copy(qc), map_j);
1128 0 : comp = isl_map_union(comp, map_j);
1129 : }
1130 :
1131 0 : comp = isl_map_compute_divs(comp);
1132 0 : comp = isl_map_coalesce(comp);
1133 :
1134 0 : isl_map_free(qc);
1135 :
1136 0 : return comp;
1137 : }
1138 :
1139 : /* Compute the transitive closure of "map" incrementally by
1140 : * computing
1141 : *
1142 : * map_i^+ \cup qc^+
1143 : *
1144 : * or
1145 : *
1146 : * map_i^+ \cup ((id \cup map_i^) \circ qc^+)
1147 : *
1148 : * or
1149 : *
1150 : * map_i^+ \cup (qc^+ \circ (id \cup map_i^))
1151 : *
1152 : * depending on whether left or right are NULL.
1153 : */
1154 0 : static __isl_give isl_map *compute_incremental(
1155 : __isl_take isl_space *dim, __isl_keep isl_map *map,
1156 : int i, __isl_take isl_map *qc, int *left, int *right, int *exact)
1157 : {
1158 : isl_map *map_i;
1159 : isl_map *tc;
1160 0 : isl_map *rtc = NULL;
1161 :
1162 0 : if (!map)
1163 0 : goto error;
1164 0 : isl_assert(map->ctx, left || right, goto error);
1165 :
1166 0 : map_i = isl_map_from_basic_map(isl_basic_map_copy(map->p[i]));
1167 0 : tc = construct_projected_component(isl_space_copy(dim), map_i,
1168 : exact, 1);
1169 0 : isl_map_free(map_i);
1170 :
1171 0 : if (*exact)
1172 0 : qc = isl_map_transitive_closure(qc, exact);
1173 :
1174 0 : if (!*exact) {
1175 0 : isl_space_free(dim);
1176 0 : isl_map_free(tc);
1177 0 : isl_map_free(qc);
1178 0 : return isl_map_universe(isl_map_get_space(map));
1179 : }
1180 :
1181 0 : if (!left || !right)
1182 0 : rtc = isl_map_union(isl_map_copy(tc),
1183 : isl_map_identity(isl_map_get_space(tc)));
1184 0 : if (!right)
1185 0 : qc = isl_map_apply_range(rtc, qc);
1186 0 : if (!left)
1187 0 : qc = isl_map_apply_range(qc, rtc);
1188 0 : qc = isl_map_union(tc, qc);
1189 :
1190 0 : isl_space_free(dim);
1191 :
1192 0 : return qc;
1193 : error:
1194 0 : isl_space_free(dim);
1195 0 : isl_map_free(qc);
1196 0 : return NULL;
1197 : }
1198 :
1199 : /* Given a map "map", try to find a basic map such that
1200 : * map^+ can be computed as
1201 : *
1202 : * map^+ = map_i^+ \cup
1203 : * \bigcup_j ((map_i^+ \cup Id_C)^+ \circ map_j \circ (map_i^+ \cup Id_C))^+
1204 : *
1205 : * with C the simple hull of the domain and range of the input map.
1206 : * map_i^ \cup Id_C is computed by allowing the path lengths to be zero
1207 : * and by intersecting domain and range with C.
1208 : * Of course, we need to check that this is actually equal to map_i^ \cup Id_C.
1209 : * Also, we only use the incremental computation if all the transitive
1210 : * closures are exact and if the number of basic maps in the union,
1211 : * after computing the integer divisions, is smaller than the number
1212 : * of basic maps in the input map.
1213 : */
1214 0 : static int incemental_on_entire_domain(__isl_keep isl_space *dim,
1215 : __isl_keep isl_map *map,
1216 : isl_set **dom, isl_set **ran, int *left, int *right,
1217 : __isl_give isl_map **res)
1218 : {
1219 : int i;
1220 : isl_set *C;
1221 : unsigned d;
1222 :
1223 0 : *res = NULL;
1224 :
1225 0 : C = isl_set_union(isl_map_domain(isl_map_copy(map)),
1226 : isl_map_range(isl_map_copy(map)));
1227 0 : C = isl_set_from_basic_set(isl_set_simple_hull(C));
1228 0 : if (!C)
1229 0 : return -1;
1230 0 : if (C->n != 1) {
1231 0 : isl_set_free(C);
1232 0 : return 0;
1233 : }
1234 :
1235 0 : d = isl_map_dim(map, isl_dim_in);
1236 :
1237 0 : for (i = 0; i < map->n; ++i) {
1238 : isl_map *qc;
1239 : int exact_i, spurious;
1240 : int j;
1241 0 : dom[i] = isl_set_from_basic_set(isl_basic_map_domain(
1242 0 : isl_basic_map_copy(map->p[i])));
1243 0 : ran[i] = isl_set_from_basic_set(isl_basic_map_range(
1244 0 : isl_basic_map_copy(map->p[i])));
1245 0 : qc = q_closure(isl_space_copy(dim), isl_set_copy(C),
1246 0 : map->p[i], &exact_i);
1247 0 : if (!qc)
1248 0 : goto error;
1249 0 : if (!exact_i) {
1250 0 : isl_map_free(qc);
1251 0 : continue;
1252 : }
1253 0 : spurious = has_spurious_elements(qc, dom[i], ran[i]);
1254 0 : if (spurious) {
1255 0 : isl_map_free(qc);
1256 0 : if (spurious < 0)
1257 0 : goto error;
1258 0 : continue;
1259 : }
1260 0 : qc = isl_map_project_out(qc, isl_dim_in, d, 1);
1261 0 : qc = isl_map_project_out(qc, isl_dim_out, d, 1);
1262 0 : qc = isl_map_compute_divs(qc);
1263 0 : for (j = 0; j < map->n; ++j)
1264 0 : left[j] = right[j] = 1;
1265 0 : qc = compose(map, i, qc, left, right);
1266 0 : if (!qc)
1267 0 : goto error;
1268 0 : if (qc->n >= map->n) {
1269 0 : isl_map_free(qc);
1270 0 : continue;
1271 : }
1272 0 : *res = compute_incremental(isl_space_copy(dim), map, i, qc,
1273 : left, right, &exact_i);
1274 0 : if (!*res)
1275 0 : goto error;
1276 0 : if (exact_i)
1277 0 : break;
1278 0 : isl_map_free(*res);
1279 0 : *res = NULL;
1280 : }
1281 :
1282 0 : isl_set_free(C);
1283 :
1284 0 : return *res != NULL;
1285 : error:
1286 0 : isl_set_free(C);
1287 0 : return -1;
1288 : }
1289 :
1290 : /* Try and compute the transitive closure of "map" as
1291 : *
1292 : * map^+ = map_i^+ \cup
1293 : * \bigcup_j ((map_i^+ \cup Id_C)^+ \circ map_j \circ (map_i^+ \cup Id_C))^+
1294 : *
1295 : * with C either the simple hull of the domain and range of the entire
1296 : * map or the simple hull of domain and range of map_i.
1297 : */
1298 0 : static __isl_give isl_map *incremental_closure(__isl_take isl_space *dim,
1299 : __isl_keep isl_map *map, int *exact, int project)
1300 : {
1301 : int i;
1302 0 : isl_set **dom = NULL;
1303 0 : isl_set **ran = NULL;
1304 0 : int *left = NULL;
1305 0 : int *right = NULL;
1306 : isl_set *C;
1307 : unsigned d;
1308 0 : isl_map *res = NULL;
1309 :
1310 0 : if (!project)
1311 0 : return construct_projected_component(dim, map, exact, project);
1312 :
1313 0 : if (!map)
1314 0 : goto error;
1315 0 : if (map->n <= 1)
1316 0 : return construct_projected_component(dim, map, exact, project);
1317 :
1318 0 : d = isl_map_dim(map, isl_dim_in);
1319 :
1320 0 : dom = isl_calloc_array(map->ctx, isl_set *, map->n);
1321 0 : ran = isl_calloc_array(map->ctx, isl_set *, map->n);
1322 0 : left = isl_calloc_array(map->ctx, int, map->n);
1323 0 : right = isl_calloc_array(map->ctx, int, map->n);
1324 0 : if (!ran || !dom || !left || !right)
1325 : goto error;
1326 :
1327 0 : if (incemental_on_entire_domain(dim, map, dom, ran, left, right, &res) < 0)
1328 0 : goto error;
1329 :
1330 0 : for (i = 0; !res && i < map->n; ++i) {
1331 : isl_map *qc;
1332 : int exact_i, spurious, comp;
1333 0 : if (!dom[i])
1334 0 : dom[i] = isl_set_from_basic_set(
1335 : isl_basic_map_domain(
1336 0 : isl_basic_map_copy(map->p[i])));
1337 0 : if (!dom[i])
1338 0 : goto error;
1339 0 : if (!ran[i])
1340 0 : ran[i] = isl_set_from_basic_set(
1341 : isl_basic_map_range(
1342 0 : isl_basic_map_copy(map->p[i])));
1343 0 : if (!ran[i])
1344 0 : goto error;
1345 0 : C = isl_set_union(isl_set_copy(dom[i]),
1346 0 : isl_set_copy(ran[i]));
1347 0 : C = isl_set_from_basic_set(isl_set_simple_hull(C));
1348 0 : if (!C)
1349 0 : goto error;
1350 0 : if (C->n != 1) {
1351 0 : isl_set_free(C);
1352 0 : continue;
1353 : }
1354 0 : comp = composability(C, i, dom, ran, left, right, map);
1355 0 : if (!comp || comp < 0) {
1356 0 : isl_set_free(C);
1357 0 : if (comp < 0)
1358 0 : goto error;
1359 0 : continue;
1360 : }
1361 0 : qc = q_closure(isl_space_copy(dim), C, map->p[i], &exact_i);
1362 0 : if (!qc)
1363 0 : goto error;
1364 0 : if (!exact_i) {
1365 0 : isl_map_free(qc);
1366 0 : continue;
1367 : }
1368 0 : spurious = has_spurious_elements(qc, dom[i], ran[i]);
1369 0 : if (spurious) {
1370 0 : isl_map_free(qc);
1371 0 : if (spurious < 0)
1372 0 : goto error;
1373 0 : continue;
1374 : }
1375 0 : qc = isl_map_project_out(qc, isl_dim_in, d, 1);
1376 0 : qc = isl_map_project_out(qc, isl_dim_out, d, 1);
1377 0 : qc = isl_map_compute_divs(qc);
1378 0 : qc = compose(map, i, qc, (comp & LEFT) ? left : NULL,
1379 0 : (comp & RIGHT) ? right : NULL);
1380 0 : if (!qc)
1381 0 : goto error;
1382 0 : if (qc->n >= map->n) {
1383 0 : isl_map_free(qc);
1384 0 : continue;
1385 : }
1386 0 : res = compute_incremental(isl_space_copy(dim), map, i, qc,
1387 0 : (comp & LEFT) ? left : NULL,
1388 0 : (comp & RIGHT) ? right : NULL, &exact_i);
1389 0 : if (!res)
1390 0 : goto error;
1391 0 : if (exact_i)
1392 0 : break;
1393 0 : isl_map_free(res);
1394 0 : res = NULL;
1395 : }
1396 :
1397 0 : for (i = 0; i < map->n; ++i) {
1398 0 : isl_set_free(dom[i]);
1399 0 : isl_set_free(ran[i]);
1400 : }
1401 0 : free(dom);
1402 0 : free(ran);
1403 0 : free(left);
1404 0 : free(right);
1405 :
1406 0 : if (res) {
1407 0 : isl_space_free(dim);
1408 0 : return res;
1409 : }
1410 :
1411 0 : return construct_projected_component(dim, map, exact, project);
1412 : error:
1413 0 : if (dom)
1414 0 : for (i = 0; i < map->n; ++i)
1415 0 : isl_set_free(dom[i]);
1416 0 : free(dom);
1417 0 : if (ran)
1418 0 : for (i = 0; i < map->n; ++i)
1419 0 : isl_set_free(ran[i]);
1420 0 : free(ran);
1421 0 : free(left);
1422 0 : free(right);
1423 0 : isl_space_free(dim);
1424 0 : return NULL;
1425 : }
1426 :
1427 : /* Given an array of sets "set", add "dom" at position "pos"
1428 : * and search for elements at earlier positions that overlap with "dom".
1429 : * If any can be found, then merge all of them, together with "dom", into
1430 : * a single set and assign the union to the first in the array,
1431 : * which becomes the new group leader for all groups involved in the merge.
1432 : * During the search, we only consider group leaders, i.e., those with
1433 : * group[i] = i, as the other sets have already been combined
1434 : * with one of the group leaders.
1435 : */
1436 0 : static int merge(isl_set **set, int *group, __isl_take isl_set *dom, int pos)
1437 : {
1438 : int i;
1439 :
1440 0 : group[pos] = pos;
1441 0 : set[pos] = isl_set_copy(dom);
1442 :
1443 0 : for (i = pos - 1; i >= 0; --i) {
1444 : isl_bool o;
1445 :
1446 0 : if (group[i] != i)
1447 0 : continue;
1448 :
1449 0 : o = isl_set_overlaps(set[i], dom);
1450 0 : if (o < 0)
1451 0 : goto error;
1452 0 : if (!o)
1453 0 : continue;
1454 :
1455 0 : set[i] = isl_set_union(set[i], set[group[pos]]);
1456 0 : set[group[pos]] = NULL;
1457 0 : if (!set[i])
1458 0 : goto error;
1459 0 : group[group[pos]] = i;
1460 0 : group[pos] = i;
1461 : }
1462 :
1463 0 : isl_set_free(dom);
1464 0 : return 0;
1465 : error:
1466 0 : isl_set_free(dom);
1467 0 : return -1;
1468 : }
1469 :
1470 : /* Replace each entry in the n by n grid of maps by the cross product
1471 : * with the relation { [i] -> [i + 1] }.
1472 : */
1473 0 : static int add_length(__isl_keep isl_map *map, isl_map ***grid, int n)
1474 : {
1475 : int i, j, k;
1476 : isl_space *dim;
1477 : isl_basic_map *bstep;
1478 : isl_map *step;
1479 : unsigned nparam;
1480 :
1481 0 : if (!map)
1482 0 : return -1;
1483 :
1484 0 : dim = isl_map_get_space(map);
1485 0 : nparam = isl_space_dim(dim, isl_dim_param);
1486 0 : dim = isl_space_drop_dims(dim, isl_dim_in, 0, isl_space_dim(dim, isl_dim_in));
1487 0 : dim = isl_space_drop_dims(dim, isl_dim_out, 0, isl_space_dim(dim, isl_dim_out));
1488 0 : dim = isl_space_add_dims(dim, isl_dim_in, 1);
1489 0 : dim = isl_space_add_dims(dim, isl_dim_out, 1);
1490 0 : bstep = isl_basic_map_alloc_space(dim, 0, 1, 0);
1491 0 : k = isl_basic_map_alloc_equality(bstep);
1492 0 : if (k < 0) {
1493 0 : isl_basic_map_free(bstep);
1494 0 : return -1;
1495 : }
1496 0 : isl_seq_clr(bstep->eq[k], 1 + isl_basic_map_total_dim(bstep));
1497 0 : isl_int_set_si(bstep->eq[k][0], 1);
1498 0 : isl_int_set_si(bstep->eq[k][1 + nparam], 1);
1499 0 : isl_int_set_si(bstep->eq[k][1 + nparam + 1], -1);
1500 0 : bstep = isl_basic_map_finalize(bstep);
1501 0 : step = isl_map_from_basic_map(bstep);
1502 :
1503 0 : for (i = 0; i < n; ++i)
1504 0 : for (j = 0; j < n; ++j)
1505 0 : grid[i][j] = isl_map_product(grid[i][j],
1506 : isl_map_copy(step));
1507 :
1508 0 : isl_map_free(step);
1509 :
1510 0 : return 0;
1511 : }
1512 :
1513 : /* The core of the Floyd-Warshall algorithm.
1514 : * Updates the given n x x matrix of relations in place.
1515 : *
1516 : * The algorithm iterates over all vertices. In each step, the whole
1517 : * matrix is updated to include all paths that go to the current vertex,
1518 : * possibly stay there a while (including passing through earlier vertices)
1519 : * and then come back. At the start of each iteration, the diagonal
1520 : * element corresponding to the current vertex is replaced by its
1521 : * transitive closure to account for all indirect paths that stay
1522 : * in the current vertex.
1523 : */
1524 0 : static void floyd_warshall_iterate(isl_map ***grid, int n, int *exact)
1525 : {
1526 : int r, p, q;
1527 :
1528 0 : for (r = 0; r < n; ++r) {
1529 : int r_exact;
1530 0 : grid[r][r] = isl_map_transitive_closure(grid[r][r],
1531 0 : (exact && *exact) ? &r_exact : NULL);
1532 0 : if (exact && *exact && !r_exact)
1533 0 : *exact = 0;
1534 :
1535 0 : for (p = 0; p < n; ++p)
1536 0 : for (q = 0; q < n; ++q) {
1537 : isl_map *loop;
1538 0 : if (p == r && q == r)
1539 0 : continue;
1540 0 : loop = isl_map_apply_range(
1541 0 : isl_map_copy(grid[p][r]),
1542 0 : isl_map_copy(grid[r][q]));
1543 0 : grid[p][q] = isl_map_union(grid[p][q], loop);
1544 0 : loop = isl_map_apply_range(
1545 0 : isl_map_copy(grid[p][r]),
1546 : isl_map_apply_range(
1547 0 : isl_map_copy(grid[r][r]),
1548 0 : isl_map_copy(grid[r][q])));
1549 0 : grid[p][q] = isl_map_union(grid[p][q], loop);
1550 0 : grid[p][q] = isl_map_coalesce(grid[p][q]);
1551 : }
1552 : }
1553 0 : }
1554 :
1555 : /* Given a partition of the domains and ranges of the basic maps in "map",
1556 : * apply the Floyd-Warshall algorithm with the elements in the partition
1557 : * as vertices.
1558 : *
1559 : * In particular, there are "n" elements in the partition and "group" is
1560 : * an array of length 2 * map->n with entries in [0,n-1].
1561 : *
1562 : * We first construct a matrix of relations based on the partition information,
1563 : * apply Floyd-Warshall on this matrix of relations and then take the
1564 : * union of all entries in the matrix as the final result.
1565 : *
1566 : * If we are actually computing the power instead of the transitive closure,
1567 : * i.e., when "project" is not set, then the result should have the
1568 : * path lengths encoded as the difference between an extra pair of
1569 : * coordinates. We therefore apply the nested transitive closures
1570 : * to relations that include these lengths. In particular, we replace
1571 : * the input relation by the cross product with the unit length relation
1572 : * { [i] -> [i + 1] }.
1573 : */
1574 0 : static __isl_give isl_map *floyd_warshall_with_groups(__isl_take isl_space *dim,
1575 : __isl_keep isl_map *map, int *exact, int project, int *group, int n)
1576 : {
1577 : int i, j, k;
1578 0 : isl_map ***grid = NULL;
1579 : isl_map *app;
1580 :
1581 0 : if (!map)
1582 0 : goto error;
1583 :
1584 0 : if (n == 1) {
1585 0 : free(group);
1586 0 : return incremental_closure(dim, map, exact, project);
1587 : }
1588 :
1589 0 : grid = isl_calloc_array(map->ctx, isl_map **, n);
1590 0 : if (!grid)
1591 0 : goto error;
1592 0 : for (i = 0; i < n; ++i) {
1593 0 : grid[i] = isl_calloc_array(map->ctx, isl_map *, n);
1594 0 : if (!grid[i])
1595 0 : goto error;
1596 0 : for (j = 0; j < n; ++j)
1597 0 : grid[i][j] = isl_map_empty(isl_map_get_space(map));
1598 : }
1599 :
1600 0 : for (k = 0; k < map->n; ++k) {
1601 0 : i = group[2 * k];
1602 0 : j = group[2 * k + 1];
1603 0 : grid[i][j] = isl_map_union(grid[i][j],
1604 : isl_map_from_basic_map(
1605 0 : isl_basic_map_copy(map->p[k])));
1606 : }
1607 :
1608 0 : if (!project && add_length(map, grid, n) < 0)
1609 0 : goto error;
1610 :
1611 0 : floyd_warshall_iterate(grid, n, exact);
1612 :
1613 0 : app = isl_map_empty(isl_map_get_space(grid[0][0]));
1614 :
1615 0 : for (i = 0; i < n; ++i) {
1616 0 : for (j = 0; j < n; ++j)
1617 0 : app = isl_map_union(app, grid[i][j]);
1618 0 : free(grid[i]);
1619 : }
1620 0 : free(grid);
1621 :
1622 0 : free(group);
1623 0 : isl_space_free(dim);
1624 :
1625 0 : return app;
1626 : error:
1627 0 : if (grid)
1628 0 : for (i = 0; i < n; ++i) {
1629 0 : if (!grid[i])
1630 0 : continue;
1631 0 : for (j = 0; j < n; ++j)
1632 0 : isl_map_free(grid[i][j]);
1633 0 : free(grid[i]);
1634 : }
1635 0 : free(grid);
1636 0 : free(group);
1637 0 : isl_space_free(dim);
1638 0 : return NULL;
1639 : }
1640 :
1641 : /* Partition the domains and ranges of the n basic relations in list
1642 : * into disjoint cells.
1643 : *
1644 : * To find the partition, we simply consider all of the domains
1645 : * and ranges in turn and combine those that overlap.
1646 : * "set" contains the partition elements and "group" indicates
1647 : * to which partition element a given domain or range belongs.
1648 : * The domain of basic map i corresponds to element 2 * i in these arrays,
1649 : * while the domain corresponds to element 2 * i + 1.
1650 : * During the construction group[k] is either equal to k,
1651 : * in which case set[k] contains the union of all the domains and
1652 : * ranges in the corresponding group, or is equal to some l < k,
1653 : * with l another domain or range in the same group.
1654 : */
1655 0 : static int *setup_groups(isl_ctx *ctx, __isl_keep isl_basic_map **list, int n,
1656 : isl_set ***set, int *n_group)
1657 : {
1658 : int i;
1659 0 : int *group = NULL;
1660 : int g;
1661 :
1662 0 : *set = isl_calloc_array(ctx, isl_set *, 2 * n);
1663 0 : group = isl_alloc_array(ctx, int, 2 * n);
1664 :
1665 0 : if (!*set || !group)
1666 : goto error;
1667 :
1668 0 : for (i = 0; i < n; ++i) {
1669 : isl_set *dom;
1670 0 : dom = isl_set_from_basic_set(isl_basic_map_domain(
1671 0 : isl_basic_map_copy(list[i])));
1672 0 : if (merge(*set, group, dom, 2 * i) < 0)
1673 0 : goto error;
1674 0 : dom = isl_set_from_basic_set(isl_basic_map_range(
1675 0 : isl_basic_map_copy(list[i])));
1676 0 : if (merge(*set, group, dom, 2 * i + 1) < 0)
1677 0 : goto error;
1678 : }
1679 :
1680 0 : g = 0;
1681 0 : for (i = 0; i < 2 * n; ++i)
1682 0 : if (group[i] == i) {
1683 0 : if (g != i) {
1684 0 : (*set)[g] = (*set)[i];
1685 0 : (*set)[i] = NULL;
1686 : }
1687 0 : group[i] = g++;
1688 : } else
1689 0 : group[i] = group[group[i]];
1690 :
1691 0 : *n_group = g;
1692 :
1693 0 : return group;
1694 : error:
1695 0 : if (*set) {
1696 0 : for (i = 0; i < 2 * n; ++i)
1697 0 : isl_set_free((*set)[i]);
1698 0 : free(*set);
1699 0 : *set = NULL;
1700 : }
1701 0 : free(group);
1702 0 : return NULL;
1703 : }
1704 :
1705 : /* Check if the domains and ranges of the basic maps in "map" can
1706 : * be partitioned, and if so, apply Floyd-Warshall on the elements
1707 : * of the partition. Note that we also apply this algorithm
1708 : * if we want to compute the power, i.e., when "project" is not set.
1709 : * However, the results are unlikely to be exact since the recursive
1710 : * calls inside the Floyd-Warshall algorithm typically result in
1711 : * non-linear path lengths quite quickly.
1712 : */
1713 0 : static __isl_give isl_map *floyd_warshall(__isl_take isl_space *dim,
1714 : __isl_keep isl_map *map, int *exact, int project)
1715 : {
1716 : int i;
1717 0 : isl_set **set = NULL;
1718 0 : int *group = NULL;
1719 : int n;
1720 :
1721 0 : if (!map)
1722 0 : goto error;
1723 0 : if (map->n <= 1)
1724 0 : return incremental_closure(dim, map, exact, project);
1725 :
1726 0 : group = setup_groups(map->ctx, map->p, map->n, &set, &n);
1727 0 : if (!group)
1728 0 : goto error;
1729 :
1730 0 : for (i = 0; i < 2 * map->n; ++i)
1731 0 : isl_set_free(set[i]);
1732 :
1733 0 : free(set);
1734 :
1735 0 : return floyd_warshall_with_groups(dim, map, exact, project, group, n);
1736 : error:
1737 0 : isl_space_free(dim);
1738 0 : return NULL;
1739 : }
1740 :
1741 : /* Structure for representing the nodes of the graph of which
1742 : * strongly connected components are being computed.
1743 : *
1744 : * list contains the actual nodes
1745 : * check_closed is set if we may have used the fact that
1746 : * a pair of basic maps can be interchanged
1747 : */
1748 : struct isl_tc_follows_data {
1749 : isl_basic_map **list;
1750 : int check_closed;
1751 : };
1752 :
1753 : /* Check whether in the computation of the transitive closure
1754 : * "list[i]" (R_1) should follow (or be part of the same component as)
1755 : * "list[j]" (R_2).
1756 : *
1757 : * That is check whether
1758 : *
1759 : * R_1 \circ R_2
1760 : *
1761 : * is a subset of
1762 : *
1763 : * R_2 \circ R_1
1764 : *
1765 : * If so, then there is no reason for R_1 to immediately follow R_2
1766 : * in any path.
1767 : *
1768 : * *check_closed is set if the subset relation holds while
1769 : * R_1 \circ R_2 is not empty.
1770 : */
1771 0 : static isl_bool basic_map_follows(int i, int j, void *user)
1772 : {
1773 0 : struct isl_tc_follows_data *data = user;
1774 0 : struct isl_map *map12 = NULL;
1775 0 : struct isl_map *map21 = NULL;
1776 : isl_bool subset;
1777 :
1778 0 : if (!isl_space_tuple_is_equal(data->list[i]->dim, isl_dim_in,
1779 0 : data->list[j]->dim, isl_dim_out))
1780 0 : return isl_bool_false;
1781 :
1782 0 : map21 = isl_map_from_basic_map(
1783 : isl_basic_map_apply_range(
1784 0 : isl_basic_map_copy(data->list[j]),
1785 0 : isl_basic_map_copy(data->list[i])));
1786 0 : subset = isl_map_is_empty(map21);
1787 0 : if (subset < 0)
1788 0 : goto error;
1789 0 : if (subset) {
1790 0 : isl_map_free(map21);
1791 0 : return isl_bool_false;
1792 : }
1793 :
1794 0 : if (!isl_space_tuple_is_equal(data->list[i]->dim, isl_dim_in,
1795 0 : data->list[i]->dim, isl_dim_out) ||
1796 0 : !isl_space_tuple_is_equal(data->list[j]->dim, isl_dim_in,
1797 0 : data->list[j]->dim, isl_dim_out)) {
1798 0 : isl_map_free(map21);
1799 0 : return isl_bool_true;
1800 : }
1801 :
1802 0 : map12 = isl_map_from_basic_map(
1803 : isl_basic_map_apply_range(
1804 0 : isl_basic_map_copy(data->list[i]),
1805 0 : isl_basic_map_copy(data->list[j])));
1806 :
1807 0 : subset = isl_map_is_subset(map21, map12);
1808 :
1809 0 : isl_map_free(map12);
1810 0 : isl_map_free(map21);
1811 :
1812 0 : if (subset)
1813 0 : data->check_closed = 1;
1814 :
1815 0 : return subset < 0 ? isl_bool_error : !subset;
1816 : error:
1817 0 : isl_map_free(map21);
1818 0 : return isl_bool_error;
1819 : }
1820 :
1821 : /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
1822 : * and a dimension specification (Z^{n+1} -> Z^{n+1}),
1823 : * construct a map that is an overapproximation of the map
1824 : * that takes an element from the dom R \times Z to an
1825 : * element from ran R \times Z, such that the first n coordinates of the
1826 : * difference between them is a sum of differences between images
1827 : * and pre-images in one of the R_i and such that the last coordinate
1828 : * is equal to the number of steps taken.
1829 : * If "project" is set, then these final coordinates are not included,
1830 : * i.e., a relation of type Z^n -> Z^n is returned.
1831 : * That is, let
1832 : *
1833 : * \Delta_i = { y - x | (x, y) in R_i }
1834 : *
1835 : * then the constructed map is an overapproximation of
1836 : *
1837 : * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
1838 : * d = (\sum_i k_i \delta_i, \sum_i k_i) and
1839 : * x in dom R and x + d in ran R }
1840 : *
1841 : * or
1842 : *
1843 : * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
1844 : * d = (\sum_i k_i \delta_i) and
1845 : * x in dom R and x + d in ran R }
1846 : *
1847 : * if "project" is set.
1848 : *
1849 : * We first split the map into strongly connected components, perform
1850 : * the above on each component and then join the results in the correct
1851 : * order, at each join also taking in the union of both arguments
1852 : * to allow for paths that do not go through one of the two arguments.
1853 : */
1854 0 : static __isl_give isl_map *construct_power_components(__isl_take isl_space *dim,
1855 : __isl_keep isl_map *map, int *exact, int project)
1856 : {
1857 : int i, n, c;
1858 0 : struct isl_map *path = NULL;
1859 : struct isl_tc_follows_data data;
1860 0 : struct isl_tarjan_graph *g = NULL;
1861 : int *orig_exact;
1862 : int local_exact;
1863 :
1864 0 : if (!map)
1865 0 : goto error;
1866 0 : if (map->n <= 1)
1867 0 : return floyd_warshall(dim, map, exact, project);
1868 :
1869 0 : data.list = map->p;
1870 0 : data.check_closed = 0;
1871 0 : g = isl_tarjan_graph_init(map->ctx, map->n, &basic_map_follows, &data);
1872 0 : if (!g)
1873 0 : goto error;
1874 :
1875 0 : orig_exact = exact;
1876 0 : if (data.check_closed && !exact)
1877 0 : exact = &local_exact;
1878 :
1879 0 : c = 0;
1880 0 : i = 0;
1881 0 : n = map->n;
1882 0 : if (project)
1883 0 : path = isl_map_empty(isl_map_get_space(map));
1884 : else
1885 0 : path = isl_map_empty(isl_space_copy(dim));
1886 0 : path = anonymize(path);
1887 0 : while (n) {
1888 : struct isl_map *comp;
1889 : isl_map *path_comp, *path_comb;
1890 0 : comp = isl_map_alloc_space(isl_map_get_space(map), n, 0);
1891 0 : while (g->order[i] != -1) {
1892 0 : comp = isl_map_add_basic_map(comp,
1893 0 : isl_basic_map_copy(map->p[g->order[i]]));
1894 0 : --n;
1895 0 : ++i;
1896 : }
1897 0 : path_comp = floyd_warshall(isl_space_copy(dim),
1898 : comp, exact, project);
1899 0 : path_comp = anonymize(path_comp);
1900 0 : path_comb = isl_map_apply_range(isl_map_copy(path),
1901 : isl_map_copy(path_comp));
1902 0 : path = isl_map_union(path, path_comp);
1903 0 : path = isl_map_union(path, path_comb);
1904 0 : isl_map_free(comp);
1905 0 : ++i;
1906 0 : ++c;
1907 : }
1908 :
1909 0 : if (c > 1 && data.check_closed && !*exact) {
1910 : int closed;
1911 :
1912 0 : closed = isl_map_is_transitively_closed(path);
1913 0 : if (closed < 0)
1914 0 : goto error;
1915 0 : if (!closed) {
1916 0 : isl_tarjan_graph_free(g);
1917 0 : isl_map_free(path);
1918 0 : return floyd_warshall(dim, map, orig_exact, project);
1919 : }
1920 : }
1921 :
1922 0 : isl_tarjan_graph_free(g);
1923 0 : isl_space_free(dim);
1924 :
1925 0 : return path;
1926 : error:
1927 0 : isl_tarjan_graph_free(g);
1928 0 : isl_space_free(dim);
1929 0 : isl_map_free(path);
1930 0 : return NULL;
1931 : }
1932 :
1933 : /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D,
1934 : * construct a map that is an overapproximation of the map
1935 : * that takes an element from the space D to another
1936 : * element from the same space, such that the difference between
1937 : * them is a strictly positive sum of differences between images
1938 : * and pre-images in one of the R_i.
1939 : * The number of differences in the sum is equated to parameter "param".
1940 : * That is, let
1941 : *
1942 : * \Delta_i = { y - x | (x, y) in R_i }
1943 : *
1944 : * then the constructed map is an overapproximation of
1945 : *
1946 : * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
1947 : * d = \sum_i k_i \delta_i and k = \sum_i k_i > 0 }
1948 : * or
1949 : *
1950 : * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
1951 : * d = \sum_i k_i \delta_i and \sum_i k_i > 0 }
1952 : *
1953 : * if "project" is set.
1954 : *
1955 : * If "project" is not set, then
1956 : * we construct an extended mapping with an extra coordinate
1957 : * that indicates the number of steps taken. In particular,
1958 : * the difference in the last coordinate is equal to the number
1959 : * of steps taken to move from a domain element to the corresponding
1960 : * image element(s).
1961 : */
1962 0 : static __isl_give isl_map *construct_power(__isl_keep isl_map *map,
1963 : int *exact, int project)
1964 : {
1965 0 : struct isl_map *app = NULL;
1966 0 : isl_space *dim = NULL;
1967 :
1968 0 : if (!map)
1969 0 : return NULL;
1970 :
1971 0 : dim = isl_map_get_space(map);
1972 :
1973 0 : dim = isl_space_add_dims(dim, isl_dim_in, 1);
1974 0 : dim = isl_space_add_dims(dim, isl_dim_out, 1);
1975 :
1976 0 : app = construct_power_components(isl_space_copy(dim), map,
1977 : exact, project);
1978 :
1979 0 : isl_space_free(dim);
1980 :
1981 0 : return app;
1982 : }
1983 :
1984 : /* Compute the positive powers of "map", or an overapproximation.
1985 : * If the result is exact, then *exact is set to 1.
1986 : *
1987 : * If project is set, then we are actually interested in the transitive
1988 : * closure, so we can use a more relaxed exactness check.
1989 : * The lengths of the paths are also projected out instead of being
1990 : * encoded as the difference between an extra pair of final coordinates.
1991 : */
1992 0 : static __isl_give isl_map *map_power(__isl_take isl_map *map,
1993 : int *exact, int project)
1994 : {
1995 0 : struct isl_map *app = NULL;
1996 :
1997 0 : if (exact)
1998 0 : *exact = 1;
1999 :
2000 0 : if (!map)
2001 0 : return NULL;
2002 :
2003 0 : isl_assert(map->ctx,
2004 : isl_map_dim(map, isl_dim_in) == isl_map_dim(map, isl_dim_out),
2005 : goto error);
2006 :
2007 0 : app = construct_power(map, exact, project);
2008 :
2009 0 : isl_map_free(map);
2010 0 : return app;
2011 : error:
2012 0 : isl_map_free(map);
2013 0 : isl_map_free(app);
2014 0 : return NULL;
2015 : }
2016 :
2017 : /* Compute the positive powers of "map", or an overapproximation.
2018 : * The result maps the exponent to a nested copy of the corresponding power.
2019 : * If the result is exact, then *exact is set to 1.
2020 : * map_power constructs an extended relation with the path lengths
2021 : * encoded as the difference between the final coordinates.
2022 : * In the final step, this difference is equated to an extra parameter
2023 : * and made positive. The extra coordinates are subsequently projected out
2024 : * and the parameter is turned into the domain of the result.
2025 : */
2026 0 : __isl_give isl_map *isl_map_power(__isl_take isl_map *map, int *exact)
2027 : {
2028 : isl_space *target_dim;
2029 : isl_space *dim;
2030 : isl_map *diff;
2031 : unsigned d;
2032 : unsigned param;
2033 :
2034 0 : if (!map)
2035 0 : return NULL;
2036 :
2037 0 : d = isl_map_dim(map, isl_dim_in);
2038 0 : param = isl_map_dim(map, isl_dim_param);
2039 :
2040 0 : map = isl_map_compute_divs(map);
2041 0 : map = isl_map_coalesce(map);
2042 :
2043 0 : if (isl_map_plain_is_empty(map)) {
2044 0 : map = isl_map_from_range(isl_map_wrap(map));
2045 0 : map = isl_map_add_dims(map, isl_dim_in, 1);
2046 0 : map = isl_map_set_dim_name(map, isl_dim_in, 0, "k");
2047 0 : return map;
2048 : }
2049 :
2050 0 : target_dim = isl_map_get_space(map);
2051 0 : target_dim = isl_space_from_range(isl_space_wrap(target_dim));
2052 0 : target_dim = isl_space_add_dims(target_dim, isl_dim_in, 1);
2053 0 : target_dim = isl_space_set_dim_name(target_dim, isl_dim_in, 0, "k");
2054 :
2055 0 : map = map_power(map, exact, 0);
2056 :
2057 0 : map = isl_map_add_dims(map, isl_dim_param, 1);
2058 0 : dim = isl_map_get_space(map);
2059 0 : diff = equate_parameter_to_length(dim, param);
2060 0 : map = isl_map_intersect(map, diff);
2061 0 : map = isl_map_project_out(map, isl_dim_in, d, 1);
2062 0 : map = isl_map_project_out(map, isl_dim_out, d, 1);
2063 0 : map = isl_map_from_range(isl_map_wrap(map));
2064 0 : map = isl_map_move_dims(map, isl_dim_in, 0, isl_dim_param, param, 1);
2065 :
2066 0 : map = isl_map_reset_space(map, target_dim);
2067 :
2068 0 : return map;
2069 : }
2070 :
2071 : /* Compute a relation that maps each element in the range of the input
2072 : * relation to the lengths of all paths composed of edges in the input
2073 : * relation that end up in the given range element.
2074 : * The result may be an overapproximation, in which case *exact is set to 0.
2075 : * The resulting relation is very similar to the power relation.
2076 : * The difference are that the domain has been projected out, the
2077 : * range has become the domain and the exponent is the range instead
2078 : * of a parameter.
2079 : */
2080 0 : __isl_give isl_map *isl_map_reaching_path_lengths(__isl_take isl_map *map,
2081 : int *exact)
2082 : {
2083 : isl_space *dim;
2084 : isl_map *diff;
2085 : unsigned d;
2086 : unsigned param;
2087 :
2088 0 : if (!map)
2089 0 : return NULL;
2090 :
2091 0 : d = isl_map_dim(map, isl_dim_in);
2092 0 : param = isl_map_dim(map, isl_dim_param);
2093 :
2094 0 : map = isl_map_compute_divs(map);
2095 0 : map = isl_map_coalesce(map);
2096 :
2097 0 : if (isl_map_plain_is_empty(map)) {
2098 0 : if (exact)
2099 0 : *exact = 1;
2100 0 : map = isl_map_project_out(map, isl_dim_out, 0, d);
2101 0 : map = isl_map_add_dims(map, isl_dim_out, 1);
2102 0 : return map;
2103 : }
2104 :
2105 0 : map = map_power(map, exact, 0);
2106 :
2107 0 : map = isl_map_add_dims(map, isl_dim_param, 1);
2108 0 : dim = isl_map_get_space(map);
2109 0 : diff = equate_parameter_to_length(dim, param);
2110 0 : map = isl_map_intersect(map, diff);
2111 0 : map = isl_map_project_out(map, isl_dim_in, 0, d + 1);
2112 0 : map = isl_map_project_out(map, isl_dim_out, d, 1);
2113 0 : map = isl_map_reverse(map);
2114 0 : map = isl_map_move_dims(map, isl_dim_out, 0, isl_dim_param, param, 1);
2115 :
2116 0 : return map;
2117 : }
2118 :
2119 : /* Given a map, compute the smallest superset of this map that is of the form
2120 : *
2121 : * { i -> j : L <= j - i <= U and exists a_p: j_p - i_p = M_p a_p }
2122 : *
2123 : * (where p ranges over the (non-parametric) dimensions),
2124 : * compute the transitive closure of this map, i.e.,
2125 : *
2126 : * { i -> j : exists k > 0:
2127 : * k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p }
2128 : *
2129 : * and intersect domain and range of this transitive closure with
2130 : * the given domain and range.
2131 : *
2132 : * If with_id is set, then try to include as much of the identity mapping
2133 : * as possible, by computing
2134 : *
2135 : * { i -> j : exists k >= 0:
2136 : * k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p }
2137 : *
2138 : * instead (i.e., allow k = 0).
2139 : *
2140 : * In practice, we compute the difference set
2141 : *
2142 : * delta = { j - i | i -> j in map },
2143 : *
2144 : * look for stride constraint on the individual dimensions and compute
2145 : * (constant) lower and upper bounds for each individual dimension,
2146 : * adding a constraint for each bound not equal to infinity.
2147 : */
2148 0 : static __isl_give isl_map *box_closure_on_domain(__isl_take isl_map *map,
2149 : __isl_take isl_set *dom, __isl_take isl_set *ran, int with_id)
2150 : {
2151 : int i;
2152 : int k;
2153 : unsigned d;
2154 : unsigned nparam;
2155 : unsigned total;
2156 : isl_space *dim;
2157 : isl_set *delta;
2158 0 : isl_map *app = NULL;
2159 0 : isl_basic_set *aff = NULL;
2160 0 : isl_basic_map *bmap = NULL;
2161 0 : isl_vec *obj = NULL;
2162 : isl_int opt;
2163 :
2164 0 : isl_int_init(opt);
2165 :
2166 0 : delta = isl_map_deltas(isl_map_copy(map));
2167 :
2168 0 : aff = isl_set_affine_hull(isl_set_copy(delta));
2169 0 : if (!aff)
2170 0 : goto error;
2171 0 : dim = isl_map_get_space(map);
2172 0 : d = isl_space_dim(dim, isl_dim_in);
2173 0 : nparam = isl_space_dim(dim, isl_dim_param);
2174 0 : total = isl_space_dim(dim, isl_dim_all);
2175 0 : bmap = isl_basic_map_alloc_space(dim,
2176 0 : aff->n_div + 1, aff->n_div, 2 * d + 1);
2177 0 : for (i = 0; i < aff->n_div + 1; ++i) {
2178 0 : k = isl_basic_map_alloc_div(bmap);
2179 0 : if (k < 0)
2180 0 : goto error;
2181 0 : isl_int_set_si(bmap->div[k][0], 0);
2182 : }
2183 0 : for (i = 0; i < aff->n_eq; ++i) {
2184 0 : if (!isl_basic_set_eq_is_stride(aff, i))
2185 0 : continue;
2186 0 : k = isl_basic_map_alloc_equality(bmap);
2187 0 : if (k < 0)
2188 0 : goto error;
2189 0 : isl_seq_clr(bmap->eq[k], 1 + nparam);
2190 0 : isl_seq_cpy(bmap->eq[k] + 1 + nparam + d,
2191 0 : aff->eq[i] + 1 + nparam, d);
2192 0 : isl_seq_neg(bmap->eq[k] + 1 + nparam,
2193 0 : aff->eq[i] + 1 + nparam, d);
2194 0 : isl_seq_cpy(bmap->eq[k] + 1 + nparam + 2 * d,
2195 0 : aff->eq[i] + 1 + nparam + d, aff->n_div);
2196 0 : isl_int_set_si(bmap->eq[k][1 + total + aff->n_div], 0);
2197 : }
2198 0 : obj = isl_vec_alloc(map->ctx, 1 + nparam + d);
2199 0 : if (!obj)
2200 0 : goto error;
2201 0 : isl_seq_clr(obj->el, 1 + nparam + d);
2202 0 : for (i = 0; i < d; ++ i) {
2203 : enum isl_lp_result res;
2204 :
2205 0 : isl_int_set_si(obj->el[1 + nparam + i], 1);
2206 :
2207 0 : res = isl_set_solve_lp(delta, 0, obj->el, map->ctx->one, &opt,
2208 : NULL, NULL);
2209 0 : if (res == isl_lp_error)
2210 0 : goto error;
2211 0 : if (res == isl_lp_ok) {
2212 0 : k = isl_basic_map_alloc_inequality(bmap);
2213 0 : if (k < 0)
2214 0 : goto error;
2215 0 : isl_seq_clr(bmap->ineq[k],
2216 0 : 1 + nparam + 2 * d + bmap->n_div);
2217 0 : isl_int_set_si(bmap->ineq[k][1 + nparam + i], -1);
2218 0 : isl_int_set_si(bmap->ineq[k][1 + nparam + d + i], 1);
2219 0 : isl_int_neg(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], opt);
2220 : }
2221 :
2222 0 : res = isl_set_solve_lp(delta, 1, obj->el, map->ctx->one, &opt,
2223 : NULL, NULL);
2224 0 : if (res == isl_lp_error)
2225 0 : goto error;
2226 0 : if (res == isl_lp_ok) {
2227 0 : k = isl_basic_map_alloc_inequality(bmap);
2228 0 : if (k < 0)
2229 0 : goto error;
2230 0 : isl_seq_clr(bmap->ineq[k],
2231 0 : 1 + nparam + 2 * d + bmap->n_div);
2232 0 : isl_int_set_si(bmap->ineq[k][1 + nparam + i], 1);
2233 0 : isl_int_set_si(bmap->ineq[k][1 + nparam + d + i], -1);
2234 0 : isl_int_set(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], opt);
2235 : }
2236 :
2237 0 : isl_int_set_si(obj->el[1 + nparam + i], 0);
2238 : }
2239 0 : k = isl_basic_map_alloc_inequality(bmap);
2240 0 : if (k < 0)
2241 0 : goto error;
2242 0 : isl_seq_clr(bmap->ineq[k],
2243 0 : 1 + nparam + 2 * d + bmap->n_div);
2244 0 : if (!with_id)
2245 0 : isl_int_set_si(bmap->ineq[k][0], -1);
2246 0 : isl_int_set_si(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], 1);
2247 :
2248 0 : app = isl_map_from_domain_and_range(dom, ran);
2249 :
2250 0 : isl_vec_free(obj);
2251 0 : isl_basic_set_free(aff);
2252 0 : isl_map_free(map);
2253 0 : bmap = isl_basic_map_finalize(bmap);
2254 0 : isl_set_free(delta);
2255 0 : isl_int_clear(opt);
2256 :
2257 0 : map = isl_map_from_basic_map(bmap);
2258 0 : map = isl_map_intersect(map, app);
2259 :
2260 0 : return map;
2261 : error:
2262 0 : isl_vec_free(obj);
2263 0 : isl_basic_map_free(bmap);
2264 0 : isl_basic_set_free(aff);
2265 0 : isl_set_free(dom);
2266 0 : isl_set_free(ran);
2267 0 : isl_map_free(map);
2268 0 : isl_set_free(delta);
2269 0 : isl_int_clear(opt);
2270 0 : return NULL;
2271 : }
2272 :
2273 : /* Given a map, compute the smallest superset of this map that is of the form
2274 : *
2275 : * { i -> j : L <= j - i <= U and exists a_p: j_p - i_p = M_p a_p }
2276 : *
2277 : * (where p ranges over the (non-parametric) dimensions),
2278 : * compute the transitive closure of this map, i.e.,
2279 : *
2280 : * { i -> j : exists k > 0:
2281 : * k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p }
2282 : *
2283 : * and intersect domain and range of this transitive closure with
2284 : * domain and range of the original map.
2285 : */
2286 0 : static __isl_give isl_map *box_closure(__isl_take isl_map *map)
2287 : {
2288 : isl_set *domain;
2289 : isl_set *range;
2290 :
2291 0 : domain = isl_map_domain(isl_map_copy(map));
2292 0 : domain = isl_set_coalesce(domain);
2293 0 : range = isl_map_range(isl_map_copy(map));
2294 0 : range = isl_set_coalesce(range);
2295 :
2296 0 : return box_closure_on_domain(map, domain, range, 0);
2297 : }
2298 :
2299 : /* Given a map, compute the smallest superset of this map that is of the form
2300 : *
2301 : * { i -> j : L <= j - i <= U and exists a_p: j_p - i_p = M_p a_p }
2302 : *
2303 : * (where p ranges over the (non-parametric) dimensions),
2304 : * compute the transitive and partially reflexive closure of this map, i.e.,
2305 : *
2306 : * { i -> j : exists k >= 0:
2307 : * k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p }
2308 : *
2309 : * and intersect domain and range of this transitive closure with
2310 : * the given domain.
2311 : */
2312 0 : static __isl_give isl_map *box_closure_with_identity(__isl_take isl_map *map,
2313 : __isl_take isl_set *dom)
2314 : {
2315 0 : return box_closure_on_domain(map, dom, isl_set_copy(dom), 1);
2316 : }
2317 :
2318 : /* Check whether app is the transitive closure of map.
2319 : * In particular, check that app is acyclic and, if so,
2320 : * check that
2321 : *
2322 : * app \subset (map \cup (map \circ app))
2323 : */
2324 0 : static int check_exactness_omega(__isl_keep isl_map *map,
2325 : __isl_keep isl_map *app)
2326 : {
2327 : isl_set *delta;
2328 : int i;
2329 : int is_empty, is_exact;
2330 : unsigned d;
2331 : isl_map *test;
2332 :
2333 0 : delta = isl_map_deltas(isl_map_copy(app));
2334 0 : d = isl_set_dim(delta, isl_dim_set);
2335 0 : for (i = 0; i < d; ++i)
2336 0 : delta = isl_set_fix_si(delta, isl_dim_set, i, 0);
2337 0 : is_empty = isl_set_is_empty(delta);
2338 0 : isl_set_free(delta);
2339 0 : if (is_empty < 0)
2340 0 : return -1;
2341 0 : if (!is_empty)
2342 0 : return 0;
2343 :
2344 0 : test = isl_map_apply_range(isl_map_copy(app), isl_map_copy(map));
2345 0 : test = isl_map_union(test, isl_map_copy(map));
2346 0 : is_exact = isl_map_is_subset(app, test);
2347 0 : isl_map_free(test);
2348 :
2349 0 : return is_exact;
2350 : }
2351 :
2352 : /* Check if basic map M_i can be combined with all the other
2353 : * basic maps such that
2354 : *
2355 : * (\cup_j M_j)^+
2356 : *
2357 : * can be computed as
2358 : *
2359 : * M_i \cup (\cup_{j \ne i} M_i^* \circ M_j \circ M_i^*)^+
2360 : *
2361 : * In particular, check if we can compute a compact representation
2362 : * of
2363 : *
2364 : * M_i^* \circ M_j \circ M_i^*
2365 : *
2366 : * for each j != i.
2367 : * Let M_i^? be an extension of M_i^+ that allows paths
2368 : * of length zero, i.e., the result of box_closure(., 1).
2369 : * The criterion, as proposed by Kelly et al., is that
2370 : * id = M_i^? - M_i^+ can be represented as a basic map
2371 : * and that
2372 : *
2373 : * id \circ M_j \circ id = M_j
2374 : *
2375 : * for each j != i.
2376 : *
2377 : * If this function returns 1, then tc and qc are set to
2378 : * M_i^+ and M_i^?, respectively.
2379 : */
2380 0 : static int can_be_split_off(__isl_keep isl_map *map, int i,
2381 : __isl_give isl_map **tc, __isl_give isl_map **qc)
2382 : {
2383 0 : isl_map *map_i, *id = NULL;
2384 0 : int j = -1;
2385 : isl_set *C;
2386 :
2387 0 : *tc = NULL;
2388 0 : *qc = NULL;
2389 :
2390 0 : C = isl_set_union(isl_map_domain(isl_map_copy(map)),
2391 : isl_map_range(isl_map_copy(map)));
2392 0 : C = isl_set_from_basic_set(isl_set_simple_hull(C));
2393 0 : if (!C)
2394 0 : goto error;
2395 :
2396 0 : map_i = isl_map_from_basic_map(isl_basic_map_copy(map->p[i]));
2397 0 : *tc = box_closure(isl_map_copy(map_i));
2398 0 : *qc = box_closure_with_identity(map_i, C);
2399 0 : id = isl_map_subtract(isl_map_copy(*qc), isl_map_copy(*tc));
2400 :
2401 0 : if (!id || !*qc)
2402 : goto error;
2403 0 : if (id->n != 1 || (*qc)->n != 1)
2404 : goto done;
2405 :
2406 0 : for (j = 0; j < map->n; ++j) {
2407 : isl_map *map_j, *test;
2408 : int is_ok;
2409 :
2410 0 : if (i == j)
2411 0 : continue;
2412 0 : map_j = isl_map_from_basic_map(
2413 0 : isl_basic_map_copy(map->p[j]));
2414 0 : test = isl_map_apply_range(isl_map_copy(id),
2415 : isl_map_copy(map_j));
2416 0 : test = isl_map_apply_range(test, isl_map_copy(id));
2417 0 : is_ok = isl_map_is_equal(test, map_j);
2418 0 : isl_map_free(map_j);
2419 0 : isl_map_free(test);
2420 0 : if (is_ok < 0)
2421 0 : goto error;
2422 0 : if (!is_ok)
2423 0 : break;
2424 : }
2425 :
2426 : done:
2427 0 : isl_map_free(id);
2428 0 : if (j == map->n)
2429 0 : return 1;
2430 :
2431 0 : isl_map_free(*qc);
2432 0 : isl_map_free(*tc);
2433 0 : *qc = NULL;
2434 0 : *tc = NULL;
2435 :
2436 0 : return 0;
2437 : error:
2438 0 : isl_map_free(id);
2439 0 : isl_map_free(*qc);
2440 0 : isl_map_free(*tc);
2441 0 : *qc = NULL;
2442 0 : *tc = NULL;
2443 0 : return -1;
2444 : }
2445 :
2446 0 : static __isl_give isl_map *box_closure_with_check(__isl_take isl_map *map,
2447 : int *exact)
2448 : {
2449 : isl_map *app;
2450 :
2451 0 : app = box_closure(isl_map_copy(map));
2452 0 : if (exact)
2453 0 : *exact = check_exactness_omega(map, app);
2454 :
2455 0 : isl_map_free(map);
2456 0 : return app;
2457 : }
2458 :
2459 : /* Compute an overapproximation of the transitive closure of "map"
2460 : * using a variation of the algorithm from
2461 : * "Transitive Closure of Infinite Graphs and its Applications"
2462 : * by Kelly et al.
2463 : *
2464 : * We first check whether we can can split of any basic map M_i and
2465 : * compute
2466 : *
2467 : * (\cup_j M_j)^+
2468 : *
2469 : * as
2470 : *
2471 : * M_i \cup (\cup_{j \ne i} M_i^* \circ M_j \circ M_i^*)^+
2472 : *
2473 : * using a recursive call on the remaining map.
2474 : *
2475 : * If not, we simply call box_closure on the whole map.
2476 : */
2477 0 : static __isl_give isl_map *transitive_closure_omega(__isl_take isl_map *map,
2478 : int *exact)
2479 : {
2480 : int i, j;
2481 : int exact_i;
2482 : isl_map *app;
2483 :
2484 0 : if (!map)
2485 0 : return NULL;
2486 0 : if (map->n == 1)
2487 0 : return box_closure_with_check(map, exact);
2488 :
2489 0 : for (i = 0; i < map->n; ++i) {
2490 : int ok;
2491 : isl_map *qc, *tc;
2492 0 : ok = can_be_split_off(map, i, &tc, &qc);
2493 0 : if (ok < 0)
2494 0 : goto error;
2495 0 : if (!ok)
2496 0 : continue;
2497 :
2498 0 : app = isl_map_alloc_space(isl_map_get_space(map), map->n - 1, 0);
2499 :
2500 0 : for (j = 0; j < map->n; ++j) {
2501 0 : if (j == i)
2502 0 : continue;
2503 0 : app = isl_map_add_basic_map(app,
2504 0 : isl_basic_map_copy(map->p[j]));
2505 : }
2506 :
2507 0 : app = isl_map_apply_range(isl_map_copy(qc), app);
2508 0 : app = isl_map_apply_range(app, qc);
2509 :
2510 0 : app = isl_map_union(tc, transitive_closure_omega(app, NULL));
2511 0 : exact_i = check_exactness_omega(map, app);
2512 0 : if (exact_i == 1) {
2513 0 : if (exact)
2514 0 : *exact = exact_i;
2515 0 : isl_map_free(map);
2516 0 : return app;
2517 : }
2518 0 : isl_map_free(app);
2519 0 : if (exact_i < 0)
2520 0 : goto error;
2521 : }
2522 :
2523 0 : return box_closure_with_check(map, exact);
2524 : error:
2525 0 : isl_map_free(map);
2526 0 : return NULL;
2527 : }
2528 :
2529 : /* Compute the transitive closure of "map", or an overapproximation.
2530 : * If the result is exact, then *exact is set to 1.
2531 : * Simply use map_power to compute the powers of map, but tell
2532 : * it to project out the lengths of the paths instead of equating
2533 : * the length to a parameter.
2534 : */
2535 0 : __isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
2536 : int *exact)
2537 : {
2538 : isl_space *target_dim;
2539 : int closed;
2540 :
2541 0 : if (!map)
2542 0 : goto error;
2543 :
2544 0 : if (map->ctx->opt->closure == ISL_CLOSURE_BOX)
2545 0 : return transitive_closure_omega(map, exact);
2546 :
2547 0 : map = isl_map_compute_divs(map);
2548 0 : map = isl_map_coalesce(map);
2549 0 : closed = isl_map_is_transitively_closed(map);
2550 0 : if (closed < 0)
2551 0 : goto error;
2552 0 : if (closed) {
2553 0 : if (exact)
2554 0 : *exact = 1;
2555 0 : return map;
2556 : }
2557 :
2558 0 : target_dim = isl_map_get_space(map);
2559 0 : map = map_power(map, exact, 1);
2560 0 : map = isl_map_reset_space(map, target_dim);
2561 :
2562 0 : return map;
2563 : error:
2564 0 : isl_map_free(map);
2565 0 : return NULL;
2566 : }
2567 :
2568 0 : static isl_stat inc_count(__isl_take isl_map *map, void *user)
2569 : {
2570 0 : int *n = user;
2571 :
2572 0 : *n += map->n;
2573 :
2574 0 : isl_map_free(map);
2575 :
2576 0 : return isl_stat_ok;
2577 : }
2578 :
2579 0 : static isl_stat collect_basic_map(__isl_take isl_map *map, void *user)
2580 : {
2581 : int i;
2582 0 : isl_basic_map ***next = user;
2583 :
2584 0 : for (i = 0; i < map->n; ++i) {
2585 0 : **next = isl_basic_map_copy(map->p[i]);
2586 0 : if (!**next)
2587 0 : goto error;
2588 0 : (*next)++;
2589 : }
2590 :
2591 0 : isl_map_free(map);
2592 0 : return isl_stat_ok;
2593 : error:
2594 0 : isl_map_free(map);
2595 0 : return isl_stat_error;
2596 : }
2597 :
2598 : /* Perform Floyd-Warshall on the given list of basic relations.
2599 : * The basic relations may live in different dimensions,
2600 : * but basic relations that get assigned to the diagonal of the
2601 : * grid have domains and ranges of the same dimension and so
2602 : * the standard algorithm can be used because the nested transitive
2603 : * closures are only applied to diagonal elements and because all
2604 : * compositions are peformed on relations with compatible domains and ranges.
2605 : */
2606 0 : static __isl_give isl_union_map *union_floyd_warshall_on_list(isl_ctx *ctx,
2607 : __isl_keep isl_basic_map **list, int n, int *exact)
2608 : {
2609 : int i, j, k;
2610 : int n_group;
2611 0 : int *group = NULL;
2612 0 : isl_set **set = NULL;
2613 0 : isl_map ***grid = NULL;
2614 : isl_union_map *app;
2615 :
2616 0 : group = setup_groups(ctx, list, n, &set, &n_group);
2617 0 : if (!group)
2618 0 : goto error;
2619 :
2620 0 : grid = isl_calloc_array(ctx, isl_map **, n_group);
2621 0 : if (!grid)
2622 0 : goto error;
2623 0 : for (i = 0; i < n_group; ++i) {
2624 0 : grid[i] = isl_calloc_array(ctx, isl_map *, n_group);
2625 0 : if (!grid[i])
2626 0 : goto error;
2627 0 : for (j = 0; j < n_group; ++j) {
2628 : isl_space *dim1, *dim2, *dim;
2629 0 : dim1 = isl_space_reverse(isl_set_get_space(set[i]));
2630 0 : dim2 = isl_set_get_space(set[j]);
2631 0 : dim = isl_space_join(dim1, dim2);
2632 0 : grid[i][j] = isl_map_empty(dim);
2633 : }
2634 : }
2635 :
2636 0 : for (k = 0; k < n; ++k) {
2637 0 : i = group[2 * k];
2638 0 : j = group[2 * k + 1];
2639 0 : grid[i][j] = isl_map_union(grid[i][j],
2640 : isl_map_from_basic_map(
2641 0 : isl_basic_map_copy(list[k])));
2642 : }
2643 :
2644 0 : floyd_warshall_iterate(grid, n_group, exact);
2645 :
2646 0 : app = isl_union_map_empty(isl_map_get_space(grid[0][0]));
2647 :
2648 0 : for (i = 0; i < n_group; ++i) {
2649 0 : for (j = 0; j < n_group; ++j)
2650 0 : app = isl_union_map_add_map(app, grid[i][j]);
2651 0 : free(grid[i]);
2652 : }
2653 0 : free(grid);
2654 :
2655 0 : for (i = 0; i < 2 * n; ++i)
2656 0 : isl_set_free(set[i]);
2657 0 : free(set);
2658 :
2659 0 : free(group);
2660 0 : return app;
2661 : error:
2662 0 : if (grid)
2663 0 : for (i = 0; i < n_group; ++i) {
2664 0 : if (!grid[i])
2665 0 : continue;
2666 0 : for (j = 0; j < n_group; ++j)
2667 0 : isl_map_free(grid[i][j]);
2668 0 : free(grid[i]);
2669 : }
2670 0 : free(grid);
2671 0 : if (set) {
2672 0 : for (i = 0; i < 2 * n; ++i)
2673 0 : isl_set_free(set[i]);
2674 0 : free(set);
2675 : }
2676 0 : free(group);
2677 0 : return NULL;
2678 : }
2679 :
2680 : /* Perform Floyd-Warshall on the given union relation.
2681 : * The implementation is very similar to that for non-unions.
2682 : * The main difference is that it is applied unconditionally.
2683 : * We first extract a list of basic maps from the union map
2684 : * and then perform the algorithm on this list.
2685 : */
2686 0 : static __isl_give isl_union_map *union_floyd_warshall(
2687 : __isl_take isl_union_map *umap, int *exact)
2688 : {
2689 : int i, n;
2690 : isl_ctx *ctx;
2691 0 : isl_basic_map **list = NULL;
2692 : isl_basic_map **next;
2693 : isl_union_map *res;
2694 :
2695 0 : n = 0;
2696 0 : if (isl_union_map_foreach_map(umap, inc_count, &n) < 0)
2697 0 : goto error;
2698 :
2699 0 : ctx = isl_union_map_get_ctx(umap);
2700 0 : list = isl_calloc_array(ctx, isl_basic_map *, n);
2701 0 : if (!list)
2702 0 : goto error;
2703 :
2704 0 : next = list;
2705 0 : if (isl_union_map_foreach_map(umap, collect_basic_map, &next) < 0)
2706 0 : goto error;
2707 :
2708 0 : res = union_floyd_warshall_on_list(ctx, list, n, exact);
2709 :
2710 0 : if (list) {
2711 0 : for (i = 0; i < n; ++i)
2712 0 : isl_basic_map_free(list[i]);
2713 0 : free(list);
2714 : }
2715 :
2716 0 : isl_union_map_free(umap);
2717 0 : return res;
2718 : error:
2719 0 : if (list) {
2720 0 : for (i = 0; i < n; ++i)
2721 0 : isl_basic_map_free(list[i]);
2722 0 : free(list);
2723 : }
2724 0 : isl_union_map_free(umap);
2725 0 : return NULL;
2726 : }
2727 :
2728 : /* Decompose the give union relation into strongly connected components.
2729 : * The implementation is essentially the same as that of
2730 : * construct_power_components with the major difference that all
2731 : * operations are performed on union maps.
2732 : */
2733 0 : static __isl_give isl_union_map *union_components(
2734 : __isl_take isl_union_map *umap, int *exact)
2735 : {
2736 : int i;
2737 : int n;
2738 : isl_ctx *ctx;
2739 0 : isl_basic_map **list = NULL;
2740 : isl_basic_map **next;
2741 0 : isl_union_map *path = NULL;
2742 : struct isl_tc_follows_data data;
2743 0 : struct isl_tarjan_graph *g = NULL;
2744 : int c, l;
2745 0 : int recheck = 0;
2746 :
2747 0 : n = 0;
2748 0 : if (isl_union_map_foreach_map(umap, inc_count, &n) < 0)
2749 0 : goto error;
2750 :
2751 0 : if (n == 0)
2752 0 : return umap;
2753 0 : if (n <= 1)
2754 0 : return union_floyd_warshall(umap, exact);
2755 :
2756 0 : ctx = isl_union_map_get_ctx(umap);
2757 0 : list = isl_calloc_array(ctx, isl_basic_map *, n);
2758 0 : if (!list)
2759 0 : goto error;
2760 :
2761 0 : next = list;
2762 0 : if (isl_union_map_foreach_map(umap, collect_basic_map, &next) < 0)
2763 0 : goto error;
2764 :
2765 0 : data.list = list;
2766 0 : data.check_closed = 0;
2767 0 : g = isl_tarjan_graph_init(ctx, n, &basic_map_follows, &data);
2768 0 : if (!g)
2769 0 : goto error;
2770 :
2771 0 : c = 0;
2772 0 : i = 0;
2773 0 : l = n;
2774 0 : path = isl_union_map_empty(isl_union_map_get_space(umap));
2775 0 : while (l) {
2776 : isl_union_map *comp;
2777 : isl_union_map *path_comp, *path_comb;
2778 0 : comp = isl_union_map_empty(isl_union_map_get_space(umap));
2779 0 : while (g->order[i] != -1) {
2780 0 : comp = isl_union_map_add_map(comp,
2781 : isl_map_from_basic_map(
2782 0 : isl_basic_map_copy(list[g->order[i]])));
2783 0 : --l;
2784 0 : ++i;
2785 : }
2786 0 : path_comp = union_floyd_warshall(comp, exact);
2787 0 : path_comb = isl_union_map_apply_range(isl_union_map_copy(path),
2788 : isl_union_map_copy(path_comp));
2789 0 : path = isl_union_map_union(path, path_comp);
2790 0 : path = isl_union_map_union(path, path_comb);
2791 0 : ++i;
2792 0 : ++c;
2793 : }
2794 :
2795 0 : if (c > 1 && data.check_closed && !*exact) {
2796 : int closed;
2797 :
2798 0 : closed = isl_union_map_is_transitively_closed(path);
2799 0 : if (closed < 0)
2800 0 : goto error;
2801 0 : recheck = !closed;
2802 : }
2803 :
2804 0 : isl_tarjan_graph_free(g);
2805 :
2806 0 : for (i = 0; i < n; ++i)
2807 0 : isl_basic_map_free(list[i]);
2808 0 : free(list);
2809 :
2810 0 : if (recheck) {
2811 0 : isl_union_map_free(path);
2812 0 : return union_floyd_warshall(umap, exact);
2813 : }
2814 :
2815 0 : isl_union_map_free(umap);
2816 :
2817 0 : return path;
2818 : error:
2819 0 : isl_tarjan_graph_free(g);
2820 0 : if (list) {
2821 0 : for (i = 0; i < n; ++i)
2822 0 : isl_basic_map_free(list[i]);
2823 0 : free(list);
2824 : }
2825 0 : isl_union_map_free(umap);
2826 0 : isl_union_map_free(path);
2827 0 : return NULL;
2828 : }
2829 :
2830 : /* Compute the transitive closure of "umap", or an overapproximation.
2831 : * If the result is exact, then *exact is set to 1.
2832 : */
2833 0 : __isl_give isl_union_map *isl_union_map_transitive_closure(
2834 : __isl_take isl_union_map *umap, int *exact)
2835 : {
2836 : int closed;
2837 :
2838 0 : if (!umap)
2839 0 : return NULL;
2840 :
2841 0 : if (exact)
2842 0 : *exact = 1;
2843 :
2844 0 : umap = isl_union_map_compute_divs(umap);
2845 0 : umap = isl_union_map_coalesce(umap);
2846 0 : closed = isl_union_map_is_transitively_closed(umap);
2847 0 : if (closed < 0)
2848 0 : goto error;
2849 0 : if (closed)
2850 0 : return umap;
2851 0 : umap = union_components(umap, exact);
2852 0 : return umap;
2853 : error:
2854 0 : isl_union_map_free(umap);
2855 0 : return NULL;
2856 : }
2857 :
2858 : struct isl_union_power {
2859 : isl_union_map *pow;
2860 : int *exact;
2861 : };
2862 :
2863 0 : static isl_stat power(__isl_take isl_map *map, void *user)
2864 : {
2865 0 : struct isl_union_power *up = user;
2866 :
2867 0 : map = isl_map_power(map, up->exact);
2868 0 : up->pow = isl_union_map_from_map(map);
2869 :
2870 0 : return isl_stat_error;
2871 : }
2872 :
2873 : /* Construct a map [x] -> [x+1], with parameters prescribed by "dim".
2874 : */
2875 0 : static __isl_give isl_union_map *increment(__isl_take isl_space *dim)
2876 : {
2877 : int k;
2878 : isl_basic_map *bmap;
2879 :
2880 0 : dim = isl_space_add_dims(dim, isl_dim_in, 1);
2881 0 : dim = isl_space_add_dims(dim, isl_dim_out, 1);
2882 0 : bmap = isl_basic_map_alloc_space(dim, 0, 1, 0);
2883 0 : k = isl_basic_map_alloc_equality(bmap);
2884 0 : if (k < 0)
2885 0 : goto error;
2886 0 : isl_seq_clr(bmap->eq[k], 1 + isl_basic_map_total_dim(bmap));
2887 0 : isl_int_set_si(bmap->eq[k][0], 1);
2888 0 : isl_int_set_si(bmap->eq[k][isl_basic_map_offset(bmap, isl_dim_in)], 1);
2889 0 : isl_int_set_si(bmap->eq[k][isl_basic_map_offset(bmap, isl_dim_out)], -1);
2890 0 : return isl_union_map_from_map(isl_map_from_basic_map(bmap));
2891 : error:
2892 0 : isl_basic_map_free(bmap);
2893 0 : return NULL;
2894 : }
2895 :
2896 : /* Construct a map [[x]->[y]] -> [y-x], with parameters prescribed by "dim".
2897 : */
2898 0 : static __isl_give isl_union_map *deltas_map(__isl_take isl_space *dim)
2899 : {
2900 : isl_basic_map *bmap;
2901 :
2902 0 : dim = isl_space_add_dims(dim, isl_dim_in, 1);
2903 0 : dim = isl_space_add_dims(dim, isl_dim_out, 1);
2904 0 : bmap = isl_basic_map_universe(dim);
2905 0 : bmap = isl_basic_map_deltas_map(bmap);
2906 :
2907 0 : return isl_union_map_from_map(isl_map_from_basic_map(bmap));
2908 : }
2909 :
2910 : /* Compute the positive powers of "map", or an overapproximation.
2911 : * The result maps the exponent to a nested copy of the corresponding power.
2912 : * If the result is exact, then *exact is set to 1.
2913 : */
2914 0 : __isl_give isl_union_map *isl_union_map_power(__isl_take isl_union_map *umap,
2915 : int *exact)
2916 : {
2917 : int n;
2918 : isl_union_map *inc;
2919 : isl_union_map *dm;
2920 :
2921 0 : if (!umap)
2922 0 : return NULL;
2923 0 : n = isl_union_map_n_map(umap);
2924 0 : if (n == 0)
2925 0 : return umap;
2926 0 : if (n == 1) {
2927 0 : struct isl_union_power up = { NULL, exact };
2928 0 : isl_union_map_foreach_map(umap, &power, &up);
2929 0 : isl_union_map_free(umap);
2930 0 : return up.pow;
2931 : }
2932 0 : inc = increment(isl_union_map_get_space(umap));
2933 0 : umap = isl_union_map_product(inc, umap);
2934 0 : umap = isl_union_map_transitive_closure(umap, exact);
2935 0 : umap = isl_union_map_zip(umap);
2936 0 : dm = deltas_map(isl_union_map_get_space(umap));
2937 0 : umap = isl_union_map_apply_domain(umap, dm);
2938 :
2939 0 : return umap;
2940 : }
2941 :
2942 : #undef TYPE
2943 : #define TYPE isl_map
2944 : #include "isl_power_templ.c"
2945 :
2946 : #undef TYPE
2947 : #define TYPE isl_union_map
2948 : #include "isl_power_templ.c"
|
