+ * map_i^+ \cup qc^+
+ *
+ * or
+ *
+ * map_i^+ \cup ((id \cup map_i^) \circ qc^+)
+ *
+ * or
+ *
+ * map_i^+ \cup (qc^+ \circ (id \cup map_i^))
+ *
+ * depending on whether left or right are NULL.
+ */
+static __isl_give isl_map *compute_incremental(
+ __isl_take isl_space *dim, __isl_keep isl_map *map,
+ int i, __isl_take isl_map *qc, int *left, int *right, int *exact)
+{
+ isl_map *map_i;
+ isl_map *tc;
+ isl_map *rtc = NULL;
+
+ if (!map)
+ goto error;
+ isl_assert(map->ctx, left || right, goto error);
+
+ map_i = isl_map_from_basic_map(isl_basic_map_copy(map->p[i]));
+ tc = construct_projected_component(isl_space_copy(dim), map_i,
+ exact, 1);
+ isl_map_free(map_i);
+
+ if (*exact)
+ qc = isl_map_transitive_closure(qc, exact);
+
+ if (!*exact) {
+ isl_space_free(dim);
+ isl_map_free(tc);
+ isl_map_free(qc);
+ return isl_map_universe(isl_map_get_space(map));
+ }
+
+ if (!left || !right)
+ rtc = isl_map_union(isl_map_copy(tc),
+ isl_map_identity(isl_map_get_space(tc)));
+ if (!right)
+ qc = isl_map_apply_range(rtc, qc);
+ if (!left)
+ qc = isl_map_apply_range(qc, rtc);
+ qc = isl_map_union(tc, qc);
+
+ isl_space_free(dim);
+
+ return qc;
+error:
+ isl_space_free(dim);
+ isl_map_free(qc);
+ return NULL;
+}
+
+/* Given a map "map", try to find a basic map such that
+ * map^+ can be computed as
+ *
+ * map^+ = map_i^+ \cup
+ * \bigcup_j ((map_i^+ \cup Id_C)^+ \circ map_j \circ (map_i^+ \cup Id_C))^+
+ *
+ * with C the simple hull of the domain and range of the input map.
+ * map_i^ \cup Id_C is computed by allowing the path lengths to be zero
+ * and by intersecting domain and range with C.
+ * Of course, we need to check that this is actually equal to map_i^ \cup Id_C.
+ * Also, we only use the incremental computation if all the transitive
+ * closures are exact and if the number of basic maps in the union,
+ * after computing the integer divisions, is smaller than the number
+ * of basic maps in the input map.
+ */
+static int incemental_on_entire_domain(__isl_keep isl_space *dim,
+ __isl_keep isl_map *map,
+ isl_set **dom, isl_set **ran, int *left, int *right,
+ __isl_give isl_map **res)
+{
+ int i;
+ isl_set *C;
+ unsigned d;
+
+ *res = NULL;
+
+ C = isl_set_union(isl_map_domain(isl_map_copy(map)),
+ isl_map_range(isl_map_copy(map)));
+ C = isl_set_from_basic_set(isl_set_simple_hull(C));
+ if (!C)
+ return -1;
+ if (C->n != 1) {
+ isl_set_free(C);
+ return 0;
+ }
+
+ d = isl_map_dim(map, isl_dim_in);
+
+ for (i = 0; i < map->n; ++i) {
+ isl_map *qc;
+ int exact_i, spurious;
+ int j;
+ dom[i] = isl_set_from_basic_set(isl_basic_map_domain(
+ isl_basic_map_copy(map->p[i])));
+ ran[i] = isl_set_from_basic_set(isl_basic_map_range(
+ isl_basic_map_copy(map->p[i])));
+ qc = q_closure(isl_space_copy(dim), isl_set_copy(C),
+ map->p[i], &exact_i);
+ if (!qc)
+ goto error;
+ if (!exact_i) {
+ isl_map_free(qc);
+ continue;
+ }
+ spurious = has_spurious_elements(qc, dom[i], ran[i]);
+ if (spurious) {
+ isl_map_free(qc);
+ if (spurious < 0)
+ goto error;
+ continue;
+ }
+ qc = isl_map_project_out(qc, isl_dim_in, d, 1);
+ qc = isl_map_project_out(qc, isl_dim_out, d, 1);
+ qc = isl_map_compute_divs(qc);
+ for (j = 0; j < map->n; ++j)
+ left[j] = right[j] = 1;
+ qc = compose(map, i, qc, left, right);
+ if (!qc)
+ goto error;
+ if (qc->n >= map->n) {
+ isl_map_free(qc);
+ continue;
+ }
+ *res = compute_incremental(isl_space_copy(dim), map, i, qc,
+ left, right, &exact_i);
+ if (!*res)
+ goto error;
+ if (exact_i)
+ break;
+ isl_map_free(*res);
+ *res = NULL;
+ }
+
+ isl_set_free(C);
+
+ return *res != NULL;
+error:
+ isl_set_free(C);
+ return -1;
+}
+
+/* Try and compute the transitive closure of "map" as
+ *
+ * map^+ = map_i^+ \cup
+ * \bigcup_j ((map_i^+ \cup Id_C)^+ \circ map_j \circ (map_i^+ \cup Id_C))^+
+ *
+ * with C either the simple hull of the domain and range of the entire
+ * map or the simple hull of domain and range of map_i.
+ */
+static __isl_give isl_map *incremental_closure(__isl_take isl_space *dim,
+ __isl_keep isl_map *map, int *exact, int project)
+{
+ int i;
+ isl_set **dom = NULL;
+ isl_set **ran = NULL;
+ int *left = NULL;
+ int *right = NULL;
+ isl_set *C;
+ unsigned d;
+ isl_map *res = NULL;
+
+ if (!project)
+ return construct_projected_component(dim, map, exact, project);
+
+ if (!map)
+ goto error;
+ if (map->n <= 1)
+ return construct_projected_component(dim, map, exact, project);
+
+ d = isl_map_dim(map, isl_dim_in);
+
+ dom = isl_calloc_array(map->ctx, isl_set *, map->n);
+ ran = isl_calloc_array(map->ctx, isl_set *, map->n);
+ left = isl_calloc_array(map->ctx, int, map->n);
+ right = isl_calloc_array(map->ctx, int, map->n);
+ if (!ran || !dom || !left || !right)
+ goto error;
+
+ if (incemental_on_entire_domain(dim, map, dom, ran, left, right, &res) < 0)
+ goto error;
+
+ for (i = 0; !res && i < map->n; ++i) {
+ isl_map *qc;
+ int exact_i, spurious, comp;
+ if (!dom[i])
+ dom[i] = isl_set_from_basic_set(
+ isl_basic_map_domain(
+ isl_basic_map_copy(map->p[i])));
+ if (!dom[i])
+ goto error;
+ if (!ran[i])
+ ran[i] = isl_set_from_basic_set(
+ isl_basic_map_range(
+ isl_basic_map_copy(map->p[i])));
+ if (!ran[i])
+ goto error;
+ C = isl_set_union(isl_set_copy(dom[i]),
+ isl_set_copy(ran[i]));
+ C = isl_set_from_basic_set(isl_set_simple_hull(C));
+ if (!C)
+ goto error;
+ if (C->n != 1) {
+ isl_set_free(C);
+ continue;
+ }
+ comp = composability(C, i, dom, ran, left, right, map);
+ if (!comp || comp < 0) {
+ isl_set_free(C);
+ if (comp < 0)
+ goto error;
+ continue;
+ }
+ qc = q_closure(isl_space_copy(dim), C, map->p[i], &exact_i);
+ if (!qc)
+ goto error;
+ if (!exact_i) {
+ isl_map_free(qc);
+ continue;
+ }
+ spurious = has_spurious_elements(qc, dom[i], ran[i]);
+ if (spurious) {
+ isl_map_free(qc);
+ if (spurious < 0)
+ goto error;
+ continue;
+ }
+ qc = isl_map_project_out(qc, isl_dim_in, d, 1);
+ qc = isl_map_project_out(qc, isl_dim_out, d, 1);
+ qc = isl_map_compute_divs(qc);
+ qc = compose(map, i, qc, (comp & LEFT) ? left : NULL,
+ (comp & RIGHT) ? right : NULL);
+ if (!qc)
+ goto error;
+ if (qc->n >= map->n) {
+ isl_map_free(qc);
+ continue;
+ }
+ res = compute_incremental(isl_space_copy(dim), map, i, qc,
+ (comp & LEFT) ? left : NULL,
+ (comp & RIGHT) ? right : NULL, &exact_i);
+ if (!res)
+ goto error;
+ if (exact_i)
+ break;
+ isl_map_free(res);
+ res = NULL;
+ }
+
+ for (i = 0; i < map->n; ++i) {
+ isl_set_free(dom[i]);
+ isl_set_free(ran[i]);
+ }
+ free(dom);
+ free(ran);
+ free(left);
+ free(right);
+
+ if (res) {
+ isl_space_free(dim);
+ return res;
+ }
+
+ return construct_projected_component(dim, map, exact, project);
+error:
+ if (dom)
+ for (i = 0; i < map->n; ++i)
+ isl_set_free(dom[i]);
+ free(dom);
+ if (ran)
+ for (i = 0; i < map->n; ++i)
+ isl_set_free(ran[i]);
+ free(ran);
+ free(left);
+ free(right);
+ isl_space_free(dim);
+ return NULL;
+}
+
+/* Given an array of sets "set", add "dom" at position "pos"
+ * and search for elements at earlier positions that overlap with "dom".
+ * If any can be found, then merge all of them, together with "dom", into
+ * a single set and assign the union to the first in the array,
+ * which becomes the new group leader for all groups involved in the merge.
+ * During the search, we only consider group leaders, i.e., those with
+ * group[i] = i, as the other sets have already been combined
+ * with one of the group leaders.
+ */
+static int merge(isl_set **set, int *group, __isl_take isl_set *dom, int pos)
+{
+ int i;
+
+ group[pos] = pos;
+ set[pos] = isl_set_copy(dom);
+
+ for (i = pos - 1; i >= 0; --i) {
+ int o;
+
+ if (group[i] != i)
+ continue;
+
+ o = isl_set_overlaps(set[i], dom);
+ if (o < 0)
+ goto error;
+ if (!o)
+ continue;
+
+ set[i] = isl_set_union(set[i], set[group[pos]]);
+ set[group[pos]] = NULL;
+ if (!set[i])
+ goto error;
+ group[group[pos]] = i;
+ group[pos] = i;
+ }
+
+ isl_set_free(dom);
+ return 0;
+error:
+ isl_set_free(dom);
+ return -1;
+}
+
+/* Replace each entry in the n by n grid of maps by the cross product
+ * with the relation { [i] -> [i + 1] }.
+ */
+static int add_length(__isl_keep isl_map *map, isl_map ***grid, int n)
+{
+ int i, j, k;
+ isl_space *dim;
+ isl_basic_map *bstep;
+ isl_map *step;
+ unsigned nparam;
+
+ if (!map)
+ return -1;
+
+ dim = isl_map_get_space(map);
+ nparam = isl_space_dim(dim, isl_dim_param);
+ dim = isl_space_drop_dims(dim, isl_dim_in, 0, isl_space_dim(dim, isl_dim_in));
+ dim = isl_space_drop_dims(dim, isl_dim_out, 0, isl_space_dim(dim, isl_dim_out));
+ dim = isl_space_add_dims(dim, isl_dim_in, 1);
+ dim = isl_space_add_dims(dim, isl_dim_out, 1);
+ bstep = isl_basic_map_alloc_space(dim, 0, 1, 0);
+ k = isl_basic_map_alloc_equality(bstep);
+ if (k < 0) {
+ isl_basic_map_free(bstep);
+ return -1;
+ }
+ isl_seq_clr(bstep->eq[k], 1 + isl_basic_map_total_dim(bstep));
+ isl_int_set_si(bstep->eq[k][0], 1);
+ isl_int_set_si(bstep->eq[k][1 + nparam], 1);
+ isl_int_set_si(bstep->eq[k][1 + nparam + 1], -1);
+ bstep = isl_basic_map_finalize(bstep);
+ step = isl_map_from_basic_map(bstep);
+
+ for (i = 0; i < n; ++i)
+ for (j = 0; j < n; ++j)
+ grid[i][j] = isl_map_product(grid[i][j],
+ isl_map_copy(step));
+
+ isl_map_free(step);
+
+ return 0;
+}
+
+/* The core of the Floyd-Warshall algorithm.
+ * Updates the given n x x matrix of relations in place.
+ *
+ * The algorithm iterates over all vertices. In each step, the whole
+ * matrix is updated to include all paths that go to the current vertex,
+ * possibly stay there a while (including passing through earlier vertices)
+ * and then come back. At the start of each iteration, the diagonal
+ * element corresponding to the current vertex is replaced by its
+ * transitive closure to account for all indirect paths that stay
+ * in the current vertex.
+ */
+static void floyd_warshall_iterate(isl_map ***grid, int n, int *exact)
+{
+ int r, p, q;
+
+ for (r = 0; r < n; ++r) {
+ int r_exact;
+ grid[r][r] = isl_map_transitive_closure(grid[r][r],
+ (exact && *exact) ? &r_exact : NULL);
+ if (exact && *exact && !r_exact)
+ *exact = 0;
+
+ for (p = 0; p < n; ++p)
+ for (q = 0; q < n; ++q) {
+ isl_map *loop;
+ if (p == r && q == r)
+ continue;
+ loop = isl_map_apply_range(
+ isl_map_copy(grid[p][r]),
+ isl_map_copy(grid[r][q]));
+ grid[p][q] = isl_map_union(grid[p][q], loop);
+ loop = isl_map_apply_range(
+ isl_map_copy(grid[p][r]),
+ isl_map_apply_range(
+ isl_map_copy(grid[r][r]),
+ isl_map_copy(grid[r][q])));
+ grid[p][q] = isl_map_union(grid[p][q], loop);
+ grid[p][q] = isl_map_coalesce(grid[p][q]);
+ }
+ }
+}
+
+/* Given a partition of the domains and ranges of the basic maps in "map",
+ * apply the Floyd-Warshall algorithm with the elements in the partition
+ * as vertices.
+ *
+ * In particular, there are "n" elements in the partition and "group" is
+ * an array of length 2 * map->n with entries in [0,n-1].
+ *
+ * We first construct a matrix of relations based on the partition information,
+ * apply Floyd-Warshall on this matrix of relations and then take the
+ * union of all entries in the matrix as the final result.
+ *
+ * If we are actually computing the power instead of the transitive closure,
+ * i.e., when "project" is not set, then the result should have the
+ * path lengths encoded as the difference between an extra pair of
+ * coordinates. We therefore apply the nested transitive closures
+ * to relations that include these lengths. In particular, we replace
+ * the input relation by the cross product with the unit length relation
+ * { [i] -> [i + 1] }.
+ */
+static __isl_give isl_map *floyd_warshall_with_groups(__isl_take isl_space *dim,
+ __isl_keep isl_map *map, int *exact, int project, int *group, int n)
+{
+ int i, j, k;
+ isl_map ***grid = NULL;
+ isl_map *app;
+
+ if (!map)
+ goto error;
+
+ if (n == 1) {
+ free(group);
+ return incremental_closure(dim, map, exact, project);
+ }
+
+ grid = isl_calloc_array(map->ctx, isl_map **, n);
+ if (!grid)
+ goto error;
+ for (i = 0; i < n; ++i) {
+ grid[i] = isl_calloc_array(map->ctx, isl_map *, n);
+ if (!grid[i])
+ goto error;
+ for (j = 0; j < n; ++j)
+ grid[i][j] = isl_map_empty(isl_map_get_space(map));
+ }
+
+ for (k = 0; k < map->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(map->p[k])));
+ }
+
+ if (!project && add_length(map, grid, n) < 0)
+ goto error;
+
+ floyd_warshall_iterate(grid, n, exact);
+
+ app = isl_map_empty(isl_map_get_space(map));
+
+ for (i = 0; i < n; ++i) {
+ for (j = 0; j < n; ++j)
+ app = isl_map_union(app, grid[i][j]);
+ free(grid[i]);
+ }
+ free(grid);
+
+ free(group);
+ isl_space_free(dim);
+
+ return app;
+error:
+ if (grid)
+ for (i = 0; i < n; ++i) {
+ if (!grid[i])
+ continue;
+ for (j = 0; j < n; ++j)
+ isl_map_free(grid[i][j]);
+ free(grid[i]);
+ }
+ free(grid);
+ free(group);
+ isl_space_free(dim);
+ return NULL;
+}
+
+/* Partition the domains and ranges of the n basic relations in list
+ * into disjoint cells.
+ *
+ * To find the partition, we simply consider all of the domains
+ * and ranges in turn and combine those that overlap.
+ * "set" contains the partition elements and "group" indicates
+ * to which partition element a given domain or range belongs.