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