+
+static int inc_count(__isl_take isl_map *map, void *user)
+{
+ int *n = user;
+
+ *n += map->n;
+
+ isl_map_free(map);
+
+ return 0;
+}
+
+static int collect_basic_map(__isl_take isl_map *map, void *user)
+{
+ int i;
+ isl_basic_map ***next = user;
+
+ for (i = 0; i < map->n; ++i) {
+ **next = isl_basic_map_copy(map->p[i]);
+ if (!**next)
+ goto error;
+ (*next)++;
+ }
+
+ isl_map_free(map);
+ return 0;
+error:
+ isl_map_free(map);
+ return -1;
+}
+
+/* Perform Floyd-Warshall on the given list of basic relations.
+ * The basic relations may live in different dimensions,
+ * but basic relations that get assigned to the diagonal of the
+ * grid have domains and ranges of the same dimension and so
+ * the standard algorithm can be used because the nested transitive
+ * closures are only applied to diagonal elements and because all
+ * compositions are peformed on relations with compatible domains and ranges.
+ */
+static __isl_give isl_union_map *union_floyd_warshall_on_list(isl_ctx *ctx,
+ __isl_keep isl_basic_map **list, int n, int *exact)
+{
+ int i, j, k;
+ int n_group;
+ int *group = NULL;
+ isl_set **set = NULL;
+ isl_map ***grid = NULL;
+ isl_union_map *app;
+
+ group = setup_groups(ctx, list, n, &set, &n_group);
+ if (!group)
+ goto error;
+
+ grid = isl_calloc_array(ctx, isl_map **, n_group);
+ if (!grid)
+ goto error;
+ for (i = 0; i < n_group; ++i) {
+ grid[i] = isl_calloc_array(ctx, isl_map *, n_group);
+ if (!grid[i])
+ goto error;
+ for (j = 0; j < n_group; ++j) {
+ isl_space *dim1, *dim2, *dim;
+ dim1 = isl_space_reverse(isl_set_get_space(set[i]));
+ dim2 = isl_set_get_space(set[j]);
+ dim = isl_space_join(dim1, dim2);
+ grid[i][j] = isl_map_empty(dim);
+ }
+ }
+
+ for (k = 0; k < n; ++k) {
+ i = group[2 * k];
+ j = group[2 * k + 1];
+ grid[i][j] = isl_map_union(grid[i][j],
+ isl_map_from_basic_map(
+ isl_basic_map_copy(list[k])));
+ }
+
+ floyd_warshall_iterate(grid, n_group, exact);
+
+ app = isl_union_map_empty(isl_map_get_space(grid[0][0]));
+
+ for (i = 0; i < n_group; ++i) {
+ for (j = 0; j < n_group; ++j)
+ app = isl_union_map_add_map(app, grid[i][j]);
+ free(grid[i]);
+ }
+ free(grid);
+
+ for (i = 0; i < 2 * n; ++i)
+ isl_set_free(set[i]);
+ free(set);
+
+ free(group);
+ return app;
+error:
+ if (grid)
+ for (i = 0; i < n_group; ++i) {
+ if (!grid[i])
+ continue;
+ for (j = 0; j < n_group; ++j)
+ isl_map_free(grid[i][j]);
+ free(grid[i]);
+ }
+ free(grid);
+ if (set) {
+ for (i = 0; i < 2 * n; ++i)
+ isl_set_free(set[i]);
+ free(set);
+ }
+ free(group);
+ return NULL;
+}
+
+/* Perform Floyd-Warshall on the given union relation.
+ * The implementation is very similar to that for non-unions.
+ * The main difference is that it is applied unconditionally.
+ * We first extract a list of basic maps from the union map
+ * and then perform the algorithm on this list.
+ */
+static __isl_give isl_union_map *union_floyd_warshall(
+ __isl_take isl_union_map *umap, int *exact)
+{
+ int i, n;
+ isl_ctx *ctx;
+ isl_basic_map **list = NULL;
+ isl_basic_map **next;
+ isl_union_map *res;
+
+ n = 0;
+ if (isl_union_map_foreach_map(umap, inc_count, &n) < 0)
+ goto error;
+
+ ctx = isl_union_map_get_ctx(umap);
+ list = isl_calloc_array(ctx, isl_basic_map *, n);
+ if (!list)
+ goto error;
+
+ next = list;
+ if (isl_union_map_foreach_map(umap, collect_basic_map, &next) < 0)
+ goto error;
+
+ res = union_floyd_warshall_on_list(ctx, list, n, exact);
+
+ if (list) {
+ for (i = 0; i < n; ++i)
+ isl_basic_map_free(list[i]);
+ free(list);
+ }
+
+ isl_union_map_free(umap);
+ return res;
+error:
+ if (list) {
+ for (i = 0; i < n; ++i)
+ isl_basic_map_free(list[i]);
+ free(list);
+ }
+ isl_union_map_free(umap);
+ return NULL;
+}
+
+/* Decompose the give union relation into strongly connected components.
+ * The implementation is essentially the same as that of
+ * construct_power_components with the major difference that all
+ * operations are performed on union maps.
+ */
+static __isl_give isl_union_map *union_components(
+ __isl_take isl_union_map *umap, int *exact)
+{
+ int i;
+ int n;
+ isl_ctx *ctx;
+ isl_basic_map **list = NULL;
+ isl_basic_map **next;
+ isl_union_map *path = NULL;
+ struct isl_tc_follows_data data;
+ struct isl_tarjan_graph *g = NULL;
+ int c, l;
+ int recheck = 0;
+
+ n = 0;
+ if (isl_union_map_foreach_map(umap, inc_count, &n) < 0)
+ goto error;
+
+ if (n <= 1)
+ return union_floyd_warshall(umap, exact);
+
+ ctx = isl_union_map_get_ctx(umap);
+ list = isl_calloc_array(ctx, isl_basic_map *, n);
+ if (!list)
+ goto error;
+
+ next = list;
+ if (isl_union_map_foreach_map(umap, collect_basic_map, &next) < 0)
+ goto error;
+
+ data.list = list;
+ data.check_closed = 0;
+ g = isl_tarjan_graph_init(ctx, n, &basic_map_follows, &data);
+ if (!g)
+ goto error;
+
+ c = 0;
+ i = 0;
+ l = n;
+ path = isl_union_map_empty(isl_union_map_get_space(umap));
+ while (l) {
+ isl_union_map *comp;
+ isl_union_map *path_comp, *path_comb;
+ comp = isl_union_map_empty(isl_union_map_get_space(umap));
+ while (g->order[i] != -1) {
+ comp = isl_union_map_add_map(comp,
+ isl_map_from_basic_map(
+ isl_basic_map_copy(list[g->order[i]])));
+ --l;
+ ++i;
+ }
+ path_comp = union_floyd_warshall(comp, exact);
+ path_comb = isl_union_map_apply_range(isl_union_map_copy(path),
+ isl_union_map_copy(path_comp));
+ path = isl_union_map_union(path, path_comp);
+ path = isl_union_map_union(path, path_comb);
+ ++i;
+ ++c;
+ }
+
+ if (c > 1 && data.check_closed && !*exact) {
+ int closed;
+
+ closed = isl_union_map_is_transitively_closed(path);
+ if (closed < 0)
+ goto error;
+ recheck = !closed;
+ }
+
+ isl_tarjan_graph_free(g);
+
+ for (i = 0; i < n; ++i)
+ isl_basic_map_free(list[i]);
+ free(list);
+
+ if (recheck) {
+ isl_union_map_free(path);
+ return union_floyd_warshall(umap, exact);
+ }
+
+ isl_union_map_free(umap);
+
+ return path;
+error:
+ isl_tarjan_graph_free(g);
+ if (list) {
+ for (i = 0; i < n; ++i)
+ isl_basic_map_free(list[i]);
+ free(list);
+ }
+ isl_union_map_free(umap);
+ isl_union_map_free(path);
+ return NULL;
+}
+
+/* Compute the transitive closure of "umap", or an overapproximation.
+ * If the result is exact, then *exact is set to 1.
+ */
+__isl_give isl_union_map *isl_union_map_transitive_closure(
+ __isl_take isl_union_map *umap, int *exact)
+{
+ int closed;
+
+ if (!umap)
+ return NULL;
+
+ if (exact)
+ *exact = 1;
+
+ umap = isl_union_map_compute_divs(umap);
+ umap = isl_union_map_coalesce(umap);
+ closed = isl_union_map_is_transitively_closed(umap);
+ if (closed < 0)
+ goto error;
+ if (closed)
+ return umap;
+ umap = union_components(umap, exact);
+ return umap;
+error:
+ isl_union_map_free(umap);
+ return NULL;
+}
+
+struct isl_union_power {
+ isl_union_map *pow;
+ int *exact;
+};
+
+static int power(__isl_take isl_map *map, void *user)
+{
+ struct isl_union_power *up = user;
+
+ map = isl_map_power(map, up->exact);
+ up->pow = isl_union_map_from_map(map);
+
+ return -1;
+}
+
+/* Construct a map [x] -> [x+1], with parameters prescribed by "dim".
+ */
+static __isl_give isl_union_map *increment(__isl_take isl_space *dim)
+{
+ int k;
+ isl_basic_map *bmap;
+
+ dim = isl_space_add_dims(dim, isl_dim_in, 1);
+ dim = isl_space_add_dims(dim, isl_dim_out, 1);
+ bmap = isl_basic_map_alloc_space(dim, 0, 1, 0);
+ k = isl_basic_map_alloc_equality(bmap);
+ if (k < 0)
+ goto error;
+ isl_seq_clr(bmap->eq[k], isl_basic_map_total_dim(bmap));
+ isl_int_set_si(bmap->eq[k][0], 1);
+ isl_int_set_si(bmap->eq[k][isl_basic_map_offset(bmap, isl_dim_in)], 1);
+ isl_int_set_si(bmap->eq[k][isl_basic_map_offset(bmap, isl_dim_out)], -1);
+ return isl_union_map_from_map(isl_map_from_basic_map(bmap));
+error:
+ isl_basic_map_free(bmap);
+ return NULL;
+}
+
+/* Construct a map [[x]->[y]] -> [y-x], with parameters prescribed by "dim".
+ */
+static __isl_give isl_union_map *deltas_map(__isl_take isl_space *dim)
+{
+ isl_basic_map *bmap;
+
+ dim = isl_space_add_dims(dim, isl_dim_in, 1);
+ dim = isl_space_add_dims(dim, isl_dim_out, 1);
+ bmap = isl_basic_map_universe(dim);
+ bmap = isl_basic_map_deltas_map(bmap);
+
+ return isl_union_map_from_map(isl_map_from_basic_map(bmap));
+}
+
+/* Compute the positive powers of "map", or an overapproximation.
+ * The result maps the exponent to a nested copy of the corresponding power.
+ * If the result is exact, then *exact is set to 1.
+ */
+__isl_give isl_union_map *isl_union_map_power(__isl_take isl_union_map *umap,
+ int *exact)
+{
+ int n;
+ isl_union_map *inc;
+ isl_union_map *dm;
+
+ if (!umap)
+ return NULL;
+ n = isl_union_map_n_map(umap);
+ if (n == 0)
+ return umap;
+ if (n == 1) {
+ struct isl_union_power up = { NULL, exact };
+ isl_union_map_foreach_map(umap, &power, &up);
+ isl_union_map_free(umap);
+ return up.pow;
+ }
+ inc = increment(isl_union_map_get_space(umap));
+ umap = isl_union_map_product(inc, umap);
+ umap = isl_union_map_transitive_closure(umap, exact);
+ umap = isl_union_map_zip(umap);
+ dm = deltas_map(isl_union_map_get_space(umap));
+ umap = isl_union_map_apply_domain(umap, dm);
+
+ return umap;
+}
+
+#undef TYPE
+#define TYPE isl_map
+#include "isl_power_templ.c"
+
+#undef TYPE
+#define TYPE isl_union_map
+#include "isl_power_templ.c"