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 path with the as yet unconstrained length at position "pos",
 298  * check if setting the length to zero results in only the identity
 299  * mapping.
 300  */
 301 int empty_path_is_identity(__isl_keep isl_basic_map *path, unsigned pos)
 302 {
 303         isl_basic_map *test = NULL;
 304         isl_basic_map *id = NULL;
 305         int k;
 306         int is_id;
 307
 308         test = isl_basic_map_copy(path);
 309         test = isl_basic_map_extend_constraints(test, 1, 0);
 310         k = isl_basic_map_alloc_equality(test);
 311         if (k < 0)
 312                 goto error;
 313         isl_seq_clr(test->eq[k], 1 + isl_basic_map_total_dim(test));
 314         isl_int_set_si(test->eq[k][pos], 1);
 315         id = isl_basic_map_identity(isl_dim_domain(isl_basic_map_get_dim(path)));
 316         is_id = isl_basic_map_is_subset(test, id);
 317         isl_basic_map_free(test);
 318         isl_basic_map_free(id);
 319         return is_id;
 320 error:
 321         isl_basic_map_free(test);
 322         return -1;
 323 }
 324
 325 __isl_give isl_basic_map *add_delta_constraints(__isl_take isl_basic_map *path,
 326         __isl_keep isl_basic_set *delta, unsigned off, unsigned nparam,
 327         unsigned d, int eq)
 328 {
 329         int i, k;
 330         int n = eq ? delta->n_eq : delta->n_ineq;
 331         isl_int **delta_c = eq ? delta->eq : delta->ineq;
 332         isl_int **path_c = eq ? path->eq : path->ineq;
 333
 334         for (i = 0; i < n; ++i) {
 335                 int p = purity(delta, delta_c[i], eq);
 336                 if (p < 0)
 337                         goto error;
 338                 if (p == IMPURE)
 339                         continue;
 340                 if (eq)
 341                         k = isl_basic_map_alloc_equality(path);
 342                 else
 343                         k = isl_basic_map_alloc_inequality(path);
 344                 if (k < 0)
 345                         goto error;
 346                 isl_seq_clr(path_c[k], 1 + isl_basic_map_total_dim(path));
 347                 if (p == PURE_VAR) {
 348                         isl_seq_cpy(path_c[k] + off,
 349                                     delta_c[i] + 1 + nparam, d);
 350                         isl_int_set(path_c[k][off + d], delta_c[i][0]);
 351                 } else if (p == PURE_PARAM) {
 352                         isl_seq_cpy(path_c[k], delta_c[i], 1 + nparam);
 353                 } else {
 354                         isl_seq_cpy(path_c[k] + off,
 355                                     delta_c[i] + 1 + nparam, d);
 356                         isl_seq_cpy(path_c[k], delta_c[i], 1 + nparam);
 357                 }
 358         }
 359
 360         return path;
 361 error:
 362         isl_basic_map_free(path);
 363         return NULL;
 364 }
 365
 366 /* Given a set of offsets "delta", construct a relation of the
 367  * given dimension specification (Z^{n+1} -> Z^{n+1}) that
 368  * is an overapproximation of the relations that
 369  * maps an element x to any element that can be reached
 370  * by taking a non-negative number of steps along any of
 371  * the elements in "delta".
 372  * That is, construct an approximation of
 373  *
 374  *      { [x] -> [y] : exists f \in \delta, k \in Z :
 375  *                                      y = x + k [f, 1] and k >= 0 }
 376  *
 377  * For any element in this relation, the number of steps taken
 378  * is equal to the difference in the final coordinates.
 379  *
 380  * In particular, let delta be defined as
 381  *
 382  *      \delta = [p] -> { [x] : A x + a >= and B p + b >= 0 and
 383  *                              C x + C'p + c >= 0 and
 384  *                              D x + D'p + d >= 0 }
 385  *
 386  * where the constraints C x + C'p + c >= 0 are such that the parametric
 387  * constant term of each constraint j, "C_j x + C'_j p + c_j",
 388  * can never attain positive values, then the relation is constructed as
 389  *
 390  *      { [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and
 391  *                      A f + k a >= 0 and B p + b >= 0 and
 392  *                      C f + C'p + c >= 0 and k >= 1 }
 393  *      union { [x] -> [x] }
 394  *
 395  * If the zero-length paths happen to correspond exactly to the identity
 396  * mapping, then we return
 397  *
 398  *      { [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and
 399  *                      A f + k a >= 0 and B p + b >= 0 and
 400  *                      C f + C'p + c >= 0 and k >= 0 }
 401  *
 402  * instead.
 403  *
 404  * Existentially quantified variables in \delta are currently ignored.
 405  * This is safe, but leads to an additional overapproximation.
 406  */
 407 static __isl_give isl_map *path_along_delta(__isl_take isl_dim *dim,
 408         __isl_take isl_basic_set *delta)
 409 {
 410         isl_basic_map *path = NULL;
 411         unsigned d;
 412         unsigned n_div;
 413         unsigned nparam;
 414         unsigned off;
 415         int i, k;
 416         int is_id;
 417
 418         if (!delta)
 419                 goto error;
 420         n_div = isl_basic_set_dim(delta, isl_dim_div);
 421         d = isl_basic_set_dim(delta, isl_dim_set);
 422         nparam = isl_basic_set_dim(delta, isl_dim_param);
 423         path = isl_basic_map_alloc_dim(isl_dim_copy(dim), n_div + d + 1,
 424                         d + 1 + delta->n_eq, delta->n_ineq + 1);
 425         off = 1 + nparam + 2 * (d + 1) + n_div;
 426
 427         for (i = 0; i < n_div + d + 1; ++i) {
 428                 k = isl_basic_map_alloc_div(path);
 429                 if (k < 0)
 430                         goto error;
 431                 isl_int_set_si(path->div[k][0], 0);
 432         }
 433
 434         for (i = 0; i < d + 1; ++i) {
 435                 k = isl_basic_map_alloc_equality(path);
 436                 if (k < 0)
 437                         goto error;
 438                 isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
 439                 isl_int_set_si(path->eq[k][1 + nparam + i], 1);
 440                 isl_int_set_si(path->eq[k][1 + nparam + d + 1 + i], -1);
 441                 isl_int_set_si(path->eq[k][off + i], 1);
 442         }
 443
 444         path = add_delta_constraints(path, delta, off, nparam, d, 1);
 445         path = add_delta_constraints(path, delta, off, nparam, d, 0);
 446
 447         is_id = empty_path_is_identity(path, off + d);
 448         if (is_id < 0)
 449                 goto error;
 450
 451         k = isl_basic_map_alloc_inequality(path);
 452         if (k < 0)
 453                 goto error;
 454         isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
 455         if (!is_id)
 456                 isl_int_set_si(path->ineq[k][0], -1);
 457         isl_int_set_si(path->ineq[k][off + d], 1);
 458
 459         isl_basic_set_free(delta);
 460         path = isl_basic_map_finalize(path);
 461         if (is_id) {
 462                 isl_dim_free(dim);
 463                 return isl_map_from_basic_map(path);
 464         }
 465         return isl_basic_map_union(path,
 466                                 isl_basic_map_identity(isl_dim_domain(dim)));
 467 error:
 468         isl_dim_free(dim);
 469         isl_basic_set_free(delta);
 470         isl_basic_map_free(path);
 471         return NULL;
 472 }
 473
 474 /* Given a dimenion specification Z^{n+1} -> Z^{n+1} and a parameter "param",
 475  * construct a map that equates the parameter to the difference
 476  * in the final coordinates and imposes that this difference is positive.
 477  * That is, construct
 478  *
 479  *      { [x,x_s] -> [y,y_s] : k = y_s - x_s > 0 }
 480  */
 481 static __isl_give isl_map *equate_parameter_to_length(__isl_take isl_dim *dim,
 482         unsigned param)
 483 {
 484         struct isl_basic_map *bmap;
 485         unsigned d;
 486         unsigned nparam;
 487         int k;
 488
 489         d = isl_dim_size(dim, isl_dim_in);
 490         nparam = isl_dim_size(dim, isl_dim_param);
 491         bmap = isl_basic_map_alloc_dim(dim, 0, 1, 1);
 492         k = isl_basic_map_alloc_equality(bmap);
 493         if (k < 0)
 494                 goto error;
 495         isl_seq_clr(bmap->eq[k], 1 + isl_basic_map_total_dim(bmap));
 496         isl_int_set_si(bmap->eq[k][1 + param], -1);
 497         isl_int_set_si(bmap->eq[k][1 + nparam + d - 1], -1);
 498         isl_int_set_si(bmap->eq[k][1 + nparam + d + d - 1], 1);
 499
 500         k = isl_basic_map_alloc_inequality(bmap);
 501         if (k < 0)
 502                 goto error;
 503         isl_seq_clr(bmap->ineq[k], 1 + isl_basic_map_total_dim(bmap));
 504         isl_int_set_si(bmap->ineq[k][1 + param], 1);
 505         isl_int_set_si(bmap->ineq[k][0], -1);
 506
 507         bmap = isl_basic_map_finalize(bmap);
 508         return isl_map_from_basic_map(bmap);
 509 error:
 510         isl_basic_map_free(bmap);
 511         return NULL;
 512 }
 513
 514 /* Check whether "path" is acyclic, where the last coordinates of domain
 515  * and range of path encode the number of steps taken.
 516  * That is, check whether
 517  *
 518  *      { d | d = y - x and (x,y) in path }
 519  *
 520  * does not contain any element with positive last coordinate (positive length)
 521  * and zero remaining coordinates (cycle).
 522  */
 523 static int is_acyclic(__isl_take isl_map *path)
 524 {
 525         int i;
 526         int acyclic;
 527         unsigned dim;
 528         struct isl_set *delta;
 529
 530         delta = isl_map_deltas(path);
 531         dim = isl_set_dim(delta, isl_dim_set);
 532         for (i = 0; i < dim; ++i) {
 533                 if (i == dim -1)
 534                         delta = isl_set_lower_bound_si(delta, isl_dim_set, i, 1);
 535                 else
 536                         delta = isl_set_fix_si(delta, isl_dim_set, i, 0);
 537         }
 538
 539         acyclic = isl_set_is_empty(delta);
 540         isl_set_free(delta);
 541
 542         return acyclic;
 543 }
 544
 545 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
 546  * and a dimension specification (Z^{n+1} -> Z^{n+1}),
 547  * construct a map that is an overapproximation of the map
 548  * that takes an element from the space D \times Z to another
 549  * element from the same space, such that the first n coordinates of the
 550  * difference between them is a sum of differences between images
 551  * and pre-images in one of the R_i and such that the last coordinate
 552  * is equal to the number of steps taken.
 553  * That is, let
 554  *
 555  *      \Delta_i = { y - x | (x, y) in R_i }
 556  *
 557  * then the constructed map is an overapproximation of
 558  *
 559  *      { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
 560  *                              d = (\sum_i k_i \delta_i, \sum_i k_i) }
 561  *
 562  * The elements of the singleton \Delta_i's are collected as the
 563  * rows of the steps matrix.  For all these \Delta_i's together,
 564  * a single path is constructed.
 565  * For each of the other \Delta_i's, we compute an overapproximation
 566  * of the paths along elements of \Delta_i.
 567  * Since each of these paths performs an addition, composition is
 568  * symmetric and we can simply compose all resulting paths in any order.
 569  */
 570 static __isl_give isl_map *construct_extended_path(__isl_take isl_dim *dim,
 571         __isl_keep isl_map *map, int *project)
 572 {
 573         struct isl_mat *steps = NULL;
 574         struct isl_map *path = NULL;
 575         unsigned d;
 576         int i, j, n;
 577
 578         d = isl_map_dim(map, isl_dim_in);
 579
 580         path = isl_map_identity(isl_dim_domain(isl_dim_copy(dim)));
 581
 582         steps = isl_mat_alloc(map->ctx, map->n, d);
 583         if (!steps)
 584                 goto error;
 585
 586         n = 0;
 587         for (i = 0; i < map->n; ++i) {
 588                 struct isl_basic_set *delta;
 589
 590                 delta = isl_basic_map_deltas(isl_basic_map_copy(map->p[i]));
 591
 592                 for (j = 0; j < d; ++j) {
 593                         int fixed;
 594
 595                         fixed = isl_basic_set_fast_dim_is_fixed(delta, j,
 596                                                             &steps->row[n][j]);
 597                         if (fixed < 0) {
 598                                 isl_basic_set_free(delta);
 599                                 goto error;
 600                         }
 601                         if (!fixed)
 602                                 break;
 603                 }
 604
 605
 606                 if (j < d) {
 607                         path = isl_map_apply_range(path,
 608                                 path_along_delta(isl_dim_copy(dim), delta));
 609                         path = isl_map_coalesce(path);
 610                 } else {
 611                         isl_basic_set_free(delta);
 612                         ++n;
 613                 }
 614         }
 615
 616         if (n > 0) {
 617                 steps->n_row = n;
 618                 path = isl_map_apply_range(path,
 619                                 path_along_steps(isl_dim_copy(dim), steps));
 620         }
 621
 622         if (project && *project) {
 623                 *project = is_acyclic(isl_map_copy(path));
 624                 if (*project < 0)
 625                         goto error;
 626         }
 627
 628         isl_dim_free(dim);
 629         isl_mat_free(steps);
 630         return path;
 631 error:
 632         isl_dim_free(dim);
 633         isl_mat_free(steps);
 634         isl_map_free(path);
 635         return NULL;
 636 }
 637
 638 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
 639  * and a dimension specification (Z^{n+1} -> Z^{n+1}),
 640  * construct a map that is the union of the identity map and
 641  * an overapproximation of the map
 642  * that takes an element from the dom R \times Z to an
 643  * element from ran R \times Z, such that the first n coordinates of the
 644  * difference between them is a sum of differences between images
 645  * and pre-images in one of the R_i and such that the last coordinate
 646  * is equal to the number of steps taken.
 647  * That is, let
 648  *
 649  *      \Delta_i = { y - x | (x, y) in R_i }
 650  *
 651  * then the constructed map is an overapproximation of
 652  *
 653  *      { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
 654  *                              d = (\sum_i k_i \delta_i, \sum_i k_i) and
 655  *                              x in dom R and x + d in ran R } union
 656  *      { (x) -> (x) }
 657  */
 658 static __isl_give isl_map *construct_component(__isl_take isl_dim *dim,
 659         __isl_keep isl_map *map, int *exact, int project)
 660 {
 661         struct isl_set *domain = NULL;
 662         struct isl_set *range = NULL;
 663         struct isl_set *overlap;
 664         struct isl_map *app = NULL;
 665         struct isl_map *path = NULL;
 666
 667         domain = isl_map_domain(isl_map_copy(map));
 668         domain = isl_set_coalesce(domain);
 669         range = isl_map_range(isl_map_copy(map));
 670         range = isl_set_coalesce(range);
 671         overlap = isl_set_intersect(isl_set_copy(domain), isl_set_copy(range));
 672         if (isl_set_is_empty(overlap) == 1) {
 673                 isl_set_free(domain);
 674                 isl_set_free(range);
 675                 isl_set_free(overlap);
 676                 isl_dim_free(dim);
 677
 678                 map = isl_map_copy(map);
 679                 map = isl_map_add(map, isl_dim_in, 1);
 680                 map = isl_map_add(map, isl_dim_out, 1);
 681                 map = set_path_length(map, 1, 1);
 682                 return map;
 683         }
 684         isl_set_free(overlap);
 685         app = isl_map_from_domain_and_range(domain, range);
 686         app = isl_map_add(app, isl_dim_in, 1);
 687         app = isl_map_add(app, isl_dim_out, 1);
 688
 689         path = construct_extended_path(isl_dim_copy(dim), map,
 690                                         exact && *exact ? &project : NULL);
 691         app = isl_map_intersect(app, path);
 692
 693         if (exact && *exact &&
 694             (*exact = check_exactness(isl_map_copy(map), isl_map_copy(app),
 695                                       project)) < 0)
 696                 goto error;
 697
 698         return isl_map_union(app, isl_map_identity(isl_dim_domain(dim)));
 699 error:
 700         isl_dim_free(dim);
 701         isl_map_free(app);
 702         return NULL;
 703 }
 704
 705 /* Structure for representing the nodes in the graph being traversed
 706  * using Tarjan's algorithm.
 707  * index represents the order in which nodes are visited.
 708  * min_index is the index of the root of a (sub)component.
 709  * on_stack indicates whether the node is currently on the stack.
 710  */
 711 struct basic_map_sort_node {
 712         int index;
 713         int min_index;
 714         int on_stack;
 715 };
 716 /* Structure for representing the graph being traversed
 717  * using Tarjan's algorithm.
 718  * len is the number of nodes
 719  * node is an array of nodes
 720  * stack contains the nodes on the path from the root to the current node
 721  * sp is the stack pointer
 722  * index is the index of the last node visited
 723  * order contains the elements of the components separated by -1
 724  * op represents the current position in order
 725  */
 726 struct basic_map_sort {
 727         int len;
 728         struct basic_map_sort_node *node;
 729         int *stack;
 730         int sp;
 731         int index;
 732         int *order;
 733         int op;
 734 };
 735
 736 static void basic_map_sort_free(struct basic_map_sort *s)
 737 {
 738         if (!s)
 739                 return;
 740         free(s->node);
 741         free(s->stack);
 742         free(s->order);
 743         free(s);
 744 }
 745
 746 static struct basic_map_sort *basic_map_sort_alloc(struct isl_ctx *ctx, int len)
 747 {
 748         struct basic_map_sort *s;
 749         int i;
 750
 751         s = isl_calloc_type(ctx, struct basic_map_sort);
 752         if (!s)
 753                 return NULL;
 754         s->len = len;
 755         s->node = isl_alloc_array(ctx, struct basic_map_sort_node, len);
 756         if (!s->node)
 757                 goto error;
 758         for (i = 0; i < len; ++i)
 759                 s->node[i].index = -1;
 760         s->stack = isl_alloc_array(ctx, int, len);
 761         if (!s->stack)
 762                 goto error;
 763         s->order = isl_alloc_array(ctx, int, 2 * len);
 764         if (!s->order)
 765                 goto error;
 766
 767         s->sp = 0;
 768         s->index = 0;
 769         s->op = 0;
 770
 771         return s;
 772 error:
 773         basic_map_sort_free(s);
 774         return NULL;
 775 }
 776
 777 /* Check whether in the computation of the transitive closure
 778  * "bmap1" (R_1) should follow (or be part of the same component as)
 779  * "bmap2" (R_2).
 780  *
 781  * That is check whether
 782  *
 783  *      R_1 \circ R_2
 784  *
 785  * is a subset of
 786  *
 787  *      R_2 \circ R_1
 788  *
 789  * If so, then there is no reason for R_1 to immediately follow R_2
 790  * in any path.
 791  */
 792 static int basic_map_follows(__isl_keep isl_basic_map *bmap1,
 793         __isl_keep isl_basic_map *bmap2)
 794 {
 795         struct isl_map *map12 = NULL;
 796         struct isl_map *map21 = NULL;
 797         int subset;
 798
 799         map21 = isl_map_from_basic_map(
 800                         isl_basic_map_apply_range(
 801                                 isl_basic_map_copy(bmap2),
 802                                 isl_basic_map_copy(bmap1)));
 803         subset = isl_map_is_empty(map21);
 804         if (subset < 0)
 805                 goto error;
 806         if (subset) {
 807                 isl_map_free(map21);
 808                 return 0;
 809         }
 810
 811         map12 = isl_map_from_basic_map(
 812                         isl_basic_map_apply_range(
 813                                 isl_basic_map_copy(bmap1),
 814                                 isl_basic_map_copy(bmap2)));
 815
 816         subset = isl_map_is_subset(map21, map12);
 817
 818         isl_map_free(map12);
 819         isl_map_free(map21);
 820
 821         return subset < 0 ? -1 : !subset;
 822 error:
 823         isl_map_free(map21);
 824         return -1;
 825 }
 826
 827 /* Perform Tarjan's algorithm for computing the strongly connected components
 828  * in the graph with the disjuncts of "map" as vertices and with an
 829  * edge between any pair of disjuncts such that the first has
 830  * to be applied after the second.
 831  */
 832 static int power_components_tarjan(struct basic_map_sort *s,
 833         __isl_keep isl_map *map, int i)
 834 {
 835         int j;
 836
 837         s->node[i].index = s->index;
 838         s->node[i].min_index = s->index;
 839         s->node[i].on_stack = 1;
 840         s->index++;
 841         s->stack[s->sp++] = i;
 842
 843         for (j = s->len - 1; j >= 0; --j) {
 844                 int f;
 845
 846                 if (j == i)
 847                         continue;
 848                 if (s->node[j].index >= 0 &&
 849                         (!s->node[j].on_stack ||
 850                          s->node[j].index > s->node[i].min_index))
 851                         continue;
 852
 853                 f = basic_map_follows(map->p[i], map->p[j]);
 854                 if (f < 0)
 855                         return -1;
 856                 if (!f)
 857                         continue;
 858
 859                 if (s->node[j].index < 0) {
 860                         power_components_tarjan(s, map, j);
 861                         if (s->node[j].min_index < s->node[i].min_index)
 862                                 s->node[i].min_index = s->node[j].min_index;
 863                 } else if (s->node[j].index < s->node[i].min_index)
 864                         s->node[i].min_index = s->node[j].index;
 865         }
 866
 867         if (s->node[i].index != s->node[i].min_index)
 868                 return 0;
 869
 870         do {
 871                 j = s->stack[--s->sp];
 872                 s->node[j].on_stack = 0;
 873                 s->order[s->op++] = j;
 874         } while (j != i);
 875         s->order[s->op++] = -1;
 876
 877         return 0;
 878 }
 879
 880 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
 881  * and a dimension specification (Z^{n+1} -> Z^{n+1}),
 882  * construct a map that is the union of the identity map and
 883  * an overapproximation of the map
 884  * that takes an element from the dom R \times Z to an
 885  * element from ran R \times Z, such that the first n coordinates of the
 886  * difference between them is a sum of differences between images
 887  * and pre-images in one of the R_i and such that the last coordinate
 888  * is equal to the number of steps taken.
 889  * That is, let
 890  *
 891  *      \Delta_i = { y - x | (x, y) in R_i }
 892  *
 893  * then the constructed map is an overapproximation of
 894  *
 895  *      { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
 896  *                              d = (\sum_i k_i \delta_i, \sum_i k_i) and
 897  *                              x in dom R and x + d in ran R } union
 898  *      { (x) -> (x) }
 899  *
 900  * We first split the map into strongly connected components, perform
 901  * the above on each component and the join the results in the correct
 902  * order.  The power of each of the components needs to be extended
 903  * with the identity map because a path in the global result need
 904  * not go through every component.
 905  * The final result will then also contain the identity map, but
 906  * this part will be removed when the length of the path is forced
 907  * to be strictly positive.
 908  */
 909 static __isl_give isl_map *construct_power_components(__isl_take isl_dim *dim,
 910         __isl_keep isl_map *map, int *exact, int project)
 911 {
 912         int i, n;
 913         struct isl_map *path = NULL;
 914         struct basic_map_sort *s = NULL;
 915
 916         if (!map)
 917                 goto error;
 918         if (map->n <= 1)
 919                 return construct_component(dim, map, exact, project);
 920
 921         s = basic_map_sort_alloc(map->ctx, map->n);
 922         if (!s)
 923                 goto error;
 924         for (i = map->n - 1; i >= 0; --i) {
 925                 if (s->node[i].index >= 0)
 926                         continue;
 927                 if (power_components_tarjan(s, map, i) < 0)
 928                         goto error;
 929         }
 930
 931         i = 0;
 932         n = map->n;
 933         path = isl_map_identity(isl_dim_domain(isl_dim_copy(dim)));
 934         while (n) {
 935                 struct isl_map *comp;
 936                 comp = isl_map_alloc_dim(isl_map_get_dim(map), n, 0);
 937                 while (s->order[i] != -1) {
 938                         comp = isl_map_add_basic_map(comp,
 939                                     isl_basic_map_copy(map->p[s->order[i]]));
 940                         --n;
 941                         ++i;
 942                 }
 943                 path = isl_map_apply_range(path,
 944                             construct_component(isl_dim_copy(dim), comp,
 945                                                 exact, project));
 946                 isl_map_free(comp);
 947                 ++i;
 948         }
 949
 950         basic_map_sort_free(s);
 951         isl_dim_free(dim);
 952
 953         return path;
 954 error:
 955         basic_map_sort_free(s);
 956         isl_dim_free(dim);
 957         return NULL;
 958 }
 959
 960 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D,
 961  * construct a map that is an overapproximation of the map
 962  * that takes an element from the space D to another
 963  * element from the same space, such that the difference between
 964  * them is a strictly positive sum of differences between images
 965  * and pre-images in one of the R_i.
 966  * The number of differences in the sum is equated to parameter "param".
 967  * That is, let
 968  *
 969  *      \Delta_i = { y - x | (x, y) in R_i }
 970  *
 971  * then the constructed map is an overapproximation of
 972  *
 973  *      { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
 974  *                              d = \sum_i k_i \delta_i and k = \sum_i k_i > 0 }
 975  *
 976  * We first construct an extended mapping with an extra coordinate
 977  * that indicates the number of steps taken.  In particular,
 978  * the difference in the last coordinate is equal to the number
 979  * of steps taken to move from a domain element to the corresponding
 980  * image element(s).
 981  * In the final step, this difference is equated to the parameter "param"
 982  * and made positive.  The extra coordinates are subsequently projected out.
 983  */
 984 static __isl_give isl_map *construct_power(__isl_keep isl_map *map,
 985         unsigned param, int *exact, int project)
 986 {
 987         struct isl_map *app = NULL;
 988         struct isl_map *diff;
 989         struct isl_dim *dim = NULL;
 990         unsigned d;
 991
 992         if (!map)
 993                 return NULL;
 994
 995         dim = isl_map_get_dim(map);
 996
 997         d = isl_dim_size(dim, isl_dim_in);
 998         dim = isl_dim_add(dim, isl_dim_in, 1);
 999         dim = isl_dim_add(dim, isl_dim_out, 1);
1000
1001         app = construct_power_components(isl_dim_copy(dim), map,
1002                                         exact, project);
1003
1004         diff = equate_parameter_to_length(dim, param);
1005         app = isl_map_intersect(app, diff);
1006         app = isl_map_project_out(app, isl_dim_in, d, 1);
1007         app = isl_map_project_out(app, isl_dim_out, d, 1);
1008
1009         return app;
1010 }
1011
1012 /* Compute the positive powers of "map", or an overapproximation.
1013  * The power is given by parameter "param".  If the result is exact,
1014  * then *exact is set to 1.
1015  * If project is set, then we are actually interested in the transitive
1016  * closure, so we can use a more relaxed exactness check.
1017  */
1018 static __isl_give isl_map *map_power(__isl_take isl_map *map, unsigned param,
1019         int *exact, int project)
1020 {
1021         struct isl_map *app = NULL;
1022
1023         if (exact)
1024                 *exact = 1;
1025
1026         map = isl_map_remove_empty_parts(map);
1027         if (!map)
1028                 return NULL;
1029
1030         if (isl_map_fast_is_empty(map))
1031                 return map;
1032
1033         isl_assert(map->ctx, param < isl_map_dim(map, isl_dim_param), goto error);
1034         isl_assert(map->ctx,
1035                 isl_map_dim(map, isl_dim_in) == isl_map_dim(map, isl_dim_out),
1036                 goto error);
1037
1038         app = construct_power(map, param, exact, project);
1039
1040         isl_map_free(map);
1041         return app;
1042 error:
1043         isl_map_free(map);
1044         isl_map_free(app);
1045         return NULL;
1046 }
1047
1048 /* Compute the positive powers of "map", or an overapproximation.
1049  * The power is given by parameter "param".  If the result is exact,
1050  * then *exact is set to 1.
1051  */
1052 __isl_give isl_map *isl_map_power(__isl_take isl_map *map, unsigned param,
1053         int *exact)
1054 {
1055         return map_power(map, param, exact, 0);
1056 }
1057
1058 /* Compute the transitive closure  of "map", or an overapproximation.
1059  * If the result is exact, then *exact is set to 1.
1060  * Simply compute the powers of map and then project out the parameter
1061  * describing the power.
1062  */
1063 __isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
1064         int *exact)
1065 {
1066         unsigned param;
1067
1068         if (!map)
1069                 goto error;
1070
1071         param = isl_map_dim(map, isl_dim_param);
1072         map = isl_map_add(map, isl_dim_param, 1);
1073         map = map_power(map, param, exact, 1);
1074         map = isl_map_project_out(map, isl_dim_param, param, 1);
1075
1076         return map;
1077 error:
1078         isl_map_free(map);
1079         return NULL;
1080 }