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