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