isl_transitive_closure.c

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