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