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                 } else {
 560                         isl_basic_set_free(delta);
 561                         ++n;
 562                 }
 563         }
 564
 565         if (n > 0) {
 566                 steps->n_row = n;
 567                 path = isl_map_apply_range(path,
 568                                 path_along_steps(isl_dim_copy(dim), steps));
 569         }
 570
 571         if (project && *project) {
 572                 *project = is_acyclic(isl_map_copy(path));
 573                 if (*project < 0)
 574                         goto error;
 575         }
 576
 577         isl_dim_free(dim);
 578         isl_mat_free(steps);
 579         return path;
 580 error:
 581         isl_dim_free(dim);
 582         isl_mat_free(steps);
 583         isl_map_free(path);
 584         return NULL;
 585 }
 586
 587 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
 588  * and a dimension specification (Z^{n+1} -> Z^{n+1}),
 589  * construct a map that is the union of the identity map and
 590  * an overapproximation of the map
 591  * that takes an element from the dom R \times Z to an
 592  * element from ran R \times Z, such that the first n coordinates of the
 593  * difference between them is a sum of differences between images
 594  * and pre-images in one of the R_i and such that the last coordinate
 595  * is equal to the number of steps taken.
 596  * That is, let
 597  *
 598  *      \Delta_i = { y - x | (x, y) in R_i }
 599  *
 600  * then the constructed map is an overapproximation of
 601  *
 602  *      { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
 603  *                              d = (\sum_i k_i \delta_i, \sum_i k_i) and
 604  *                              x in dom R and x + d in ran R } union
 605  *      { (x) -> (x) }
 606  */
 607 static __isl_give isl_map *construct_component(__isl_take isl_dim *dim,
 608         __isl_keep isl_map *map, int *exact, int project)
 609 {
 610         struct isl_set *domain = NULL;
 611         struct isl_set *range = NULL;
 612         struct isl_map *app = NULL;
 613         struct isl_map *path = NULL;
 614
 615         domain = isl_map_domain(isl_map_copy(map));
 616         domain = isl_set_coalesce(domain);
 617         range = isl_map_range(isl_map_copy(map));
 618         range = isl_set_coalesce(range);
 619         app = isl_map_from_domain_and_range(domain, range);
 620         app = isl_map_add(app, isl_dim_in, 1);
 621         app = isl_map_add(app, isl_dim_out, 1);
 622
 623         path = construct_extended_path(isl_dim_copy(dim), map,
 624                                         exact && *exact ? &project : NULL);
 625         app = isl_map_intersect(app, path);
 626
 627         if (exact && *exact &&
 628             (*exact = check_exactness(isl_map_copy(map), isl_map_copy(app),
 629                                       project)) < 0)
 630                 goto error;
 631
 632         return isl_map_union(app, isl_map_identity(isl_dim_domain(dim)));
 633 error:
 634         isl_dim_free(dim);
 635         isl_map_free(app);
 636         return NULL;
 637 }
 638
 639 /* Structure for representing the nodes in the graph being traversed
 640  * using Tarjan's algorithm.
 641  * index represents the order in which nodes are visited.
 642  * min_index is the index of the root of a (sub)component.
 643  * on_stack indicates whether the node is currently on the stack.
 644  */
 645 struct basic_map_sort_node {
 646         int index;
 647         int min_index;
 648         int on_stack;
 649 };
 650 /* Structure for representing the graph being traversed
 651  * using Tarjan's algorithm.
 652  * len is the number of nodes
 653  * node is an array of nodes
 654  * stack contains the nodes on the path from the root to the current node
 655  * sp is the stack pointer
 656  * index is the index of the last node visited
 657  * order contains the elements of the components separated by -1
 658  * op represents the current position in order
 659  */
 660 struct basic_map_sort {
 661         int len;
 662         struct basic_map_sort_node *node;
 663         int *stack;
 664         int sp;
 665         int index;
 666         int *order;
 667         int op;
 668 };
 669
 670 static void basic_map_sort_free(struct basic_map_sort *s)
 671 {
 672         if (!s)
 673                 return;
 674         free(s->node);
 675         free(s->stack);
 676         free(s->order);
 677         free(s);
 678 }
 679
 680 static struct basic_map_sort *basic_map_sort_alloc(struct isl_ctx *ctx, int len)
 681 {
 682         struct basic_map_sort *s;
 683         int i;
 684
 685         s = isl_calloc_type(ctx, struct basic_map_sort);
 686         if (!s)
 687                 return NULL;
 688         s->len = len;
 689         s->node = isl_alloc_array(ctx, struct basic_map_sort_node, len);
 690         if (!s->node)
 691                 goto error;
 692         for (i = 0; i < len; ++i)
 693                 s->node[i].index = -1;
 694         s->stack = isl_alloc_array(ctx, int, len);
 695         if (!s->stack)
 696                 goto error;
 697         s->order = isl_alloc_array(ctx, int, 2 * len);
 698         if (!s->order)
 699                 goto error;
 700
 701         s->sp = 0;
 702         s->index = 0;
 703         s->op = 0;
 704
 705         return s;
 706 error:
 707         basic_map_sort_free(s);
 708         return NULL;
 709 }
 710
 711 /* Check whether in the computation of the transitive closure
 712  * "bmap1" (R_1) should follow (or be part of the same component as)
 713  * "bmap2" (R_2).
 714  *
 715  * That is check whether
 716  *
 717  *      R_1 \circ R_2
 718  *
 719  * is non-empty and that moreover, it is non-empty on the set
 720  * of elements that do not get mapped to the same set of elements
 721  * by both "R_1 \circ R_2" and "R_2 \circ R_1".
 722  * For elements that do get mapped to the same elements by these
 723  * two compositions, R_1 and R_2 are commutative, so if these
 724  * elements are the only ones for which R_1 \circ R_2 is non-empty,
 725  * then you may just as well apply R_1 first.
 726  */
 727 static int basic_map_follows(__isl_keep isl_basic_map *bmap1,
 728         __isl_keep isl_basic_map *bmap2)
 729 {
 730         struct isl_map *map12 = NULL;
 731         struct isl_map *map21 = NULL;
 732         struct isl_map *d = NULL;
 733         struct isl_set *dom = NULL;
 734         int empty;
 735
 736         map21 = isl_map_from_basic_map(
 737                         isl_basic_map_apply_range(
 738                                 isl_basic_map_copy(bmap2),
 739                                 isl_basic_map_copy(bmap1)));
 740         empty = isl_map_is_empty(map21);
 741         if (empty < 0)
 742                 goto error;
 743         if (empty) {
 744                 isl_map_free(map21);
 745                 return 0;
 746         }
 747
 748         map12 = isl_map_from_basic_map(
 749                         isl_basic_map_apply_range(
 750                                 isl_basic_map_copy(bmap1),
 751                                 isl_basic_map_copy(bmap2)));
 752         d = isl_map_subtract(isl_map_copy(map12), isl_map_copy(map21));
 753         d = isl_map_union(d,
 754                 isl_map_subtract(isl_map_copy(map21), isl_map_copy(map12)));
 755         dom = isl_map_domain(d);
 756
 757         map21 = isl_map_intersect_domain(map21, dom);
 758         empty = isl_map_is_empty(map21);
 759
 760         isl_map_free(map12);
 761         isl_map_free(map21);
 762
 763         return empty < 0 ? -1 : !empty;
 764 error:
 765         isl_map_free(map21);
 766         return -1;
 767 }
 768
 769 /* Perform Tarjan's algorithm for computing the strongly connected components
 770  * in the graph with the disjuncts of "map" as vertices and with an
 771  * edge between any pair of disjuncts such that the first has
 772  * to be applied after the second.
 773  */
 774 static int power_components_tarjan(struct basic_map_sort *s,
 775         __isl_keep isl_map *map, int i)
 776 {
 777         int j;
 778
 779         s->node[i].index = s->index;
 780         s->node[i].min_index = s->index;
 781         s->node[i].on_stack = 1;
 782         s->index++;
 783         s->stack[s->sp++] = i;
 784
 785         for (j = s->len - 1; j >= 0; --j) {
 786                 int f;
 787
 788                 if (j == i)
 789                         continue;
 790                 if (s->node[j].index >= 0 &&
 791                         (!s->node[j].on_stack ||
 792                          s->node[j].index > s->node[i].min_index))
 793                         continue;
 794
 795                 f = basic_map_follows(map->p[i], map->p[j]);
 796                 if (f < 0)
 797                         return -1;
 798                 if (!f)
 799                         continue;
 800
 801                 if (s->node[j].index < 0) {
 802                         power_components_tarjan(s, map, j);
 803                         if (s->node[j].min_index < s->node[i].min_index)
 804                                 s->node[i].min_index = s->node[j].min_index;
 805                 } else if (s->node[j].index < s->node[i].min_index)
 806                         s->node[i].min_index = s->node[j].index;
 807         }
 808
 809         if (s->node[i].index != s->node[i].min_index)
 810                 return 0;
 811
 812         do {
 813                 j = s->stack[--s->sp];
 814                 s->node[j].on_stack = 0;
 815                 s->order[s->op++] = j;
 816         } while (j != i);
 817         s->order[s->op++] = -1;
 818
 819         return 0;
 820 }
 821
 822 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
 823  * and a dimension specification (Z^{n+1} -> Z^{n+1}),
 824  * construct a map that is the union of the identity map and
 825  * an overapproximation of the map
 826  * that takes an element from the dom R \times Z to an
 827  * element from ran R \times Z, such that the first n coordinates of the
 828  * difference between them is a sum of differences between images
 829  * and pre-images in one of the R_i and such that the last coordinate
 830  * is equal to the number of steps taken.
 831  * That is, let
 832  *
 833  *      \Delta_i = { y - x | (x, y) in R_i }
 834  *
 835  * then the constructed map is an overapproximation of
 836  *
 837  *      { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
 838  *                              d = (\sum_i k_i \delta_i, \sum_i k_i) and
 839  *                              x in dom R and x + d in ran R } union
 840  *      { (x) -> (x) }
 841  *
 842  * We first split the map into strongly connected components, perform
 843  * the above on each component and the join the results in the correct
 844  * order.  The power of each of the components needs to be extended
 845  * with the identity map because a path in the global result need
 846  * not go through every component.
 847  * The final result will then also contain the identity map, but
 848  * this part will be removed when the length of the path is forced
 849  * to be strictly positive.
 850  */
 851 static __isl_give isl_map *construct_power_components(__isl_take isl_dim *dim,
 852         __isl_keep isl_map *map, int *exact, int project)
 853 {
 854         int i, n;
 855         struct isl_map *path = NULL;
 856         struct basic_map_sort *s = NULL;
 857
 858         if (!map)
 859                 goto error;
 860         if (map->n <= 1)
 861                 return construct_component(dim, map, exact, project);
 862
 863         s = basic_map_sort_alloc(map->ctx, map->n);
 864         if (!s)
 865                 goto error;
 866         for (i = map->n - 1; i >= 0; --i) {
 867                 if (s->node[i].index >= 0)
 868                         continue;
 869                 if (power_components_tarjan(s, map, i) < 0)
 870                         goto error;
 871         }
 872
 873         i = 0;
 874         n = map->n;
 875         path = isl_map_identity(isl_dim_domain(isl_dim_copy(dim)));
 876         while (n) {
 877                 struct isl_map *comp;
 878                 comp = isl_map_alloc_dim(isl_map_get_dim(map), n, 0);
 879                 while (s->order[i] != -1) {
 880                         comp = isl_map_add_basic_map(comp,
 881                                     isl_basic_map_copy(map->p[s->order[i]]));
 882                         --n;
 883                         ++i;
 884                 }
 885                 path = isl_map_apply_range(path,
 886                             construct_component(isl_dim_copy(dim), comp,
 887                                                 exact, project));
 888                 isl_map_free(comp);
 889                 ++i;
 890         }
 891
 892         basic_map_sort_free(s);
 893         isl_dim_free(dim);
 894
 895         return path;
 896 error:
 897         basic_map_sort_free(s);
 898         isl_dim_free(dim);
 899         return NULL;
 900 }
 901
 902 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D,
 903  * construct a map that is an overapproximation of the map
 904  * that takes an element from the space D to another
 905  * element from the same space, such that the difference between
 906  * them is a strictly positive sum of differences between images
 907  * and pre-images in one of the R_i.
 908  * The number of differences in the sum is equated to parameter "param".
 909  * That is, let
 910  *
 911  *      \Delta_i = { y - x | (x, y) in R_i }
 912  *
 913  * then the constructed map is an overapproximation of
 914  *
 915  *      { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
 916  *                              d = \sum_i k_i \delta_i and k = \sum_i k_i > 0 }
 917  *
 918  * We first construct an extended mapping with an extra coordinate
 919  * that indicates the number of steps taken.  In particular,
 920  * the difference in the last coordinate is equal to the number
 921  * of steps taken to move from a domain element to the corresponding
 922  * image element(s).
 923  * In the final step, this difference is equated to the parameter "param"
 924  * and made positive.  The extra coordinates are subsequently projected out.
 925  */
 926 static __isl_give isl_map *construct_power(__isl_keep isl_map *map,
 927         unsigned param, int *exact, int project)
 928 {
 929         struct isl_map *app = NULL;
 930         struct isl_map *diff;
 931         struct isl_dim *dim = NULL;
 932         unsigned d;
 933
 934         if (!map)
 935                 return NULL;
 936
 937         dim = isl_map_get_dim(map);
 938
 939         d = isl_dim_size(dim, isl_dim_in);
 940         dim = isl_dim_add(dim, isl_dim_in, 1);
 941         dim = isl_dim_add(dim, isl_dim_out, 1);
 942
 943         app = construct_power_components(isl_dim_copy(dim), map,
 944                                         exact, project);
 945
 946         diff = equate_parameter_to_length(dim, param);
 947         app = isl_map_intersect(app, diff);
 948         app = isl_map_project_out(app, isl_dim_in, d, 1);
 949         app = isl_map_project_out(app, isl_dim_out, d, 1);
 950
 951         return app;
 952 }
 953
 954 /* Compute the positive powers of "map", or an overapproximation.
 955  * The power is given by parameter "param".  If the result is exact,
 956  * then *exact is set to 1.
 957  * If project is set, then we are actually interested in the transitive
 958  * closure, so we can use a more relaxed exactness check.
 959  */
 960 static __isl_give isl_map *map_power(__isl_take isl_map *map, unsigned param,
 961         int *exact, int project)
 962 {
 963         struct isl_map *app = NULL;
 964
 965         if (exact)
 966                 *exact = 1;
 967
 968         map = isl_map_remove_empty_parts(map);
 969         if (!map)
 970                 return NULL;
 971
 972         if (isl_map_fast_is_empty(map))
 973                 return map;
 974
 975         isl_assert(map->ctx, param < isl_map_dim(map, isl_dim_param), goto error);
 976         isl_assert(map->ctx,
 977                 isl_map_dim(map, isl_dim_in) == isl_map_dim(map, isl_dim_out),
 978                 goto error);
 979
 980         app = construct_power(map, param, exact, project);
 981
 982         isl_map_free(map);
 983         return app;
 984 error:
 985         isl_map_free(map);
 986         isl_map_free(app);
 987         return NULL;
 988 }
 989
 990 /* Compute the positive powers of "map", or an overapproximation.
 991  * The power is given by parameter "param".  If the result is exact,
 992  * then *exact is set to 1.
 993  */
 994 __isl_give isl_map *isl_map_power(__isl_take isl_map *map, unsigned param,
 995         int *exact)
 996 {
 997         return map_power(map, param, exact, 0);
 998 }
 999
1000 /* Compute the transitive closure  of "map", or an overapproximation.
1001  * If the result is exact, then *exact is set to 1.
1002  * Simply compute the powers of map and then project out the parameter
1003  * describing the power.
1004  */
1005 __isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
1006         int *exact)
1007 {
1008         unsigned param;
1009
1010         if (!map)
1011                 goto error;
1012
1013         param = isl_map_dim(map, isl_dim_param);
1014         map = isl_map_add(map, isl_dim_param, 1);
1015         map = map_power(map, param, exact, 1);
1016         map = isl_map_project_out(map, isl_dim_param, param, 1);
1017
1018         return map;
1019 error:
1020         isl_map_free(map);
1021         return NULL;
1022 }