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
  15 /* Given a map that represents a path with the length of the path
  16  * encoded as the difference between the last output coordindate
  17  * and the last input coordinate, set this length to either
  18  * exactly "length" (if "exactly" is set) or at least "length"
  19  * (if "exactly" is not set).
  20  */
  21 static __isl_give isl_map *set_path_length(__isl_take isl_map *map,
  22         int exactly, int length)
  23 {
  24         struct isl_dim *dim;
  25         struct isl_basic_map *bmap;
  26         unsigned d;
  27         unsigned nparam;
  28         int k;
  29         isl_int *c;
  30
  31         if (!map)
  32                 return NULL;
  33
  34         dim = isl_map_get_dim(map);
  35         d = isl_dim_size(dim, isl_dim_in);
  36         nparam = isl_dim_size(dim, isl_dim_param);
  37         bmap = isl_basic_map_alloc_dim(dim, 0, 1, 1);
  38         if (exactly) {
  39                 k = isl_basic_map_alloc_equality(bmap);
  40                 c = bmap->eq[k];
  41         } else {
  42                 k = isl_basic_map_alloc_inequality(bmap);
  43                 c = bmap->ineq[k];
  44         }
  45         if (k < 0)
  46                 goto error;
  47         isl_seq_clr(c, 1 + isl_basic_map_total_dim(bmap));
  48         isl_int_set_si(c[0], -length);
  49         isl_int_set_si(c[1 + nparam + d - 1], -1);
  50         isl_int_set_si(c[1 + nparam + d + d - 1], 1);
  51
  52         bmap = isl_basic_map_finalize(bmap);
  53         map = isl_map_intersect(map, isl_map_from_basic_map(bmap));
  54
  55         return map;
  56 error:
  57         isl_basic_map_free(bmap);
  58         isl_map_free(map);
  59         return NULL;
  60 }
  61
  62 /* Check whether the overapproximation of the power of "map" is exactly
  63  * the power of "map".  Let R be "map" and A_k the overapproximation.
  64  * The approximation is exact if
  65  *
  66  *      A_1 = R
  67  *      A_k = A_{k-1} \circ R                   k >= 2
  68  *
  69  * Since A_k is known to be an overapproximation, we only need to check
  70  *
  71  *      A_1 \subset R
  72  *      A_k \subset A_{k-1} \circ R             k >= 2
  73  *
  74  * In practice, "app" has an extra input and output coordinate
  75  * to encode the length of the path.  So, we first need to add
  76  * this coordinate to "map" and set the length of the path to
  77  * one.
  78  */
  79 static int check_power_exactness(__isl_take isl_map *map,
  80         __isl_take isl_map *app)
  81 {
  82         int exact;
  83         isl_map *app_1;
  84         isl_map *app_2;
  85
  86         map = isl_map_add(map, isl_dim_in, 1);
  87         map = isl_map_add(map, isl_dim_out, 1);
  88         map = set_path_length(map, 1, 1);
  89
  90         app_1 = set_path_length(isl_map_copy(app), 1, 1);
  91
  92         exact = isl_map_is_subset(app_1, map);
  93         isl_map_free(app_1);
  94
  95         if (!exact || exact < 0) {
  96                 isl_map_free(app);
  97                 isl_map_free(map);
  98                 return exact;
  99         }
 100
 101         app_1 = set_path_length(isl_map_copy(app), 0, 1);
 102         app_2 = set_path_length(app, 0, 2);
 103         app_1 = isl_map_apply_range(map, app_1);
 104
 105         exact = isl_map_is_subset(app_2, app_1);
 106
 107         isl_map_free(app_1);
 108         isl_map_free(app_2);
 109
 110         return exact;
 111 }
 112
 113 /* Check whether the overapproximation of the power of "map" is exactly
 114  * the power of "map", possibly after projecting out the power (if "project"
 115  * is set).
 116  *
 117  * If "project" is set and if "steps" can only result in acyclic paths,
 118  * then we check
 119  *
 120  *      A = R \cup (A \circ R)
 121  *
 122  * where A is the overapproximation with the power projected out, i.e.,
 123  * an overapproximation of the transitive closure.
 124  * More specifically, since A is known to be an overapproximation, we check
 125  *
 126  *      A \subset R \cup (A \circ R)
 127  *
 128  * Otherwise, we check if the power is exact.
 129  *
 130  * Note that "app" has an extra input and output coordinate to encode
 131  * the length of the part.  If we are only interested in the transitive
 132  * closure, then we can simply project out these coordinates first.
 133  */
 134 static int check_exactness(__isl_take isl_map *map, __isl_take isl_map *app,
 135         int project)
 136 {
 137         isl_map *test;
 138         int exact;
 139         unsigned d;
 140
 141         if (!project)
 142                 return check_power_exactness(map, app);
 143
 144         d = isl_map_dim(map, isl_dim_in);
 145         app = set_path_length(app, 0, 1);
 146         app = isl_map_project_out(app, isl_dim_in, d, 1);
 147         app = isl_map_project_out(app, isl_dim_out, d, 1);
 148
 149         test = isl_map_apply_range(isl_map_copy(map), isl_map_copy(app));
 150         test = isl_map_union(test, isl_map_copy(map));
 151
 152         exact = isl_map_is_subset(app, test);
 153
 154         isl_map_free(app);
 155         isl_map_free(test);
 156
 157         isl_map_free(map);
 158
 159         return exact;
 160 error:
 161         isl_map_free(app);
 162         isl_map_free(map);
 163         return -1;
 164 }
 165
 166 /*
 167  * The transitive closure implementation is based on the paper
 168  * "Computing the Transitive Closure of a Union of Affine Integer
 169  * Tuple Relations" by Anna Beletska, Denis Barthou, Wlodzimierz Bielecki and
 170  * Albert Cohen.
 171  */
 172
 173 /* Given a set of n offsets v_i (the rows of "steps"), construct a relation
 174  * of the given dimension specification (Z^{n+1} -> Z^{n+1})
 175  * that maps an element x to any element that can be reached
 176  * by taking a non-negative number of steps along any of
 177  * the extended offsets v'_i = [v_i 1].
 178  * That is, construct
 179  *
 180  * { [x] -> [y] : exists k_i >= 0, y = x + \sum_i k_i v'_i }
 181  *
 182  * For any element in this relation, the number of steps taken
 183  * is equal to the difference in the final coordinates.
 184  */
 185 static __isl_give isl_map *path_along_steps(__isl_take isl_dim *dim,
 186         __isl_keep isl_mat *steps)
 187 {
 188         int i, j, k;
 189         struct isl_basic_map *path = NULL;
 190         unsigned d;
 191         unsigned n;
 192         unsigned nparam;
 193
 194         if (!dim || !steps)
 195                 goto error;
 196
 197         d = isl_dim_size(dim, isl_dim_in);
 198         n = steps->n_row;
 199         nparam = isl_dim_size(dim, isl_dim_param);
 200
 201         path = isl_basic_map_alloc_dim(isl_dim_copy(dim), n, d, n);
 202
 203         for (i = 0; i < n; ++i) {
 204                 k = isl_basic_map_alloc_div(path);
 205                 if (k < 0)
 206                         goto error;
 207                 isl_assert(steps->ctx, i == k, goto error);
 208                 isl_int_set_si(path->div[k][0], 0);
 209         }
 210
 211         for (i = 0; i < d; ++i) {
 212                 k = isl_basic_map_alloc_equality(path);
 213                 if (k < 0)
 214                         goto error;
 215                 isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
 216                 isl_int_set_si(path->eq[k][1 + nparam + i], 1);
 217                 isl_int_set_si(path->eq[k][1 + nparam + d + i], -1);
 218                 if (i == d - 1)
 219                         for (j = 0; j < n; ++j)
 220                                 isl_int_set_si(path->eq[k][1 + nparam + 2 * d + j], 1);
 221                 else
 222                         for (j = 0; j < n; ++j)
 223                                 isl_int_set(path->eq[k][1 + nparam + 2 * d + j],
 224                                             steps->row[j][i]);
 225         }
 226
 227         for (i = 0; i < n; ++i) {
 228                 k = isl_basic_map_alloc_inequality(path);
 229                 if (k < 0)
 230                         goto error;
 231                 isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
 232                 isl_int_set_si(path->ineq[k][1 + nparam + 2 * d + i], 1);
 233         }
 234
 235         isl_dim_free(dim);
 236
 237         path = isl_basic_map_simplify(path);
 238         path = isl_basic_map_finalize(path);
 239         return isl_map_from_basic_map(path);
 240 error:
 241         isl_dim_free(dim);
 242         isl_basic_map_free(path);
 243         return NULL;
 244 }
 245
 246 #define IMPURE          0
 247 #define PURE_PARAM      1
 248 #define PURE_VAR        2
 249 #define MIXED           3
 250
 251 /* Return PURE_PARAM if only the coefficients of the parameters are non-zero.
 252  * Return PURE_VAR if only the coefficients of the set variables are non-zero.
 253  * Return MIXED if only the coefficients of the parameters and the set
 254  *      variables are non-zero and if moreover the parametric constant
 255  *      can never attain positive values.
 256  * Return IMPURE otherwise.
 257  */
 258 static int purity(__isl_keep isl_basic_set *bset, isl_int *c, int eq)
 259 {
 260         unsigned d;
 261         unsigned n_div;
 262         unsigned nparam;
 263         int k;
 264         int empty;
 265
 266         n_div = isl_basic_set_dim(bset, isl_dim_div);
 267         d = isl_basic_set_dim(bset, isl_dim_set);
 268         nparam = isl_basic_set_dim(bset, isl_dim_param);
 269
 270         if (isl_seq_first_non_zero(c + 1 + nparam + d, n_div) != -1)
 271                 return IMPURE;
 272         if (isl_seq_first_non_zero(c + 1, nparam) == -1)
 273                 return PURE_VAR;
 274         if (isl_seq_first_non_zero(c + 1 + nparam, d) == -1)
 275                 return PURE_PARAM;
 276         if (eq)
 277                 return IMPURE;
 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         isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
 288         empty = isl_basic_set_is_empty(bset);
 289         isl_basic_set_free(bset);
 290
 291         return empty < 0 ? -1 : empty ? MIXED : IMPURE;
 292 error:
 293         isl_basic_set_free(bset);
 294         return -1;
 295 }
 296
 297 /* Given a set of offsets "delta", construct a relation of the
 298  * given dimension specification (Z^{n+1} -> Z^{n+1}) that
 299  * is an overapproximation of the relations that
 300  * maps an element x to any element that can be reached
 301  * by taking a non-negative number of steps along any of
 302  * the elements in "delta".
 303  * That is, construct an approximation of
 304  *
 305  *      { [x] -> [y] : exists f \in \delta, k \in Z :
 306  *                                      y = x + k [f, 1] and k >= 0 }
 307  *
 308  * For any element in this relation, the number of steps taken
 309  * is equal to the difference in the final coordinates.
 310  *
 311  * In particular, let delta be defined as
 312  *
 313  *      \delta = [p] -> { [x] : A x + a >= and B p + b >= 0 and
 314  *                              C x + C'p + c >= 0 and
 315  *                              D x + D'p + d >= 0 }
 316  *
 317  * where the constraints C x + C'p + c >= 0 are such that the parametric
 318  * constant term of each constraint j, "C_j x + C'_j p + c_j",
 319  * can never attain positive values, then the relation is constructed as
 320  *
 321  *      { [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and
 322  *                      A f + k a >= 0 and B p + b >= 0 and
 323  *                      C f + C'p + c >= 0 and k >= 1 }
 324  *      union { [x] -> [x] }
 325  *
 326  * Existentially quantified variables in \delta are currently ignored.
 327  * This is safe, but leads to an additional overapproximation.
 328  */
 329 static __isl_give isl_map *path_along_delta(__isl_take isl_dim *dim,
 330         __isl_take isl_basic_set *delta)
 331 {
 332         isl_basic_map *path = NULL;
 333         unsigned d;
 334         unsigned n_div;
 335         unsigned nparam;
 336         unsigned off;
 337         int i, k;
 338
 339         if (!delta)
 340                 goto error;
 341         n_div = isl_basic_set_dim(delta, isl_dim_div);
 342         d = isl_basic_set_dim(delta, isl_dim_set);
 343         nparam = isl_basic_set_dim(delta, isl_dim_param);
 344         path = isl_basic_map_alloc_dim(isl_dim_copy(dim), n_div + d + 1,
 345                         d + 1 + delta->n_eq, delta->n_ineq + 1);
 346         off = 1 + nparam + 2 * (d + 1) + n_div;
 347
 348         for (i = 0; i < n_div + d + 1; ++i) {
 349                 k = isl_basic_map_alloc_div(path);
 350                 if (k < 0)
 351                         goto error;
 352                 isl_int_set_si(path->div[k][0], 0);
 353         }
 354
 355         for (i = 0; i < d + 1; ++i) {
 356                 k = isl_basic_map_alloc_equality(path);
 357                 if (k < 0)
 358                         goto error;
 359                 isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
 360                 isl_int_set_si(path->eq[k][1 + nparam + i], 1);
 361                 isl_int_set_si(path->eq[k][1 + nparam + d + 1 + i], -1);
 362                 isl_int_set_si(path->eq[k][off + i], 1);
 363         }
 364
 365         for (i = 0; i < delta->n_eq; ++i) {
 366                 int p = purity(delta, delta->eq[i], 1);
 367                 if (p < 0)
 368                         goto error;
 369                 if (p == IMPURE)
 370                         continue;
 371                 k = isl_basic_map_alloc_equality(path);
 372                 if (k < 0)
 373                         goto error;
 374                 isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
 375                 if (p == PURE_VAR) {
 376                         isl_seq_cpy(path->eq[k] + off,
 377                                     delta->eq[i] + 1 + nparam, d);
 378                         isl_int_set(path->eq[k][off + d], delta->eq[i][0]);
 379                 } else
 380                         isl_seq_cpy(path->eq[k], delta->eq[i], 1 + nparam);
 381         }
 382
 383         for (i = 0; i < delta->n_ineq; ++i) {
 384                 int p = purity(delta, delta->ineq[i], 0);
 385                 if (p < 0)
 386                         goto error;
 387                 if (p == IMPURE)
 388                         continue;
 389                 k = isl_basic_map_alloc_inequality(path);
 390                 if (k < 0)
 391                         goto error;
 392                 isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
 393                 if (p == PURE_VAR) {
 394                         isl_seq_cpy(path->ineq[k] + off,
 395                                     delta->ineq[i] + 1 + nparam, d);
 396                         isl_int_set(path->ineq[k][off + d], delta->ineq[i][0]);
 397                 } else if (p == PURE_PARAM) {
 398                         isl_seq_cpy(path->ineq[k], delta->ineq[i], 1 + nparam);
 399                 } else {
 400                         isl_seq_cpy(path->ineq[k] + off,
 401                                     delta->ineq[i] + 1 + nparam, d);
 402                         isl_seq_cpy(path->ineq[k], delta->ineq[i], 1 + nparam);
 403                 }
 404         }
 405
 406         k = isl_basic_map_alloc_inequality(path);
 407         if (k < 0)
 408                 goto error;
 409         isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
 410         isl_int_set_si(path->ineq[k][0], -1);
 411         isl_int_set_si(path->ineq[k][off + d], 1);
 412
 413         isl_basic_set_free(delta);
 414         path = isl_basic_map_finalize(path);
 415         return isl_basic_map_union(path,
 416                                 isl_basic_map_identity(isl_dim_domain(dim)));
 417 error:
 418         isl_dim_free(dim);
 419         isl_basic_set_free(delta);
 420         isl_basic_map_free(path);
 421         return NULL;
 422 }
 423
 424 /* Given a dimenion specification Z^{n+1} -> Z^{n+1} and a parameter "param",
 425  * construct a map that equates the parameter to the difference
 426  * in the final coordinates and imposes that this difference is positive.
 427  * That is, construct
 428  *
 429  *      { [x,x_s] -> [y,y_s] : k = y_s - x_s > 0 }
 430  */
 431 static __isl_give isl_map *equate_parameter_to_length(__isl_take isl_dim *dim,
 432         unsigned param)
 433 {
 434         struct isl_basic_map *bmap;
 435         unsigned d;
 436         unsigned nparam;
 437         int k;
 438
 439         d = isl_dim_size(dim, isl_dim_in);
 440         nparam = isl_dim_size(dim, isl_dim_param);
 441         bmap = isl_basic_map_alloc_dim(dim, 0, 1, 1);
 442         k = isl_basic_map_alloc_equality(bmap);
 443         if (k < 0)
 444                 goto error;
 445         isl_seq_clr(bmap->eq[k], 1 + isl_basic_map_total_dim(bmap));
 446         isl_int_set_si(bmap->eq[k][1 + param], -1);
 447         isl_int_set_si(bmap->eq[k][1 + nparam + d - 1], -1);
 448         isl_int_set_si(bmap->eq[k][1 + nparam + d + d - 1], 1);
 449
 450         k = isl_basic_map_alloc_inequality(bmap);
 451         if (k < 0)
 452                 goto error;
 453         isl_seq_clr(bmap->ineq[k], 1 + isl_basic_map_total_dim(bmap));
 454         isl_int_set_si(bmap->ineq[k][1 + param], 1);
 455         isl_int_set_si(bmap->ineq[k][0], -1);
 456
 457         bmap = isl_basic_map_finalize(bmap);
 458         return isl_map_from_basic_map(bmap);
 459 error:
 460         isl_basic_map_free(bmap);
 461         return NULL;
 462 }
 463
 464 /* Check whether "path" is acyclic, where the last coordinates of domain
 465  * and range of path encode the number of steps taken.
 466  * That is, check whether
 467  *
 468  *      { d | d = y - x and (x,y) in path }
 469  *
 470  * does not contain any element with positive last coordinate (positive length)
 471  * and zero remaining coordinates (cycle).
 472  */
 473 static int is_acyclic(__isl_take isl_map *path)
 474 {
 475         int i;
 476         int acyclic;
 477         unsigned dim;
 478         struct isl_set *delta;
 479
 480         delta = isl_map_deltas(path);
 481         dim = isl_set_dim(delta, isl_dim_set);
 482         for (i = 0; i < dim; ++i) {
 483                 if (i == dim -1)
 484                         delta = isl_set_lower_bound_si(delta, isl_dim_set, i, 1);
 485                 else
 486                         delta = isl_set_fix_si(delta, isl_dim_set, i, 0);
 487         }
 488
 489         acyclic = isl_set_is_empty(delta);
 490         isl_set_free(delta);
 491
 492         return acyclic;
 493 }
 494
 495 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
 496  * and a dimension specification (Z^{n+1} -> Z^{n+1}),
 497  * construct a map that is an overapproximation of the map
 498  * that takes an element from the space D \times Z to another
 499  * element from the same space, such that the first n coordinates of the
 500  * difference between them is a sum of differences between images
 501  * and pre-images in one of the R_i and such that the last coordinate
 502  * is equal to the number of steps taken.
 503  * That is, let
 504  *
 505  *      \Delta_i = { y - x | (x, y) in R_i }
 506  *
 507  * then the constructed map is an overapproximation of
 508  *
 509  *      { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
 510  *                              d = (\sum_i k_i \delta_i, \sum_i k_i) }
 511  *
 512  * The elements of the singleton \Delta_i's are collected as the
 513  * rows of the steps matrix.  For all these \Delta_i's together,
 514  * a single path is constructed.
 515  * For each of the other \Delta_i's, we compute an overapproximation
 516  * of the paths along elements of \Delta_i.
 517  * Since each of these paths performs an addition, composition is
 518  * symmetric and we can simply compose all resulting paths in any order.
 519  */
 520 static __isl_give isl_map *construct_extended_path(__isl_take isl_dim *dim,
 521         __isl_keep isl_map *map, int *project)
 522 {
 523         struct isl_mat *steps = NULL;
 524         struct isl_map *path = NULL;
 525         unsigned d;
 526         int i, j, n;
 527
 528         d = isl_map_dim(map, isl_dim_in);
 529
 530         path = isl_map_identity(isl_dim_domain(isl_dim_copy(dim)));
 531
 532         steps = isl_mat_alloc(map->ctx, map->n, d);
 533         if (!steps)
 534                 goto error;
 535
 536         n = 0;
 537         for (i = 0; i < map->n; ++i) {
 538                 struct isl_basic_set *delta;
 539
 540                 delta = isl_basic_map_deltas(isl_basic_map_copy(map->p[i]));
 541
 542                 for (j = 0; j < d; ++j) {
 543                         int fixed;
 544
 545                         fixed = isl_basic_set_fast_dim_is_fixed(delta, j,
 546                                                             &steps->row[n][j]);
 547                         if (fixed < 0) {
 548                                 isl_basic_set_free(delta);
 549                                 goto error;
 550                         }
 551                         if (!fixed)
 552                                 break;
 553                 }
 554
 555
 556                 if (j < d) {
 557                         path = isl_map_apply_range(path,
 558                                 path_along_delta(isl_dim_copy(dim), delta));
 559                         path = isl_map_coalesce(path);
 560                 } else {
 561                         isl_basic_set_free(delta);
 562                         ++n;
 563                 }
 564         }
 565
 566         if (n > 0) {
 567                 steps->n_row = n;
 568                 path = isl_map_apply_range(path,
 569                                 path_along_steps(isl_dim_copy(dim), steps));
 570         }
 571
 572         if (project && *project) {
 573                 *project = is_acyclic(isl_map_copy(path));
 574                 if (*project < 0)
 575                         goto error;
 576         }
 577
 578         isl_dim_free(dim);
 579         isl_mat_free(steps);
 580         return path;
 581 error:
 582         isl_dim_free(dim);
 583         isl_mat_free(steps);
 584         isl_map_free(path);
 585         return NULL;
 586 }
 587
 588 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
 589  * and a dimension specification (Z^{n+1} -> Z^{n+1}),
 590  * construct a map that is the union of the identity map and
 591  * an overapproximation of the map
 592  * that takes an element from the dom R \times Z to an
 593  * element from ran R \times Z, such that the first n coordinates of the
 594  * difference between them is a sum of differences between images
 595  * and pre-images in one of the R_i and such that the last coordinate
 596  * is equal to the number of steps taken.
 597  * That is, let
 598  *
 599  *      \Delta_i = { y - x | (x, y) in R_i }
 600  *
 601  * then the constructed map is an overapproximation of
 602  *
 603  *      { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
 604  *                              d = (\sum_i k_i \delta_i, \sum_i k_i) and
 605  *                              x in dom R and x + d in ran R } union
 606  *      { (x) -> (x) }
 607  */
 608 static __isl_give isl_map *construct_component(__isl_take isl_dim *dim,
 609         __isl_keep isl_map *map, int *exact, int project)
 610 {
 611         struct isl_set *domain = NULL;
 612         struct isl_set *range = NULL;
 613         struct isl_set *overlap;
 614         struct isl_map *app = NULL;
 615         struct isl_map *path = NULL;
 616
 617         domain = isl_map_domain(isl_map_copy(map));
 618         domain = isl_set_coalesce(domain);
 619         range = isl_map_range(isl_map_copy(map));
 620         range = isl_set_coalesce(range);
 621         overlap = isl_set_intersect(isl_set_copy(domain), isl_set_copy(range));
 622         if (isl_set_is_empty(overlap) == 1) {
 623                 isl_set_free(domain);
 624                 isl_set_free(range);
 625                 isl_set_free(overlap);
 626                 isl_dim_free(dim);
 627
 628                 map = isl_map_copy(map);
 629                 map = isl_map_add(map, isl_dim_in, 1);
 630                 map = isl_map_add(map, isl_dim_out, 1);
 631                 map = set_path_length(map, 1, 1);
 632                 return map;
 633         }
 634         isl_set_free(overlap);
 635         app = isl_map_from_domain_and_range(domain, range);
 636         app = isl_map_add(app, isl_dim_in, 1);
 637         app = isl_map_add(app, isl_dim_out, 1);
 638
 639         path = construct_extended_path(isl_dim_copy(dim), map,
 640                                         exact && *exact ? &project : NULL);
 641         app = isl_map_intersect(app, path);
 642
 643         if (exact && *exact &&
 644             (*exact = check_exactness(isl_map_copy(map), isl_map_copy(app),
 645                                       project)) < 0)
 646                 goto error;
 647
 648         return isl_map_union(app, isl_map_identity(isl_dim_domain(dim)));
 649 error:
 650         isl_dim_free(dim);
 651         isl_map_free(app);
 652         return NULL;
 653 }
 654
 655 /* Structure for representing the nodes in the graph being traversed
 656  * using Tarjan's algorithm.
 657  * index represents the order in which nodes are visited.
 658  * min_index is the index of the root of a (sub)component.
 659  * on_stack indicates whether the node is currently on the stack.
 660  */
 661 struct basic_map_sort_node {
 662         int index;
 663         int min_index;
 664         int on_stack;
 665 };
 666 /* Structure for representing the graph being traversed
 667  * using Tarjan's algorithm.
 668  * len is the number of nodes
 669  * node is an array of nodes
 670  * stack contains the nodes on the path from the root to the current node
 671  * sp is the stack pointer
 672  * index is the index of the last node visited
 673  * order contains the elements of the components separated by -1
 674  * op represents the current position in order
 675  */
 676 struct basic_map_sort {
 677         int len;
 678         struct basic_map_sort_node *node;
 679         int *stack;
 680         int sp;
 681         int index;
 682         int *order;
 683         int op;
 684 };
 685
 686 static void basic_map_sort_free(struct basic_map_sort *s)
 687 {
 688         if (!s)
 689                 return;
 690         free(s->node);
 691         free(s->stack);
 692         free(s->order);
 693         free(s);
 694 }
 695
 696 static struct basic_map_sort *basic_map_sort_alloc(struct isl_ctx *ctx, int len)
 697 {
 698         struct basic_map_sort *s;
 699         int i;
 700
 701         s = isl_calloc_type(ctx, struct basic_map_sort);
 702         if (!s)
 703                 return NULL;
 704         s->len = len;
 705         s->node = isl_alloc_array(ctx, struct basic_map_sort_node, len);
 706         if (!s->node)
 707                 goto error;
 708         for (i = 0; i < len; ++i)
 709                 s->node[i].index = -1;
 710         s->stack = isl_alloc_array(ctx, int, len);
 711         if (!s->stack)
 712                 goto error;
 713         s->order = isl_alloc_array(ctx, int, 2 * len);
 714         if (!s->order)
 715                 goto error;
 716
 717         s->sp = 0;
 718         s->index = 0;
 719         s->op = 0;
 720
 721         return s;
 722 error:
 723         basic_map_sort_free(s);
 724         return NULL;
 725 }
 726
 727 /* Check whether in the computation of the transitive closure
 728  * "bmap1" (R_1) should follow (or be part of the same component as)
 729  * "bmap2" (R_2).
 730  *
 731  * That is check whether
 732  *
 733  *      R_1 \circ R_2
 734  *
 735  * is a subset of
 736  *
 737  *      R_2 \circ R_1
 738  *
 739  * If so, then there is no reason for R_1 to immediately follow R_2
 740  * in any path.
 741  */
 742 static int basic_map_follows(__isl_keep isl_basic_map *bmap1,
 743         __isl_keep isl_basic_map *bmap2)
 744 {
 745         struct isl_map *map12 = NULL;
 746         struct isl_map *map21 = NULL;
 747         int subset;
 748
 749         map21 = isl_map_from_basic_map(
 750                         isl_basic_map_apply_range(
 751                                 isl_basic_map_copy(bmap2),
 752                                 isl_basic_map_copy(bmap1)));
 753         subset = isl_map_is_empty(map21);
 754         if (subset < 0)
 755                 goto error;
 756         if (subset) {
 757                 isl_map_free(map21);
 758                 return 0;
 759         }
 760
 761         map12 = isl_map_from_basic_map(
 762                         isl_basic_map_apply_range(
 763                                 isl_basic_map_copy(bmap1),
 764                                 isl_basic_map_copy(bmap2)));
 765
 766         subset = isl_map_is_subset(map21, map12);
 767
 768         isl_map_free(map12);
 769         isl_map_free(map21);
 770
 771         return subset < 0 ? -1 : !subset;
 772 error:
 773         isl_map_free(map21);
 774         return -1;
 775 }
 776
 777 /* Perform Tarjan's algorithm for computing the strongly connected components
 778  * in the graph with the disjuncts of "map" as vertices and with an
 779  * edge between any pair of disjuncts such that the first has
 780  * to be applied after the second.
 781  */
 782 static int power_components_tarjan(struct basic_map_sort *s,
 783         __isl_keep isl_map *map, int i)
 784 {
 785         int j;
 786
 787         s->node[i].index = s->index;
 788         s->node[i].min_index = s->index;
 789         s->node[i].on_stack = 1;
 790         s->index++;
 791         s->stack[s->sp++] = i;
 792
 793         for (j = s->len - 1; j >= 0; --j) {
 794                 int f;
 795
 796                 if (j == i)
 797                         continue;
 798                 if (s->node[j].index >= 0 &&
 799                         (!s->node[j].on_stack ||
 800                          s->node[j].index > s->node[i].min_index))
 801                         continue;
 802
 803                 f = basic_map_follows(map->p[i], map->p[j]);
 804                 if (f < 0)
 805                         return -1;
 806                 if (!f)
 807                         continue;
 808
 809                 if (s->node[j].index < 0) {
 810                         power_components_tarjan(s, map, j);
 811                         if (s->node[j].min_index < s->node[i].min_index)
 812                                 s->node[i].min_index = s->node[j].min_index;
 813                 } else if (s->node[j].index < s->node[i].min_index)
 814                         s->node[i].min_index = s->node[j].index;
 815         }
 816
 817         if (s->node[i].index != s->node[i].min_index)
 818                 return 0;
 819
 820         do {
 821                 j = s->stack[--s->sp];
 822                 s->node[j].on_stack = 0;
 823                 s->order[s->op++] = j;
 824         } while (j != i);
 825         s->order[s->op++] = -1;
 826
 827         return 0;
 828 }
 829
 830 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
 831  * and a dimension specification (Z^{n+1} -> Z^{n+1}),
 832  * construct a map that is the union of the identity map and
 833  * an overapproximation of the map
 834  * that takes an element from the dom R \times Z to an
 835  * element from ran R \times Z, such that the first n coordinates of the
 836  * difference between them is a sum of differences between images
 837  * and pre-images in one of the R_i and such that the last coordinate
 838  * is equal to the number of steps taken.
 839  * That is, let
 840  *
 841  *      \Delta_i = { y - x | (x, y) in R_i }
 842  *
 843  * then the constructed map is an overapproximation of
 844  *
 845  *      { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
 846  *                              d = (\sum_i k_i \delta_i, \sum_i k_i) and
 847  *                              x in dom R and x + d in ran R } union
 848  *      { (x) -> (x) }
 849  *
 850  * We first split the map into strongly connected components, perform
 851  * the above on each component and the join the results in the correct
 852  * order.  The power of each of the components needs to be extended
 853  * with the identity map because a path in the global result need
 854  * not go through every component.
 855  * The final result will then also contain the identity map, but
 856  * this part will be removed when the length of the path is forced
 857  * to be strictly positive.
 858  */
 859 static __isl_give isl_map *construct_power_components(__isl_take isl_dim *dim,
 860         __isl_keep isl_map *map, int *exact, int project)
 861 {
 862         int i, n;
 863         struct isl_map *path = NULL;
 864         struct basic_map_sort *s = NULL;
 865
 866         if (!map)
 867                 goto error;
 868         if (map->n <= 1)
 869                 return construct_component(dim, map, exact, project);
 870
 871         s = basic_map_sort_alloc(map->ctx, map->n);
 872         if (!s)
 873                 goto error;
 874         for (i = map->n - 1; i >= 0; --i) {
 875                 if (s->node[i].index >= 0)
 876                         continue;
 877                 if (power_components_tarjan(s, map, i) < 0)
 878                         goto error;
 879         }
 880
 881         i = 0;
 882         n = map->n;
 883         path = isl_map_identity(isl_dim_domain(isl_dim_copy(dim)));
 884         while (n) {
 885                 struct isl_map *comp;
 886                 comp = isl_map_alloc_dim(isl_map_get_dim(map), n, 0);
 887                 while (s->order[i] != -1) {
 888                         comp = isl_map_add_basic_map(comp,
 889                                     isl_basic_map_copy(map->p[s->order[i]]));
 890                         --n;
 891                         ++i;
 892                 }
 893                 path = isl_map_apply_range(path,
 894                             construct_component(isl_dim_copy(dim), comp,
 895                                                 exact, project));
 896                 isl_map_free(comp);
 897                 ++i;
 898         }
 899
 900         basic_map_sort_free(s);
 901         isl_dim_free(dim);
 902
 903         return path;
 904 error:
 905         basic_map_sort_free(s);
 906         isl_dim_free(dim);
 907         return NULL;
 908 }
 909
 910 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D,
 911  * construct a map that is an overapproximation of the map
 912  * that takes an element from the space D to another
 913  * element from the same space, such that the difference between
 914  * them is a strictly positive sum of differences between images
 915  * and pre-images in one of the R_i.
 916  * The number of differences in the sum is equated to parameter "param".
 917  * That is, let
 918  *
 919  *      \Delta_i = { y - x | (x, y) in R_i }
 920  *
 921  * then the constructed map is an overapproximation of
 922  *
 923  *      { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
 924  *                              d = \sum_i k_i \delta_i and k = \sum_i k_i > 0 }
 925  *
 926  * We first construct an extended mapping with an extra coordinate
 927  * that indicates the number of steps taken.  In particular,
 928  * the difference in the last coordinate is equal to the number
 929  * of steps taken to move from a domain element to the corresponding
 930  * image element(s).
 931  * In the final step, this difference is equated to the parameter "param"
 932  * and made positive.  The extra coordinates are subsequently projected out.
 933  */
 934 static __isl_give isl_map *construct_power(__isl_keep isl_map *map,
 935         unsigned param, int *exact, int project)
 936 {
 937         struct isl_map *app = NULL;
 938         struct isl_map *diff;
 939         struct isl_dim *dim = NULL;
 940         unsigned d;
 941
 942         if (!map)
 943                 return NULL;
 944
 945         dim = isl_map_get_dim(map);
 946
 947         d = isl_dim_size(dim, isl_dim_in);
 948         dim = isl_dim_add(dim, isl_dim_in, 1);
 949         dim = isl_dim_add(dim, isl_dim_out, 1);
 950
 951         app = construct_power_components(isl_dim_copy(dim), map,
 952                                         exact, project);
 953
 954         diff = equate_parameter_to_length(dim, param);
 955         app = isl_map_intersect(app, diff);
 956         app = isl_map_project_out(app, isl_dim_in, d, 1);
 957         app = isl_map_project_out(app, isl_dim_out, d, 1);
 958
 959         return app;
 960 }
 961
 962 /* Compute the positive powers of "map", or an overapproximation.
 963  * The power is given by parameter "param".  If the result is exact,
 964  * then *exact is set to 1.
 965  * If project is set, then we are actually interested in the transitive
 966  * closure, so we can use a more relaxed exactness check.
 967  */
 968 static __isl_give isl_map *map_power(__isl_take isl_map *map, unsigned param,
 969         int *exact, int project)
 970 {
 971         struct isl_map *app = NULL;
 972
 973         if (exact)
 974                 *exact = 1;
 975
 976         map = isl_map_remove_empty_parts(map);
 977         if (!map)
 978                 return NULL;
 979
 980         if (isl_map_fast_is_empty(map))
 981                 return map;
 982
 983         isl_assert(map->ctx, param < isl_map_dim(map, isl_dim_param), goto error);
 984         isl_assert(map->ctx,
 985                 isl_map_dim(map, isl_dim_in) == isl_map_dim(map, isl_dim_out),
 986                 goto error);
 987
 988         app = construct_power(map, param, exact, project);
 989
 990         isl_map_free(map);
 991         return app;
 992 error:
 993         isl_map_free(map);
 994         isl_map_free(app);
 995         return NULL;
 996 }
 997
 998 /* Compute the positive powers of "map", or an overapproximation.
 999  * The power is given by parameter "param".  If the result is exact,
1000  * then *exact is set to 1.
1001  */
1002 __isl_give isl_map *isl_map_power(__isl_take isl_map *map, unsigned param,
1003         int *exact)
1004 {
1005         return map_power(map, param, exact, 0);
1006 }
1007
1008 /* Compute the transitive closure  of "map", or an overapproximation.
1009  * If the result is exact, then *exact is set to 1.
1010  * Simply compute the powers of map and then project out the parameter
1011  * describing the power.
1012  */
1013 __isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
1014         int *exact)
1015 {
1016         unsigned param;
1017
1018         if (!map)
1019                 goto error;
1020
1021         param = isl_map_dim(map, isl_dim_param);
1022         map = isl_map_add(map, isl_dim_param, 1);
1023         map = map_power(map, param, exact, 1);
1024         map = isl_map_project_out(map, isl_dim_param, param, 1);
1025
1026         return map;
1027 error:
1028         isl_map_free(map);
1029         return NULL;
1030 }