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