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 #include <isl_lp.h>
  15
  16 /* Given a map that represents a path with the length of the path
  17  * encoded as the difference between the last output coordindate
  18  * and the last input coordinate, set this length to either
  19  * exactly "length" (if "exactly" is set) or at least "length"
  20  * (if "exactly" is not set).
  21  */
  22 static __isl_give isl_map *set_path_length(__isl_take isl_map *map,
  23         int exactly, int length)
  24 {
  25         struct isl_dim *dim;
  26         struct isl_basic_map *bmap;
  27         unsigned d;
  28         unsigned nparam;
  29         int k;
  30         isl_int *c;
  31
  32         if (!map)
  33                 return NULL;
  34
  35         dim = isl_map_get_dim(map);
  36         d = isl_dim_size(dim, isl_dim_in);
  37         nparam = isl_dim_size(dim, isl_dim_param);
  38         bmap = isl_basic_map_alloc_dim(dim, 0, 1, 1);
  39         if (exactly) {
  40                 k = isl_basic_map_alloc_equality(bmap);
  41                 c = bmap->eq[k];
  42         } else {
  43                 k = isl_basic_map_alloc_inequality(bmap);
  44                 c = bmap->ineq[k];
  45         }
  46         if (k < 0)
  47                 goto error;
  48         isl_seq_clr(c, 1 + isl_basic_map_total_dim(bmap));
  49         isl_int_set_si(c[0], -length);
  50         isl_int_set_si(c[1 + nparam + d - 1], -1);
  51         isl_int_set_si(c[1 + nparam + d + d - 1], 1);
  52
  53         bmap = isl_basic_map_finalize(bmap);
  54         map = isl_map_intersect(map, isl_map_from_basic_map(bmap));
  55
  56         return map;
  57 error:
  58         isl_basic_map_free(bmap);
  59         isl_map_free(map);
  60         return NULL;
  61 }
  62
  63 /* Check whether the overapproximation of the power of "map" is exactly
  64  * the power of "map".  Let R be "map" and A_k the overapproximation.
  65  * The approximation is exact if
  66  *
  67  *      A_1 = R
  68  *      A_k = A_{k-1} \circ R                   k >= 2
  69  *
  70  * Since A_k is known to be an overapproximation, we only need to check
  71  *
  72  *      A_1 \subset R
  73  *      A_k \subset A_{k-1} \circ R             k >= 2
  74  *
  75  * In practice, "app" has an extra input and output coordinate
  76  * to encode the length of the path.  So, we first need to add
  77  * this coordinate to "map" and set the length of the path to
  78  * one.
  79  */
  80 static int check_power_exactness(__isl_take isl_map *map,
  81         __isl_take isl_map *app)
  82 {
  83         int exact;
  84         isl_map *app_1;
  85         isl_map *app_2;
  86
  87         map = isl_map_add(map, isl_dim_in, 1);
  88         map = isl_map_add(map, isl_dim_out, 1);
  89         map = set_path_length(map, 1, 1);
  90
  91         app_1 = set_path_length(isl_map_copy(app), 1, 1);
  92
  93         exact = isl_map_is_subset(app_1, map);
  94         isl_map_free(app_1);
  95
  96         if (!exact || exact < 0) {
  97                 isl_map_free(app);
  98                 isl_map_free(map);
  99                 return exact;
 100         }
 101
 102         app_1 = set_path_length(isl_map_copy(app), 0, 1);
 103         app_2 = set_path_length(app, 0, 2);
 104         app_1 = isl_map_apply_range(map, app_1);
 105
 106         exact = isl_map_is_subset(app_2, app_1);
 107
 108         isl_map_free(app_1);
 109         isl_map_free(app_2);
 110
 111         return exact;
 112 }
 113
 114 /* Check whether the overapproximation of the power of "map" is exactly
 115  * the power of "map", possibly after projecting out the power (if "project"
 116  * is set).
 117  *
 118  * If "project" is set and if "steps" can only result in acyclic paths,
 119  * then we check
 120  *
 121  *      A = R \cup (A \circ R)
 122  *
 123  * where A is the overapproximation with the power projected out, i.e.,
 124  * an overapproximation of the transitive closure.
 125  * More specifically, since A is known to be an overapproximation, we check
 126  *
 127  *      A \subset R \cup (A \circ R)
 128  *
 129  * Otherwise, we check if the power is exact.
 130  *
 131  * Note that "app" has an extra input and output coordinate to encode
 132  * the length of the part.  If we are only interested in the transitive
 133  * closure, then we can simply project out these coordinates first.
 134  */
 135 static int check_exactness(__isl_take isl_map *map, __isl_take isl_map *app,
 136         int project)
 137 {
 138         isl_map *test;
 139         int exact;
 140         unsigned d;
 141
 142         if (!project)
 143                 return check_power_exactness(map, app);
 144
 145         d = isl_map_dim(map, isl_dim_in);
 146         app = set_path_length(app, 0, 1);
 147         app = isl_map_project_out(app, isl_dim_in, d, 1);
 148         app = isl_map_project_out(app, isl_dim_out, d, 1);
 149
 150         test = isl_map_apply_range(isl_map_copy(map), isl_map_copy(app));
 151         test = isl_map_union(test, isl_map_copy(map));
 152
 153         exact = isl_map_is_subset(app, test);
 154
 155         isl_map_free(app);
 156         isl_map_free(test);
 157
 158         isl_map_free(map);
 159
 160         return exact;
 161 }
 162
 163 /*
 164  * The transitive closure implementation is based on the paper
 165  * "Computing the Transitive Closure of a Union of Affine Integer
 166  * Tuple Relations" by Anna Beletska, Denis Barthou, Wlodzimierz Bielecki and
 167  * Albert Cohen.
 168  */
 169
 170 /* Given a set of n offsets v_i (the rows of "steps"), construct a relation
 171  * of the given dimension specification (Z^{n+1} -> Z^{n+1})
 172  * that maps an element x to any element that can be reached
 173  * by taking a non-negative number of steps along any of
 174  * the extended offsets v'_i = [v_i 1].
 175  * That is, construct
 176  *
 177  * { [x] -> [y] : exists k_i >= 0, y = x + \sum_i k_i v'_i }
 178  *
 179  * For any element in this relation, the number of steps taken
 180  * is equal to the difference in the final coordinates.
 181  */
 182 static __isl_give isl_map *path_along_steps(__isl_take isl_dim *dim,
 183         __isl_keep isl_mat *steps)
 184 {
 185         int i, j, k;
 186         struct isl_basic_map *path = NULL;
 187         unsigned d;
 188         unsigned n;
 189         unsigned nparam;
 190
 191         if (!dim || !steps)
 192                 goto error;
 193
 194         d = isl_dim_size(dim, isl_dim_in);
 195         n = steps->n_row;
 196         nparam = isl_dim_size(dim, isl_dim_param);
 197
 198         path = isl_basic_map_alloc_dim(isl_dim_copy(dim), n, d, n);
 199
 200         for (i = 0; i < n; ++i) {
 201                 k = isl_basic_map_alloc_div(path);
 202                 if (k < 0)
 203                         goto error;
 204                 isl_assert(steps->ctx, i == k, goto error);
 205                 isl_int_set_si(path->div[k][0], 0);
 206         }
 207
 208         for (i = 0; i < d; ++i) {
 209                 k = isl_basic_map_alloc_equality(path);
 210                 if (k < 0)
 211                         goto error;
 212                 isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
 213                 isl_int_set_si(path->eq[k][1 + nparam + i], 1);
 214                 isl_int_set_si(path->eq[k][1 + nparam + d + i], -1);
 215                 if (i == d - 1)
 216                         for (j = 0; j < n; ++j)
 217                                 isl_int_set_si(path->eq[k][1 + nparam + 2 * d + j], 1);
 218                 else
 219                         for (j = 0; j < n; ++j)
 220                                 isl_int_set(path->eq[k][1 + nparam + 2 * d + j],
 221                                             steps->row[j][i]);
 222         }
 223
 224         for (i = 0; i < n; ++i) {
 225                 k = isl_basic_map_alloc_inequality(path);
 226                 if (k < 0)
 227                         goto error;
 228                 isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
 229                 isl_int_set_si(path->ineq[k][1 + nparam + 2 * d + i], 1);
 230         }
 231
 232         isl_dim_free(dim);
 233
 234         path = isl_basic_map_simplify(path);
 235         path = isl_basic_map_finalize(path);
 236         return isl_map_from_basic_map(path);
 237 error:
 238         isl_dim_free(dim);
 239         isl_basic_map_free(path);
 240         return NULL;
 241 }
 242
 243 #define IMPURE          0
 244 #define PURE_PARAM      1
 245 #define PURE_VAR        2
 246 #define MIXED           3
 247
 248 /* Check whether the parametric constant term of constraint c is never
 249  * positive in "bset".
 250  */
 251 static int parametric_constant_never_positive(__isl_keep isl_basic_set *bset,
 252         isl_int *c, int *div_purity)
 253 {
 254         unsigned d;
 255         unsigned n_div;
 256         unsigned nparam;
 257         int i;
 258         int k;
 259         int empty;
 260
 261         n_div = isl_basic_set_dim(bset, isl_dim_div);
 262         d = isl_basic_set_dim(bset, isl_dim_set);
 263         nparam = isl_basic_set_dim(bset, isl_dim_param);
 264
 265         bset = isl_basic_set_copy(bset);
 266         bset = isl_basic_set_cow(bset);
 267         bset = isl_basic_set_extend_constraints(bset, 0, 1);
 268         k = isl_basic_set_alloc_inequality(bset);
 269         if (k < 0)
 270                 goto error;
 271         isl_seq_clr(bset->ineq[k], 1 + isl_basic_set_total_dim(bset));
 272         isl_seq_cpy(bset->ineq[k], c, 1 + nparam);
 273         for (i = 0; i < n_div; ++i) {
 274                 if (div_purity[i] != PURE_PARAM)
 275                         continue;
 276                 isl_int_set(bset->ineq[k][1 + nparam + d + i],
 277                             c[1 + nparam + d + i]);
 278         }
 279         isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
 280         empty = isl_basic_set_is_empty(bset);
 281         isl_basic_set_free(bset);
 282
 283         return empty;
 284 error:
 285         isl_basic_set_free(bset);
 286         return -1;
 287 }
 288
 289 /* Return PURE_PARAM if only the coefficients of the parameters are non-zero.
 290  * Return PURE_VAR if only the coefficients of the set variables are non-zero.
 291  * Return MIXED if only the coefficients of the parameters and the set
 292  *      variables are non-zero and if moreover the parametric constant
 293  *      can never attain positive values.
 294  * Return IMPURE otherwise.
 295  */
 296 static int purity(__isl_keep isl_basic_set *bset, isl_int *c, int *div_purity,
 297         int eq)
 298 {
 299         unsigned d;
 300         unsigned n_div;
 301         unsigned nparam;
 302         int empty;
 303         int i;
 304         int p = 0, v = 0;
 305
 306         n_div = isl_basic_set_dim(bset, isl_dim_div);
 307         d = isl_basic_set_dim(bset, isl_dim_set);
 308         nparam = isl_basic_set_dim(bset, isl_dim_param);
 309
 310         for (i = 0; i < n_div; ++i) {
 311                 if (isl_int_is_zero(c[1 + nparam + d + i]))
 312                         continue;
 313                 switch (div_purity[i]) {
 314                 case PURE_PARAM: p = 1; break;
 315                 case PURE_VAR: v = 1; break;
 316                 default: return IMPURE;
 317                 }
 318         }
 319         if (!p && isl_seq_first_non_zero(c + 1, nparam) == -1)
 320                 return PURE_VAR;
 321         if (!v && isl_seq_first_non_zero(c + 1 + nparam, d) == -1)
 322                 return PURE_PARAM;
 323
 324         empty = parametric_constant_never_positive(bset, c, div_purity);
 325         if (eq && empty >= 0 && !empty) {
 326                 isl_seq_neg(c, c, 1 + nparam + d + n_div);
 327                 empty = parametric_constant_never_positive(bset, c, div_purity);
 328         }
 329
 330         return empty < 0 ? -1 : empty ? MIXED : IMPURE;
 331 }
 332
 333 /* Return an array of integers indicating the type of each div in bset.
 334  * If the div is (recursively) defined in terms of only the parameters,
 335  * then the type is PURE_PARAM.
 336  * If the div is (recursively) defined in terms of only the set variables,
 337  * then the type is PURE_VAR.
 338  * Otherwise, the type is IMPURE.
 339  */
 340 static __isl_give int *get_div_purity(__isl_keep isl_basic_set *bset)
 341 {
 342         int i, j;
 343         int *div_purity;
 344         unsigned d;
 345         unsigned n_div;
 346         unsigned nparam;
 347
 348         if (!bset)
 349                 return NULL;
 350
 351         n_div = isl_basic_set_dim(bset, isl_dim_div);
 352         d = isl_basic_set_dim(bset, isl_dim_set);
 353         nparam = isl_basic_set_dim(bset, isl_dim_param);
 354
 355         div_purity = isl_alloc_array(bset->ctx, int, n_div);
 356         if (!div_purity)
 357                 return NULL;
 358
 359         for (i = 0; i < bset->n_div; ++i) {
 360                 int p = 0, v = 0;
 361                 if (isl_int_is_zero(bset->div[i][0])) {
 362                         div_purity[i] = IMPURE;
 363                         continue;
 364                 }
 365                 if (isl_seq_first_non_zero(bset->div[i] + 2, nparam) != -1)
 366                         p = 1;
 367                 if (isl_seq_first_non_zero(bset->div[i] + 2 + nparam, d) != -1)
 368                         v = 1;
 369                 for (j = 0; j < i; ++j) {
 370                         if (isl_int_is_zero(bset->div[i][2 + nparam + d + j]))
 371                                 continue;
 372                         switch (div_purity[j]) {
 373                         case PURE_PARAM: p = 1; break;
 374                         case PURE_VAR: v = 1; break;
 375                         default: p = v = 1; break;
 376                         }
 377                 }
 378                 div_purity[i] = v ? p ? IMPURE : PURE_VAR : PURE_PARAM;
 379         }
 380
 381         return div_purity;
 382 }
 383
 384 /* Given a path with the as yet unconstrained length at position "pos",
 385  * check if setting the length to zero results in only the identity
 386  * mapping.
 387  */
 388 int empty_path_is_identity(__isl_keep isl_basic_map *path, unsigned pos)
 389 {
 390         isl_basic_map *test = NULL;
 391         isl_basic_map *id = NULL;
 392         int k;
 393         int is_id;
 394
 395         test = isl_basic_map_copy(path);
 396         test = isl_basic_map_extend_constraints(test, 1, 0);
 397         k = isl_basic_map_alloc_equality(test);
 398         if (k < 0)
 399                 goto error;
 400         isl_seq_clr(test->eq[k], 1 + isl_basic_map_total_dim(test));
 401         isl_int_set_si(test->eq[k][pos], 1);
 402         id = isl_basic_map_identity(isl_dim_domain(isl_basic_map_get_dim(path)));
 403         is_id = isl_basic_map_is_equal(test, id);
 404         isl_basic_map_free(test);
 405         isl_basic_map_free(id);
 406         return is_id;
 407 error:
 408         isl_basic_map_free(test);
 409         return -1;
 410 }
 411
 412 __isl_give isl_basic_map *add_delta_constraints(__isl_take isl_basic_map *path,
 413         __isl_keep isl_basic_set *delta, unsigned off, unsigned nparam,
 414         unsigned d, int *div_purity, int eq)
 415 {
 416         int i, k;
 417         int n = eq ? delta->n_eq : delta->n_ineq;
 418         isl_int **delta_c = eq ? delta->eq : delta->ineq;
 419         unsigned n_div;
 420
 421         n_div = isl_basic_set_dim(delta, isl_dim_div);
 422
 423         for (i = 0; i < n; ++i) {
 424                 isl_int *path_c;
 425                 int p = purity(delta, delta_c[i], div_purity, eq);
 426                 if (p < 0)
 427                         goto error;
 428                 if (p == IMPURE)
 429                         continue;
 430                 if (eq && p != MIXED) {
 431                         k = isl_basic_map_alloc_equality(path);
 432                         path_c = path->eq[k];
 433                 } else {
 434                         k = isl_basic_map_alloc_inequality(path);
 435                         path_c = path->ineq[k];
 436                 }
 437                 if (k < 0)
 438                         goto error;
 439                 isl_seq_clr(path_c, 1 + isl_basic_map_total_dim(path));
 440                 if (p == PURE_VAR) {
 441                         isl_seq_cpy(path_c + off,
 442                                     delta_c[i] + 1 + nparam, d);
 443                         isl_int_set(path_c[off + d], delta_c[i][0]);
 444                 } else if (p == PURE_PARAM) {
 445                         isl_seq_cpy(path_c, delta_c[i], 1 + nparam);
 446                 } else {
 447                         isl_seq_cpy(path_c + off,
 448                                     delta_c[i] + 1 + nparam, d);
 449                         isl_seq_cpy(path_c, delta_c[i], 1 + nparam);
 450                 }
 451                 isl_seq_cpy(path_c + off - n_div,
 452                             delta_c[i] + 1 + nparam + d, n_div);
 453         }
 454
 455         return path;
 456 error:
 457         isl_basic_map_free(path);
 458         return NULL;
 459 }
 460
 461 /* Given a set of offsets "delta", construct a relation of the
 462  * given dimension specification (Z^{n+1} -> Z^{n+1}) that
 463  * is an overapproximation of the relations that
 464  * maps an element x to any element that can be reached
 465  * by taking a non-negative number of steps along any of
 466  * the elements in "delta".
 467  * That is, construct an approximation of
 468  *
 469  *      { [x] -> [y] : exists f \in \delta, k \in Z :
 470  *                                      y = x + k [f, 1] and k >= 0 }
 471  *
 472  * For any element in this relation, the number of steps taken
 473  * is equal to the difference in the final coordinates.
 474  *
 475  * In particular, let delta be defined as
 476  *
 477  *      \delta = [p] -> { [x] : A x + a >= and B p + b >= 0 and
 478  *                              C x + C'p + c >= 0 and
 479  *                              D x + D'p + d >= 0 }
 480  *
 481  * where the constraints C x + C'p + c >= 0 are such that the parametric
 482  * constant term of each constraint j, "C_j x + C'_j p + c_j",
 483  * can never attain positive values, then the relation is constructed as
 484  *
 485  *      { [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and
 486  *                      A f + k a >= 0 and B p + b >= 0 and
 487  *                      C f + C'p + c >= 0 and k >= 1 }
 488  *      union { [x] -> [x] }
 489  *
 490  * If the zero-length paths happen to correspond exactly to the identity
 491  * mapping, then we return
 492  *
 493  *      { [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and
 494  *                      A f + k a >= 0 and B p + b >= 0 and
 495  *                      C f + C'p + c >= 0 and k >= 0 }
 496  *
 497  * instead.
 498  *
 499  * Existentially quantified variables in \delta are handled by
 500  * classifying them as independent of the parameters, purely
 501  * parameter dependent and others.  Constraints containing
 502  * any of the other existentially quantified variables are removed.
 503  * This is safe, but leads to an additional overapproximation.
 504  */
 505 static __isl_give isl_map *path_along_delta(__isl_take isl_dim *dim,
 506         __isl_take isl_basic_set *delta)
 507 {
 508         isl_basic_map *path = NULL;
 509         unsigned d;
 510         unsigned n_div;
 511         unsigned nparam;
 512         unsigned off;
 513         int i, k;
 514         int is_id;
 515         int *div_purity = NULL;
 516
 517         if (!delta)
 518                 goto error;
 519         n_div = isl_basic_set_dim(delta, isl_dim_div);
 520         d = isl_basic_set_dim(delta, isl_dim_set);
 521         nparam = isl_basic_set_dim(delta, isl_dim_param);
 522         path = isl_basic_map_alloc_dim(isl_dim_copy(dim), n_div + d + 1,
 523                         d + 1 + delta->n_eq, delta->n_eq + delta->n_ineq + 1);
 524         off = 1 + nparam + 2 * (d + 1) + n_div;
 525
 526         for (i = 0; i < n_div + d + 1; ++i) {
 527                 k = isl_basic_map_alloc_div(path);
 528                 if (k < 0)
 529                         goto error;
 530                 isl_int_set_si(path->div[k][0], 0);
 531         }
 532
 533         for (i = 0; i < d + 1; ++i) {
 534                 k = isl_basic_map_alloc_equality(path);
 535                 if (k < 0)
 536                         goto error;
 537                 isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
 538                 isl_int_set_si(path->eq[k][1 + nparam + i], 1);
 539                 isl_int_set_si(path->eq[k][1 + nparam + d + 1 + i], -1);
 540                 isl_int_set_si(path->eq[k][off + i], 1);
 541         }
 542
 543         div_purity = get_div_purity(delta);
 544         if (!div_purity)
 545                 goto error;
 546
 547         path = add_delta_constraints(path, delta, off, nparam, d, div_purity, 1);
 548         path = add_delta_constraints(path, delta, off, nparam, d, div_purity, 0);
 549
 550         is_id = empty_path_is_identity(path, off + d);
 551         if (is_id < 0)
 552                 goto error;
 553
 554         k = isl_basic_map_alloc_inequality(path);
 555         if (k < 0)
 556                 goto error;
 557         isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
 558         if (!is_id)
 559                 isl_int_set_si(path->ineq[k][0], -1);
 560         isl_int_set_si(path->ineq[k][off + d], 1);
 561
 562         free(div_purity);
 563         isl_basic_set_free(delta);
 564         path = isl_basic_map_finalize(path);
 565         if (is_id) {
 566                 isl_dim_free(dim);
 567                 return isl_map_from_basic_map(path);
 568         }
 569         return isl_basic_map_union(path,
 570                                 isl_basic_map_identity(isl_dim_domain(dim)));
 571 error:
 572         free(div_purity);
 573         isl_dim_free(dim);
 574         isl_basic_set_free(delta);
 575         isl_basic_map_free(path);
 576         return NULL;
 577 }
 578
 579 /* Given a dimenion specification Z^{n+1} -> Z^{n+1} and a parameter "param",
 580  * construct a map that equates the parameter to the difference
 581  * in the final coordinates and imposes that this difference is positive.
 582  * That is, construct
 583  *
 584  *      { [x,x_s] -> [y,y_s] : k = y_s - x_s > 0 }
 585  */
 586 static __isl_give isl_map *equate_parameter_to_length(__isl_take isl_dim *dim,
 587         unsigned param)
 588 {
 589         struct isl_basic_map *bmap;
 590         unsigned d;
 591         unsigned nparam;
 592         int k;
 593
 594         d = isl_dim_size(dim, isl_dim_in);
 595         nparam = isl_dim_size(dim, isl_dim_param);
 596         bmap = isl_basic_map_alloc_dim(dim, 0, 1, 1);
 597         k = isl_basic_map_alloc_equality(bmap);
 598         if (k < 0)
 599                 goto error;
 600         isl_seq_clr(bmap->eq[k], 1 + isl_basic_map_total_dim(bmap));
 601         isl_int_set_si(bmap->eq[k][1 + param], -1);
 602         isl_int_set_si(bmap->eq[k][1 + nparam + d - 1], -1);
 603         isl_int_set_si(bmap->eq[k][1 + nparam + d + d - 1], 1);
 604
 605         k = isl_basic_map_alloc_inequality(bmap);
 606         if (k < 0)
 607                 goto error;
 608         isl_seq_clr(bmap->ineq[k], 1 + isl_basic_map_total_dim(bmap));
 609         isl_int_set_si(bmap->ineq[k][1 + param], 1);
 610         isl_int_set_si(bmap->ineq[k][0], -1);
 611
 612         bmap = isl_basic_map_finalize(bmap);
 613         return isl_map_from_basic_map(bmap);
 614 error:
 615         isl_basic_map_free(bmap);
 616         return NULL;
 617 }
 618
 619 /* Check whether "path" is acyclic, where the last coordinates of domain
 620  * and range of path encode the number of steps taken.
 621  * That is, check whether
 622  *
 623  *      { d | d = y - x and (x,y) in path }
 624  *
 625  * does not contain any element with positive last coordinate (positive length)
 626  * and zero remaining coordinates (cycle).
 627  */
 628 static int is_acyclic(__isl_take isl_map *path)
 629 {
 630         int i;
 631         int acyclic;
 632         unsigned dim;
 633         struct isl_set *delta;
 634
 635         delta = isl_map_deltas(path);
 636         dim = isl_set_dim(delta, isl_dim_set);
 637         for (i = 0; i < dim; ++i) {
 638                 if (i == dim -1)
 639                         delta = isl_set_lower_bound_si(delta, isl_dim_set, i, 1);
 640                 else
 641                         delta = isl_set_fix_si(delta, isl_dim_set, i, 0);
 642         }
 643
 644         acyclic = isl_set_is_empty(delta);
 645         isl_set_free(delta);
 646
 647         return acyclic;
 648 }
 649
 650 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
 651  * and a dimension specification (Z^{n+1} -> Z^{n+1}),
 652  * construct a map that is an overapproximation of the map
 653  * that takes an element from the space D \times Z to another
 654  * element from the same space, such that the first n coordinates of the
 655  * difference between them is a sum of differences between images
 656  * and pre-images in one of the R_i and such that the last coordinate
 657  * is equal to the number of steps taken.
 658  * That is, let
 659  *
 660  *      \Delta_i = { y - x | (x, y) in R_i }
 661  *
 662  * then the constructed map is an overapproximation of
 663  *
 664  *      { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
 665  *                              d = (\sum_i k_i \delta_i, \sum_i k_i) }
 666  *
 667  * The elements of the singleton \Delta_i's are collected as the
 668  * rows of the steps matrix.  For all these \Delta_i's together,
 669  * a single path is constructed.
 670  * For each of the other \Delta_i's, we compute an overapproximation
 671  * of the paths along elements of \Delta_i.
 672  * Since each of these paths performs an addition, composition is
 673  * symmetric and we can simply compose all resulting paths in any order.
 674  */
 675 static __isl_give isl_map *construct_extended_path(__isl_take isl_dim *dim,
 676         __isl_keep isl_map *map, int *project)
 677 {
 678         struct isl_mat *steps = NULL;
 679         struct isl_map *path = NULL;
 680         unsigned d;
 681         int i, j, n;
 682
 683         d = isl_map_dim(map, isl_dim_in);
 684
 685         path = isl_map_identity(isl_dim_domain(isl_dim_copy(dim)));
 686
 687         steps = isl_mat_alloc(map->ctx, map->n, d);
 688         if (!steps)
 689                 goto error;
 690
 691         n = 0;
 692         for (i = 0; i < map->n; ++i) {
 693                 struct isl_basic_set *delta;
 694
 695                 delta = isl_basic_map_deltas(isl_basic_map_copy(map->p[i]));
 696
 697                 for (j = 0; j < d; ++j) {
 698                         int fixed;
 699
 700                         fixed = isl_basic_set_fast_dim_is_fixed(delta, j,
 701                                                             &steps->row[n][j]);
 702                         if (fixed < 0) {
 703                                 isl_basic_set_free(delta);
 704                                 goto error;
 705                         }
 706                         if (!fixed)
 707                                 break;
 708                 }
 709
 710
 711                 if (j < d) {
 712                         path = isl_map_apply_range(path,
 713                                 path_along_delta(isl_dim_copy(dim), delta));
 714                         path = isl_map_coalesce(path);
 715                 } else {
 716                         isl_basic_set_free(delta);
 717                         ++n;
 718                 }
 719         }
 720
 721         if (n > 0) {
 722                 steps->n_row = n;
 723                 path = isl_map_apply_range(path,
 724                                 path_along_steps(isl_dim_copy(dim), steps));
 725         }
 726
 727         if (project && *project) {
 728                 *project = is_acyclic(isl_map_copy(path));
 729                 if (*project < 0)
 730                         goto error;
 731         }
 732
 733         isl_dim_free(dim);
 734         isl_mat_free(steps);
 735         return path;
 736 error:
 737         isl_dim_free(dim);
 738         isl_mat_free(steps);
 739         isl_map_free(path);
 740         return NULL;
 741 }
 742
 743 static int isl_set_overlaps(__isl_keep isl_set *set1, __isl_keep isl_set *set2)
 744 {
 745         isl_set *i;
 746         int no_overlap;
 747
 748         i = isl_set_intersect(isl_set_copy(set1), isl_set_copy(set2));
 749         no_overlap = isl_set_is_empty(i);
 750         isl_set_free(i);
 751
 752         return no_overlap < 0 ? -1 : !no_overlap;
 753 }
 754
 755 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
 756  * and a dimension specification (Z^{n+1} -> Z^{n+1}),
 757  * construct a map that is an overapproximation of the map
 758  * that takes an element from the dom R \times Z to an
 759  * element from ran R \times Z, such that the first n coordinates of the
 760  * difference between them is a sum of differences between images
 761  * and pre-images in one of the R_i and such that the last coordinate
 762  * is equal to the number of steps taken.
 763  * That is, let
 764  *
 765  *      \Delta_i = { y - x | (x, y) in R_i }
 766  *
 767  * then the constructed map is an overapproximation of
 768  *
 769  *      { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
 770  *                              d = (\sum_i k_i \delta_i, \sum_i k_i) and
 771  *                              x in dom R and x + d in ran R and
 772  *                              \sum_i k_i >= 1 }
 773  */
 774 static __isl_give isl_map *construct_component(__isl_take isl_dim *dim,
 775         __isl_keep isl_map *map, int *exact, int project)
 776 {
 777         struct isl_set *domain = NULL;
 778         struct isl_set *range = NULL;
 779         struct isl_map *app = NULL;
 780         struct isl_map *path = NULL;
 781
 782         domain = isl_map_domain(isl_map_copy(map));
 783         domain = isl_set_coalesce(domain);
 784         range = isl_map_range(isl_map_copy(map));
 785         range = isl_set_coalesce(range);
 786         if (!isl_set_overlaps(domain, range)) {
 787                 isl_set_free(domain);
 788                 isl_set_free(range);
 789                 isl_dim_free(dim);
 790
 791                 map = isl_map_copy(map);
 792                 map = isl_map_add(map, isl_dim_in, 1);
 793                 map = isl_map_add(map, isl_dim_out, 1);
 794                 map = set_path_length(map, 1, 1);
 795                 return map;
 796         }
 797         app = isl_map_from_domain_and_range(domain, range);
 798         app = isl_map_add(app, isl_dim_in, 1);
 799         app = isl_map_add(app, isl_dim_out, 1);
 800
 801         path = construct_extended_path(isl_dim_copy(dim), map,
 802                                         exact && *exact ? &project : NULL);
 803         app = isl_map_intersect(app, path);
 804
 805         if (exact && *exact &&
 806             (*exact = check_exactness(isl_map_copy(map), isl_map_copy(app),
 807                                       project)) < 0)
 808                 goto error;
 809
 810         isl_dim_free(dim);
 811         app = set_path_length(app, 0, 1);
 812         return app;
 813 error:
 814         isl_dim_free(dim);
 815         isl_map_free(app);
 816         return NULL;
 817 }
 818
 819 /* Call construct_component and, if "project" is set, project out
 820  * the final coordinates.
 821  */
 822 static __isl_give isl_map *construct_projected_component(
 823         __isl_take isl_dim *dim,
 824         __isl_keep isl_map *map, int *exact, int project)
 825 {
 826         isl_map *app;
 827         unsigned d;
 828
 829         if (!dim)
 830                 return NULL;
 831         d = isl_dim_size(dim, isl_dim_in);
 832
 833         app = construct_component(dim, map, exact, project);
 834         if (project) {
 835                 app = isl_map_project_out(app, isl_dim_in, d - 1, 1);
 836                 app = isl_map_project_out(app, isl_dim_out, d - 1, 1);
 837         }
 838         return app;
 839 }
 840
 841 /* Compute an extended version, i.e., with path lengths, of
 842  * an overapproximation of the transitive closure of "bmap"
 843  * with path lengths greater than or equal to zero and with
 844  * domain and range equal to "dom".
 845  */
 846 static __isl_give isl_map *q_closure(__isl_take isl_dim *dim,
 847         __isl_take isl_set *dom, __isl_keep isl_basic_map *bmap, int *exact)
 848 {
 849         int project = 1;
 850         isl_map *path;
 851         isl_map *map;
 852         isl_map *app;
 853
 854         dom = isl_set_add(dom, isl_dim_set, 1);
 855         app = isl_map_from_domain_and_range(dom, isl_set_copy(dom));
 856         map = isl_map_from_basic_map(isl_basic_map_copy(bmap));
 857         path = construct_extended_path(dim, map, &project);
 858         app = isl_map_intersect(app, path);
 859
 860         if ((*exact = check_exactness(map, isl_map_copy(app), project)) < 0)
 861                 goto error;
 862
 863         return app;
 864 error:
 865         isl_map_free(app);
 866         return NULL;
 867 }
 868
 869 /* Check whether qc has any elements of length at least one
 870  * with domain and/or range outside of dom and ran.
 871  */
 872 static int has_spurious_elements(__isl_keep isl_map *qc,
 873         __isl_keep isl_set *dom, __isl_keep isl_set *ran)
 874 {
 875         isl_set *s;
 876         int subset;
 877         unsigned d;
 878
 879         if (!qc || !dom || !ran)
 880                 return -1;
 881
 882         d = isl_map_dim(qc, isl_dim_in);
 883
 884         qc = isl_map_copy(qc);
 885         qc = set_path_length(qc, 0, 1);
 886         qc = isl_map_project_out(qc, isl_dim_in, d - 1, 1);
 887         qc = isl_map_project_out(qc, isl_dim_out, d - 1, 1);
 888
 889         s = isl_map_domain(isl_map_copy(qc));
 890         subset = isl_set_is_subset(s, dom);
 891         isl_set_free(s);
 892         if (subset < 0)
 893                 goto error;
 894         if (!subset) {
 895                 isl_map_free(qc);
 896                 return 1;
 897         }
 898
 899         s = isl_map_range(qc);
 900         subset = isl_set_is_subset(s, ran);
 901         isl_set_free(s);
 902
 903         return subset < 0 ? -1 : !subset;
 904 error:
 905         isl_map_free(qc);
 906         return -1;
 907 }
 908
 909 #define LEFT    2
 910 #define RIGHT   1
 911
 912 /* For each basic map in "map", except i, check whether it combines
 913  * with the transitive closure that is reflexive on C combines
 914  * to the left and to the right.
 915  *
 916  * In particular, if
 917  *
 918  *      dom map_j \subseteq C
 919  *
 920  * then right[j] is set to 1.  Otherwise, if
 921  *
 922  *      ran map_i \cap dom map_j = \emptyset
 923  *
 924  * then right[j] is set to 0.  Otherwise, composing to the right
 925  * is impossible.
 926  *
 927  * Similar, for composing to the left, we have if
 928  *
 929  *      ran map_j \subseteq C
 930  *
 931  * then left[j] is set to 1.  Otherwise, if
 932  *
 933  *      dom map_i \cap ran map_j = \emptyset
 934  *
 935  * then left[j] is set to 0.  Otherwise, composing to the left
 936  * is impossible.
 937  *
 938  * The return value is or'd with LEFT if composing to the left
 939  * is possible and with RIGHT if composing to the right is possible.
 940  */
 941 static int composability(__isl_keep isl_set *C, int i,
 942         isl_set **dom, isl_set **ran, int *left, int *right,
 943         __isl_keep isl_map *map)
 944 {
 945         int j;
 946         int ok;
 947
 948         ok = LEFT | RIGHT;
 949         for (j = 0; j < map->n && ok; ++j) {
 950                 int overlaps, subset;
 951                 if (j == i)
 952                         continue;
 953
 954                 if (ok & RIGHT) {
 955                         if (!dom[j])
 956                                 dom[j] = isl_set_from_basic_set(
 957                                         isl_basic_map_domain(
 958                                                 isl_basic_map_copy(map->p[j])));
 959                         if (!dom[j])
 960                                 return -1;
 961                         overlaps = isl_set_overlaps(ran[i], dom[j]);
 962                         if (overlaps < 0)
 963                                 return -1;
 964                         if (!overlaps)
 965                                 right[j] = 0;
 966                         else {
 967                                 subset = isl_set_is_subset(dom[j], C);
 968                                 if (subset < 0)
 969                                         return -1;
 970                                 if (subset)
 971                                         right[j] = 1;
 972                                 else
 973                                         ok &= ~RIGHT;
 974                         }
 975                 }
 976
 977                 if (ok & LEFT) {
 978                         if (!ran[j])
 979                                 ran[j] = isl_set_from_basic_set(
 980                                         isl_basic_map_range(
 981                                                 isl_basic_map_copy(map->p[j])));
 982                         if (!ran[j])
 983                                 return -1;
 984                         overlaps = isl_set_overlaps(dom[i], ran[j]);
 985                         if (overlaps < 0)
 986                                 return -1;
 987                         if (!overlaps)
 988                                 left[j] = 0;
 989                         else {
 990                                 subset = isl_set_is_subset(ran[j], C);
 991                                 if (subset < 0)
 992                                         return -1;
 993                                 if (subset)
 994                                         left[j] = 1;
 995                                 else
 996                                         ok &= ~LEFT;
 997                         }
 998                 }
 999         }
1000
1001         return ok;
1002 }
1003
1004 /* Return a map that is a union of the basic maps in "map", except i,
1005  * composed to left and right with qc based on the entries of "left"
1006  * and "right".
1007  */
1008 static __isl_give isl_map *compose(__isl_keep isl_map *map, int i,
1009         __isl_take isl_map *qc, int *left, int *right)
1010 {
1011         int j;
1012         isl_map *comp;
1013
1014         comp = isl_map_empty(isl_map_get_dim(map));
1015         for (j = 0; j < map->n; ++j) {
1016                 isl_map *map_j;
1017
1018                 if (j == i)
1019                         continue;
1020
1021                 map_j = isl_map_from_basic_map(isl_basic_map_copy(map->p[j]));
1022                 if (left && left[j])
1023                         map_j = isl_map_apply_range(map_j, isl_map_copy(qc));
1024                 if (right && right[j])
1025                         map_j = isl_map_apply_range(isl_map_copy(qc), map_j);
1026                 comp = isl_map_union(comp, map_j);
1027         }
1028
1029         comp = isl_map_compute_divs(comp);
1030         comp = isl_map_coalesce(comp);
1031
1032         isl_map_free(qc);
1033
1034         return comp;
1035 }
1036
1037 /* Compute the transitive closure of "map" incrementally by
1038  * computing
1039  *
1040  *      map_i^+ \cup qc^+
1041  *
1042  * or
1043  *
1044  *      map_i^+ \cup ((id \cup map_i^) \circ qc^+)
1045  *
1046  * or
1047  *
1048  *      map_i^+ \cup (qc^+ \circ (id \cup map_i^))
1049  *
1050  * depending on whether left or right are NULL.
1051  */
1052 static __isl_give isl_map *compute_incremental(
1053         __isl_take isl_dim *dim, __isl_keep isl_map *map,
1054         int i, __isl_take isl_map *qc, int *left, int *right, int *exact)
1055 {
1056         isl_map *map_i;
1057         isl_map *tc;
1058         isl_map *rtc = NULL;
1059
1060         if (!map)
1061                 goto error;
1062         isl_assert(map->ctx, left || right, goto error);
1063
1064         map_i = isl_map_from_basic_map(isl_basic_map_copy(map->p[i]));
1065         tc = construct_projected_component(isl_dim_copy(dim), map_i,
1066                                                 exact, 1);
1067         isl_map_free(map_i);
1068
1069         if (*exact)
1070                 qc = isl_map_transitive_closure(qc, exact);
1071
1072         if (!*exact) {
1073                 isl_dim_free(dim);
1074                 isl_map_free(tc);
1075                 isl_map_free(qc);
1076                 return isl_map_universe(isl_map_get_dim(map));
1077         }
1078
1079         if (!left || !right)
1080                 rtc = isl_map_union(isl_map_copy(tc),
1081                         isl_map_identity(isl_dim_domain(isl_map_get_dim(tc))));
1082         if (!right)
1083                 qc = isl_map_apply_range(rtc, qc);
1084         if (!left)
1085                 qc = isl_map_apply_range(qc, rtc);
1086         qc = isl_map_union(tc, qc);
1087
1088         isl_dim_free(dim);
1089
1090         return qc;
1091 error:
1092         isl_dim_free(dim);
1093         isl_map_free(qc);
1094         return NULL;
1095 }
1096
1097 /* Given a map "map", try to find a basic map such that
1098  * map^+ can be computed as
1099  *
1100  * map^+ = map_i^+ \cup
1101  *    \bigcup_j ((map_i^+ \cup Id_C)^+ \circ map_j \circ (map_i^+ \cup Id_C))^+
1102  *
1103  * with C the simple hull of the domain and range of the input map.
1104  * map_i^ \cup Id_C is computed by allowing the path lengths to be zero
1105  * and by intersecting domain and range with C.
1106  * Of course, we need to check that this is actually equal to map_i^ \cup Id_C.
1107  * Also, we only use the incremental computation if all the transitive
1108  * closures are exact and if the number of basic maps in the union,
1109  * after computing the integer divisions, is smaller than the number
1110  * of basic maps in the input map.
1111  */
1112 static int incemental_on_entire_domain(__isl_keep isl_dim *dim,
1113         __isl_keep isl_map *map,
1114         isl_set **dom, isl_set **ran, int *left, int *right,
1115         __isl_give isl_map **res)
1116 {
1117         int i;
1118         isl_set *C;
1119         unsigned d;
1120
1121         *res = NULL;
1122
1123         C = isl_set_union(isl_map_domain(isl_map_copy(map)),
1124                           isl_map_range(isl_map_copy(map)));
1125         C = isl_set_from_basic_set(isl_set_simple_hull(C));
1126         if (!C)
1127                 return -1;
1128         if (C->n != 1) {
1129                 isl_set_free(C);
1130                 return 0;
1131         }
1132
1133         d = isl_map_dim(map, isl_dim_in);
1134
1135         for (i = 0; i < map->n; ++i) {
1136                 isl_map *qc;
1137                 int exact_i, spurious;
1138                 int j;
1139                 dom[i] = isl_set_from_basic_set(isl_basic_map_domain(
1140                                         isl_basic_map_copy(map->p[i])));
1141                 ran[i] = isl_set_from_basic_set(isl_basic_map_range(
1142                                         isl_basic_map_copy(map->p[i])));
1143                 qc = q_closure(isl_dim_copy(dim), isl_set_copy(C),
1144                                 map->p[i], &exact_i);
1145                 if (!qc)
1146                         goto error;
1147                 if (!exact_i) {
1148                         isl_map_free(qc);
1149                         continue;
1150                 }
1151                 spurious = has_spurious_elements(qc, dom[i], ran[i]);
1152                 if (spurious) {
1153                         isl_map_free(qc);
1154                         if (spurious < 0)
1155                                 goto error;
1156                         continue;
1157                 }
1158                 qc = isl_map_project_out(qc, isl_dim_in, d, 1);
1159                 qc = isl_map_project_out(qc, isl_dim_out, d, 1);
1160                 qc = isl_map_compute_divs(qc);
1161                 for (j = 0; j < map->n; ++j)
1162                         left[j] = right[j] = 1;
1163                 qc = compose(map, i, qc, left, right);
1164                 if (!qc)
1165                         goto error;
1166                 if (qc->n >= map->n) {
1167                         isl_map_free(qc);
1168                         continue;
1169                 }
1170                 *res = compute_incremental(isl_dim_copy(dim), map, i, qc,
1171                                 left, right, &exact_i);
1172                 if (!*res)
1173                         goto error;
1174                 if (exact_i)
1175                         break;
1176                 isl_map_free(*res);
1177                 *res = NULL;
1178         }
1179
1180         isl_set_free(C);
1181
1182         return *res != NULL;
1183 error:
1184         isl_set_free(C);
1185         return -1;
1186 }
1187
1188 /* Try and compute the transitive closure of "map" as
1189  *
1190  * map^+ = map_i^+ \cup
1191  *    \bigcup_j ((map_i^+ \cup Id_C)^+ \circ map_j \circ (map_i^+ \cup Id_C))^+
1192  *
1193  * with C either the simple hull of the domain and range of the entire
1194  * map or the simple hull of domain and range of map_i.
1195  */
1196 static __isl_give isl_map *incremental_closure(__isl_take isl_dim *dim,
1197         __isl_keep isl_map *map, int *exact, int project)
1198 {
1199         int i;
1200         isl_set **dom = NULL;
1201         isl_set **ran = NULL;
1202         int *left = NULL;
1203         int *right = NULL;
1204         isl_set *C;
1205         unsigned d;
1206         isl_map *res = NULL;
1207
1208         if (!project)
1209                 return construct_projected_component(dim, map, exact, project);
1210
1211         if (!map)
1212                 goto error;
1213         if (map->n <= 1)
1214                 return construct_projected_component(dim, map, exact, project);
1215
1216         d = isl_map_dim(map, isl_dim_in);
1217
1218         dom = isl_calloc_array(map->ctx, isl_set *, map->n);
1219         ran = isl_calloc_array(map->ctx, isl_set *, map->n);
1220         left = isl_calloc_array(map->ctx, int, map->n);
1221         right = isl_calloc_array(map->ctx, int, map->n);
1222         if (!ran || !dom || !left || !right)
1223                 goto error;
1224
1225         if (incemental_on_entire_domain(dim, map, dom, ran, left, right, &res) < 0)
1226                 goto error;
1227
1228         for (i = 0; !res && i < map->n; ++i) {
1229                 isl_map *qc;
1230                 int exact_i, spurious, comp;
1231                 if (!dom[i])
1232                         dom[i] = isl_set_from_basic_set(
1233                                         isl_basic_map_domain(
1234                                                 isl_basic_map_copy(map->p[i])));
1235                 if (!dom[i])
1236                         goto error;
1237                 if (!ran[i])
1238                         ran[i] = isl_set_from_basic_set(
1239                                         isl_basic_map_range(
1240                                                 isl_basic_map_copy(map->p[i])));
1241                 if (!ran[i])
1242                         goto error;
1243                 C = isl_set_union(isl_set_copy(dom[i]),
1244                                       isl_set_copy(ran[i]));
1245                 C = isl_set_from_basic_set(isl_set_simple_hull(C));
1246                 if (!C)
1247                         goto error;
1248                 if (C->n != 1) {
1249                         isl_set_free(C);
1250                         continue;
1251                 }
1252                 comp = composability(C, i, dom, ran, left, right, map);
1253                 if (!comp || comp < 0) {
1254                         isl_set_free(C);
1255                         if (comp < 0)
1256                                 goto error;
1257                         continue;
1258                 }
1259                 qc = q_closure(isl_dim_copy(dim), C, map->p[i], &exact_i);
1260                 if (!qc)
1261                         goto error;
1262                 if (!exact_i) {
1263                         isl_map_free(qc);
1264                         continue;
1265                 }
1266                 spurious = has_spurious_elements(qc, dom[i], ran[i]);
1267                 if (spurious) {
1268                         isl_map_free(qc);
1269                         if (spurious < 0)
1270                                 goto error;
1271                         continue;
1272                 }
1273                 qc = isl_map_project_out(qc, isl_dim_in, d, 1);
1274                 qc = isl_map_project_out(qc, isl_dim_out, d, 1);
1275                 qc = isl_map_compute_divs(qc);
1276                 qc = compose(map, i, qc, (comp & LEFT) ? left : NULL,
1277                                 (comp & RIGHT) ? right : NULL);
1278                 if (!qc)
1279                         goto error;
1280                 if (qc->n >= map->n) {
1281                         isl_map_free(qc);
1282                         continue;
1283                 }
1284                 res = compute_incremental(isl_dim_copy(dim), map, i, qc,
1285                                 (comp & LEFT) ? left : NULL,
1286                                 (comp & RIGHT) ? right : NULL, &exact_i);
1287                 if (!res)
1288                         goto error;
1289                 if (exact_i)
1290                         break;
1291                 isl_map_free(res);
1292                 res = NULL;
1293         }
1294
1295         for (i = 0; i < map->n; ++i) {
1296                 isl_set_free(dom[i]);
1297                 isl_set_free(ran[i]);
1298         }
1299         free(dom);
1300         free(ran);
1301         free(left);
1302         free(right);
1303
1304         if (res) {
1305                 isl_dim_free(dim);
1306                 return res;
1307         }
1308
1309         return construct_projected_component(dim, map, exact, project);
1310 error:
1311         if (dom)
1312                 for (i = 0; i < map->n; ++i)
1313                         isl_set_free(dom[i]);
1314         free(dom);
1315         if (ran)
1316                 for (i = 0; i < map->n; ++i)
1317                         isl_set_free(ran[i]);
1318         free(ran);
1319         free(left);
1320         free(right);
1321         isl_dim_free(dim);
1322         return NULL;
1323 }
1324
1325 /* Given an array of sets "set", add "dom" at position "pos"
1326  * and search for elements at earlier positions that overlap with "dom".
1327  * If any can be found, then merge all of them, together with "dom", into
1328  * a single set and assign the union to the first in the array,
1329  * which becomes the new group leader for all groups involved in the merge.
1330  * During the search, we only consider group leaders, i.e., those with
1331  * group[i] = i, as the other sets have already been combined
1332  * with one of the group leaders.
1333  */
1334 static int merge(isl_set **set, int *group, __isl_take isl_set *dom, int pos)
1335 {
1336         int i;
1337
1338         group[pos] = pos;
1339         set[pos] = isl_set_copy(dom);
1340
1341         for (i = pos - 1; i >= 0; --i) {
1342                 int o;
1343
1344                 if (group[i] != i)
1345                         continue;
1346
1347                 o = isl_set_overlaps(set[i], dom);
1348                 if (o < 0)
1349                         goto error;
1350                 if (!o)
1351                         continue;
1352
1353                 set[i] = isl_set_union(set[i], set[group[pos]]);
1354                 if (!set[i])
1355                         goto error;
1356                 set[group[pos]] = NULL;
1357                 group[group[pos]] = i;
1358                 group[pos] = i;
1359         }
1360
1361         isl_set_free(dom);
1362         return 0;
1363 error:
1364         isl_set_free(dom);
1365         return -1;
1366 }
1367
1368 /* Given a partition of the domains and ranges of the basic maps in "map",
1369  * apply the Floyd-Warshall algorithm with the elements in the partition
1370  * as vertices.
1371  *
1372  * In particular, there are "n" elements in the partition and "group" is
1373  * an array of length 2 * map->n with entries in [0,n-1].
1374  *
1375  * We first construct a matrix of relations based on the partition information,
1376  * apply Floyd-Warshall on this matrix of relations and then take the
1377  * union of all entries in the matrix as the final result.
1378  *
1379  * The algorithm iterates over all vertices.  In each step, the whole
1380  * matrix is updated to include all paths that go to the current vertex,
1381  * possibly stay there a while (including passing through earlier vertices)
1382  * and then come back.  At the start of each iteration, the diagonal
1383  * element corresponding to the current vertex is replaced by its
1384  * transitive closure to account for all indirect paths that stay
1385  * in the current vertex.
1386  */
1387 static __isl_give isl_map *floyd_warshall_with_groups(__isl_take isl_dim *dim,
1388         __isl_keep isl_map *map, int *exact, int project, int *group, int n)
1389 {
1390         int i, j, k;
1391         int r, p, q;
1392         isl_map ***grid = NULL;
1393         isl_map *app;
1394
1395         if (!map)
1396                 goto error;
1397
1398         if (n == 1) {
1399                 free(group);
1400                 return incremental_closure(dim, map, exact, project);
1401         }
1402
1403         grid = isl_calloc_array(map->ctx, isl_map **, n);
1404         if (!grid)
1405                 goto error;
1406         for (i = 0; i < n; ++i) {
1407                 grid[i] = isl_calloc_array(map->ctx, isl_map *, n);
1408                 if (!grid[i])
1409                         goto error;
1410                 for (j = 0; j < n; ++j)
1411                         grid[i][j] = isl_map_empty(isl_map_get_dim(map));
1412         }
1413
1414         for (k = 0; k < map->n; ++k) {
1415                 i = group[2 * k];
1416                 j = group[2 * k + 1];
1417                 grid[i][j] = isl_map_union(grid[i][j],
1418                                 isl_map_from_basic_map(
1419                                         isl_basic_map_copy(map->p[k])));
1420         }
1421
1422         for (r = 0; r < n; ++r) {
1423                 int r_exact;
1424                 grid[r][r] = isl_map_transitive_closure(grid[r][r],
1425                                 (exact && *exact) ? &r_exact : NULL);
1426                 if (exact && *exact && !r_exact)
1427                         *exact = 0;
1428
1429                 for (p = 0; p < n; ++p)
1430                         for (q = 0; q < n; ++q) {
1431                                 isl_map *loop;
1432                                 if (p == r && q == r)
1433                                         continue;
1434                                 loop = isl_map_apply_range(
1435                                                 isl_map_copy(grid[p][r]),
1436                                                 isl_map_copy(grid[r][q]));
1437                                 grid[p][q] = isl_map_union(grid[p][q], loop);
1438                                 loop = isl_map_apply_range(
1439                                                 isl_map_copy(grid[p][r]),
1440                                         isl_map_apply_range(
1441                                                 isl_map_copy(grid[r][r]),
1442                                                 isl_map_copy(grid[r][q])));
1443                                 grid[p][q] = isl_map_union(grid[p][q], loop);
1444                                 grid[p][q] = isl_map_coalesce(grid[p][q]);
1445                         }
1446         }
1447
1448         app = isl_map_empty(isl_map_get_dim(map));
1449
1450         for (i = 0; i < n; ++i) {
1451                 for (j = 0; j < n; ++j)
1452                         app = isl_map_union(app, grid[i][j]);
1453                 free(grid[i]);
1454         }
1455         free(grid);
1456
1457         free(group);
1458         isl_dim_free(dim);
1459
1460         return app;
1461 error:
1462         if (grid)
1463                 for (i = 0; i < n; ++i) {
1464                         if (!grid[i])
1465                                 continue;
1466                         for (j = 0; j < n; ++j)
1467                                 isl_map_free(grid[i][j]);
1468                         free(grid[i]);
1469                 }
1470         free(grid);
1471         free(group);
1472         isl_dim_free(dim);
1473         return NULL;
1474 }
1475
1476 /* Check if the domains and ranges of the basic maps in "map" can
1477  * be partitioned, and if so, apply Floyd-Warshall on the elements
1478  * of the partition.  Note that we can only apply this algorithm
1479  * if we want to compute the transitive closure, i.e., when "project"
1480  * is set.  If we want to compute the power, we need to keep track
1481  * of the lengths and the recursive calls inside the Floyd-Warshall
1482  * would result in non-linear lengths.
1483  *
1484  * To find the partition, we simply consider all of the domains
1485  * and ranges in turn and combine those that overlap.
1486  * "set" contains the partition elements and "group" indicates
1487  * to which partition element a given domain or range belongs.
1488  * The domain of basic map i corresponds to element 2 * i in these arrays,
1489  * while the domain corresponds to element 2 * i + 1.
1490  * During the construction group[k] is either equal to k,
1491  * in which case set[k] contains the union of all the domains and
1492  * ranges in the corresponding group, or is equal to some l < k,
1493  * with l another domain or range in the same group.
1494  */
1495 static __isl_give isl_map *floyd_warshall(__isl_take isl_dim *dim,
1496         __isl_keep isl_map *map, int *exact, int project)
1497 {
1498         int i;
1499         isl_set **set = NULL;
1500         int *group = NULL;
1501         int n;
1502
1503         if (!map)
1504                 goto error;
1505         if (!project || map->n <= 1)
1506                 return incremental_closure(dim, map, exact, project);
1507
1508         set = isl_calloc_array(map->ctx, isl_set *, 2 * map->n);
1509         group = isl_alloc_array(map->ctx, int, 2 * map->n);
1510
1511         if (!set || !group)
1512                 goto error;
1513
1514         for (i = 0; i < map->n; ++i) {
1515                 isl_set *dom;
1516                 dom = isl_set_from_basic_set(isl_basic_map_domain(
1517                                 isl_basic_map_copy(map->p[i])));
1518                 if (merge(set, group, dom, 2 * i) < 0)
1519                         goto error;
1520                 dom = isl_set_from_basic_set(isl_basic_map_range(
1521                                 isl_basic_map_copy(map->p[i])));
1522                 if (merge(set, group, dom, 2 * i + 1) < 0)
1523                         goto error;
1524         }
1525
1526         n = 0;
1527         for (i = 0; i < 2 * map->n; ++i)
1528                 if (group[i] == i)
1529                         group[i] = n++;
1530                 else
1531                         group[i] = group[group[i]];
1532
1533         for (i = 0; i < 2 * map->n; ++i)
1534                 isl_set_free(set[i]);
1535
1536         free(set);
1537
1538         return floyd_warshall_with_groups(dim, map, exact, project, group, n);
1539 error:
1540         for (i = 0; i < 2 * map->n; ++i)
1541                 isl_set_free(set[i]);
1542         free(set);
1543         free(group);
1544         isl_dim_free(dim);
1545         return NULL;
1546 }
1547
1548 /* Structure for representing the nodes in the graph being traversed
1549  * using Tarjan's algorithm.
1550  * index represents the order in which nodes are visited.
1551  * min_index is the index of the root of a (sub)component.
1552  * on_stack indicates whether the node is currently on the stack.
1553  */
1554 struct basic_map_sort_node {
1555         int index;
1556         int min_index;
1557         int on_stack;
1558 };
1559 /* Structure for representing the graph being traversed
1560  * using Tarjan's algorithm.
1561  * len is the number of nodes
1562  * node is an array of nodes
1563  * stack contains the nodes on the path from the root to the current node
1564  * sp is the stack pointer
1565  * index is the index of the last node visited
1566  * order contains the elements of the components separated by -1
1567  * op represents the current position in order
1568  */
1569 struct basic_map_sort {
1570         int len;
1571         struct basic_map_sort_node *node;
1572         int *stack;
1573         int sp;
1574         int index;
1575         int *order;
1576         int op;
1577 };
1578
1579 static void basic_map_sort_free(struct basic_map_sort *s)
1580 {
1581         if (!s)
1582                 return;
1583         free(s->node);
1584         free(s->stack);
1585         free(s->order);
1586         free(s);
1587 }
1588
1589 static struct basic_map_sort *basic_map_sort_alloc(struct isl_ctx *ctx, int len)
1590 {
1591         struct basic_map_sort *s;
1592         int i;
1593
1594         s = isl_calloc_type(ctx, struct basic_map_sort);
1595         if (!s)
1596                 return NULL;
1597         s->len = len;
1598         s->node = isl_alloc_array(ctx, struct basic_map_sort_node, len);
1599         if (!s->node)
1600                 goto error;
1601         for (i = 0; i < len; ++i)
1602                 s->node[i].index = -1;
1603         s->stack = isl_alloc_array(ctx, int, len);
1604         if (!s->stack)
1605                 goto error;
1606         s->order = isl_alloc_array(ctx, int, 2 * len);
1607         if (!s->order)
1608                 goto error;
1609
1610         s->sp = 0;
1611         s->index = 0;
1612         s->op = 0;
1613
1614         return s;
1615 error:
1616         basic_map_sort_free(s);
1617         return NULL;
1618 }
1619
1620 /* Check whether in the computation of the transitive closure
1621  * "bmap1" (R_1) should follow (or be part of the same component as)
1622  * "bmap2" (R_2).
1623  *
1624  * That is check whether
1625  *
1626  *      R_1 \circ R_2
1627  *
1628  * is a subset of
1629  *
1630  *      R_2 \circ R_1
1631  *
1632  * If so, then there is no reason for R_1 to immediately follow R_2
1633  * in any path.
1634  */
1635 static int basic_map_follows(__isl_keep isl_basic_map *bmap1,
1636         __isl_keep isl_basic_map *bmap2)
1637 {
1638         struct isl_map *map12 = NULL;
1639         struct isl_map *map21 = NULL;
1640         int subset;
1641
1642         map21 = isl_map_from_basic_map(
1643                         isl_basic_map_apply_range(
1644                                 isl_basic_map_copy(bmap2),
1645                                 isl_basic_map_copy(bmap1)));
1646         subset = isl_map_is_empty(map21);
1647         if (subset < 0)
1648                 goto error;
1649         if (subset) {
1650                 isl_map_free(map21);
1651                 return 0;
1652         }
1653
1654         map12 = isl_map_from_basic_map(
1655                         isl_basic_map_apply_range(
1656                                 isl_basic_map_copy(bmap1),
1657                                 isl_basic_map_copy(bmap2)));
1658
1659         subset = isl_map_is_subset(map21, map12);
1660
1661         isl_map_free(map12);
1662         isl_map_free(map21);
1663
1664         return subset < 0 ? -1 : !subset;
1665 error:
1666         isl_map_free(map21);
1667         return -1;
1668 }
1669
1670 /* Perform Tarjan's algorithm for computing the strongly connected components
1671  * in the graph with the disjuncts of "map" as vertices and with an
1672  * edge between any pair of disjuncts such that the first has
1673  * to be applied after the second.
1674  */
1675 static int power_components_tarjan(struct basic_map_sort *s,
1676         __isl_keep isl_map *map, int i)
1677 {
1678         int j;
1679
1680         s->node[i].index = s->index;
1681         s->node[i].min_index = s->index;
1682         s->node[i].on_stack = 1;
1683         s->index++;
1684         s->stack[s->sp++] = i;
1685
1686         for (j = s->len - 1; j >= 0; --j) {
1687                 int f;
1688
1689                 if (j == i)
1690                         continue;
1691                 if (s->node[j].index >= 0 &&
1692                         (!s->node[j].on_stack ||
1693                          s->node[j].index > s->node[i].min_index))
1694                         continue;
1695
1696                 f = basic_map_follows(map->p[i], map->p[j]);
1697                 if (f < 0)
1698                         return -1;
1699                 if (!f)
1700                         continue;
1701
1702                 if (s->node[j].index < 0) {
1703                         power_components_tarjan(s, map, j);
1704                         if (s->node[j].min_index < s->node[i].min_index)
1705                                 s->node[i].min_index = s->node[j].min_index;
1706                 } else if (s->node[j].index < s->node[i].min_index)
1707                         s->node[i].min_index = s->node[j].index;
1708         }
1709
1710         if (s->node[i].index != s->node[i].min_index)
1711                 return 0;
1712
1713         do {
1714                 j = s->stack[--s->sp];
1715                 s->node[j].on_stack = 0;
1716                 s->order[s->op++] = j;
1717         } while (j != i);
1718         s->order[s->op++] = -1;
1719
1720         return 0;
1721 }
1722
1723 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
1724  * and a dimension specification (Z^{n+1} -> Z^{n+1}),
1725  * construct a map that is an overapproximation of the map
1726  * that takes an element from the dom R \times Z to an
1727  * element from ran R \times Z, such that the first n coordinates of the
1728  * difference between them is a sum of differences between images
1729  * and pre-images in one of the R_i and such that the last coordinate
1730  * is equal to the number of steps taken.
1731  * If "project" is set, then these final coordinates are not included,
1732  * i.e., a relation of type Z^n -> Z^n is returned.
1733  * That is, let
1734  *
1735  *      \Delta_i = { y - x | (x, y) in R_i }
1736  *
1737  * then the constructed map is an overapproximation of
1738  *
1739  *      { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
1740  *                              d = (\sum_i k_i \delta_i, \sum_i k_i) and
1741  *                              x in dom R and x + d in ran R }
1742  *
1743  * or
1744  *
1745  *      { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
1746  *                              d = (\sum_i k_i \delta_i) and
1747  *                              x in dom R and x + d in ran R }
1748  *
1749  * if "project" is set.
1750  *
1751  * We first split the map into strongly connected components, perform
1752  * the above on each component and then join the results in the correct
1753  * order, at each join also taking in the union of both arguments
1754  * to allow for paths that do not go through one of the two arguments.
1755  */
1756 static __isl_give isl_map *construct_power_components(__isl_take isl_dim *dim,
1757         __isl_keep isl_map *map, int *exact, int project)
1758 {
1759         int i, n;
1760         struct isl_map *path = NULL;
1761         struct basic_map_sort *s = NULL;
1762
1763         if (!map)
1764                 goto error;
1765         if (map->n <= 1)
1766                 return floyd_warshall(dim, map, exact, project);
1767
1768         s = basic_map_sort_alloc(map->ctx, map->n);
1769         if (!s)
1770                 goto error;
1771         for (i = map->n - 1; i >= 0; --i) {
1772                 if (s->node[i].index >= 0)
1773                         continue;
1774                 if (power_components_tarjan(s, map, i) < 0)
1775                         goto error;
1776         }
1777
1778         i = 0;
1779         n = map->n;
1780         if (project)
1781                 path = isl_map_empty(isl_map_get_dim(map));
1782         else
1783                 path = isl_map_empty(isl_dim_copy(dim));
1784         while (n) {
1785                 struct isl_map *comp;
1786                 isl_map *path_comp, *path_comb;
1787                 comp = isl_map_alloc_dim(isl_map_get_dim(map), n, 0);
1788                 while (s->order[i] != -1) {
1789                         comp = isl_map_add_basic_map(comp,
1790                                     isl_basic_map_copy(map->p[s->order[i]]));
1791                         --n;
1792                         ++i;
1793                 }
1794                 path_comp = floyd_warshall(isl_dim_copy(dim),
1795                                                 comp, exact, project);
1796                 path_comb = isl_map_apply_range(isl_map_copy(path),
1797                                                 isl_map_copy(path_comp));
1798                 path = isl_map_union(path, path_comp);
1799                 path = isl_map_union(path, path_comb);
1800                 isl_map_free(comp);
1801                 ++i;
1802         }
1803
1804         basic_map_sort_free(s);
1805         isl_dim_free(dim);
1806
1807         return path;
1808 error:
1809         basic_map_sort_free(s);
1810         isl_dim_free(dim);
1811         return NULL;
1812 }
1813
1814 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D,
1815  * construct a map that is an overapproximation of the map
1816  * that takes an element from the space D to another
1817  * element from the same space, such that the difference between
1818  * them is a strictly positive sum of differences between images
1819  * and pre-images in one of the R_i.
1820  * The number of differences in the sum is equated to parameter "param".
1821  * That is, let
1822  *
1823  *      \Delta_i = { y - x | (x, y) in R_i }
1824  *
1825  * then the constructed map is an overapproximation of
1826  *
1827  *      { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
1828  *                              d = \sum_i k_i \delta_i and k = \sum_i k_i > 0 }
1829  * or
1830  *
1831  *      { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
1832  *                              d = \sum_i k_i \delta_i and \sum_i k_i > 0 }
1833  *
1834  * if "project" is set.
1835  *
1836  * If "project" is not set, then
1837  * we first construct an extended mapping with an extra coordinate
1838  * that indicates the number of steps taken.  In particular,
1839  * the difference in the last coordinate is equal to the number
1840  * of steps taken to move from a domain element to the corresponding
1841  * image element(s).
1842  * In the final step, this difference is equated to the parameter "param"
1843  * and made positive.  The extra coordinates are subsequently projected out.
1844  */
1845 static __isl_give isl_map *construct_power(__isl_keep isl_map *map,
1846         unsigned param, int *exact, int project)
1847 {
1848         struct isl_map *app = NULL;
1849         struct isl_map *diff;
1850         struct isl_dim *dim = NULL;
1851         unsigned d;
1852
1853         if (!map)
1854                 return NULL;
1855
1856         dim = isl_map_get_dim(map);
1857
1858         d = isl_dim_size(dim, isl_dim_in);
1859         dim = isl_dim_add(dim, isl_dim_in, 1);
1860         dim = isl_dim_add(dim, isl_dim_out, 1);
1861
1862         app = construct_power_components(isl_dim_copy(dim), map,
1863                                         exact, project);
1864
1865         if (project) {
1866                 isl_dim_free(dim);
1867         } else {
1868                 diff = equate_parameter_to_length(dim, param);
1869                 app = isl_map_intersect(app, diff);
1870                 app = isl_map_project_out(app, isl_dim_in, d, 1);
1871                 app = isl_map_project_out(app, isl_dim_out, d, 1);
1872         }
1873
1874         return app;
1875 }
1876
1877 /* Compute the positive powers of "map", or an overapproximation.
1878  * The power is given by parameter "param".  If the result is exact,
1879  * then *exact is set to 1.
1880  *
1881  * If project is set, then we are actually interested in the transitive
1882  * closure, so we can use a more relaxed exactness check.
1883  * The lengths of the paths are also projected out instead of being
1884  * equated to "param" (which is then ignored in this case).
1885  */
1886 static __isl_give isl_map *map_power(__isl_take isl_map *map, unsigned param,
1887         int *exact, int project)
1888 {
1889         struct isl_map *app = NULL;
1890
1891         if (exact)
1892                 *exact = 1;
1893
1894         if (!map)
1895                 return NULL;
1896
1897         if (isl_map_fast_is_empty(map))
1898                 return map;
1899
1900         isl_assert(map->ctx, project || param < isl_map_dim(map, isl_dim_param),
1901                 goto error);
1902         isl_assert(map->ctx,
1903                 isl_map_dim(map, isl_dim_in) == isl_map_dim(map, isl_dim_out),
1904                 goto error);
1905
1906         app = construct_power(map, param, exact, project);
1907
1908         isl_map_free(map);
1909         return app;
1910 error:
1911         isl_map_free(map);
1912         isl_map_free(app);
1913         return NULL;
1914 }
1915
1916 /* Compute the positive powers of "map", or an overapproximation.
1917  * The power is given by parameter "param".  If the result is exact,
1918  * then *exact is set to 1.
1919  */
1920 __isl_give isl_map *isl_map_power(__isl_take isl_map *map, unsigned param,
1921         int *exact)
1922 {
1923         map = isl_map_compute_divs(map);
1924         map = isl_map_coalesce(map);
1925         return map_power(map, param, exact, 0);
1926 }
1927
1928 /* Check whether equality i of bset is a pure stride constraint
1929  * on a single dimensions, i.e., of the form
1930  *
1931  *      v = k e
1932  *
1933  * with k a constant and e an existentially quantified variable.
1934  */
1935 static int is_eq_stride(__isl_keep isl_basic_set *bset, int i)
1936 {
1937         int k;
1938         unsigned nparam;
1939         unsigned d;
1940         unsigned n_div;
1941         int pos1;
1942         int pos2;
1943
1944         if (!bset)
1945                 return -1;
1946
1947         if (!isl_int_is_zero(bset->eq[i][0]))
1948                 return 0;
1949
1950         nparam = isl_basic_set_dim(bset, isl_dim_param);
1951         d = isl_basic_set_dim(bset, isl_dim_set);
1952         n_div = isl_basic_set_dim(bset, isl_dim_div);
1953
1954         if (isl_seq_first_non_zero(bset->eq[i] + 1, nparam) != -1)
1955                 return 0;
1956         pos1 = isl_seq_first_non_zero(bset->eq[i] + 1 + nparam, d);
1957         if (pos1 == -1)
1958                 return 0;
1959         if (isl_seq_first_non_zero(bset->eq[i] + 1 + nparam + pos1 + 1,
1960                                         d - pos1 - 1) != -1)
1961                 return 0;
1962
1963         pos2 = isl_seq_first_non_zero(bset->eq[i] + 1 + nparam + d, n_div);
1964         if (pos2 == -1)
1965                 return 0;
1966         if (isl_seq_first_non_zero(bset->eq[i] + 1 + nparam + d  + pos2 + 1,
1967                                    n_div - pos2 - 1) != -1)
1968                 return 0;
1969         if (!isl_int_is_one(bset->eq[i][1 + nparam + pos1]) &&
1970             !isl_int_is_negone(bset->eq[i][1 + nparam + pos1]))
1971                 return 0;
1972
1973         return 1;
1974 }
1975
1976 /* Given a map, compute the smallest superset of this map that is of the form
1977  *
1978  *      { i -> j : L <= j - i <= U and exists a_p: j_p - i_p = M_p a_p }
1979  *
1980  * (where p ranges over the (non-parametric) dimensions),
1981  * compute the transitive closure of this map, i.e.,
1982  *
1983  *      { i -> j : exists k > 0:
1984  *              k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p }
1985  *
1986  * and intersect domain and range of this transitive closure with
1987  * the given domain and range.
1988  *
1989  * If with_id is set, then try to include as much of the identity mapping
1990  * as possible, by computing
1991  *
1992  *      { i -> j : exists k >= 0:
1993  *              k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p }
1994  *
1995  * instead (i.e., allow k = 0).
1996  *
1997  * In practice, we compute the difference set
1998  *
1999  *      delta  = { j - i | i -> j in map },
2000  *
2001  * look for stride constraint on the individual dimensions and compute
2002  * (constant) lower and upper bounds for each individual dimension,
2003  * adding a constraint for each bound not equal to infinity.
2004  */
2005 static __isl_give isl_map *box_closure_on_domain(__isl_take isl_map *map,
2006         __isl_take isl_set *dom, __isl_take isl_set *ran, int with_id)
2007 {
2008         int i;
2009         int k;
2010         unsigned d;
2011         unsigned nparam;
2012         unsigned total;
2013         isl_dim *dim;
2014         isl_set *delta;
2015         isl_map *app = NULL;
2016         isl_basic_set *aff = NULL;
2017         isl_basic_map *bmap = NULL;
2018         isl_vec *obj = NULL;
2019         isl_int opt;
2020
2021         isl_int_init(opt);
2022
2023         delta = isl_map_deltas(isl_map_copy(map));
2024
2025         aff = isl_set_affine_hull(isl_set_copy(delta));
2026         if (!aff)
2027                 goto error;
2028         dim = isl_map_get_dim(map);
2029         d = isl_dim_size(dim, isl_dim_in);
2030         nparam = isl_dim_size(dim, isl_dim_param);
2031         total = isl_dim_total(dim);
2032         bmap = isl_basic_map_alloc_dim(dim,
2033                                         aff->n_div + 1, aff->n_div, 2 * d + 1);
2034         for (i = 0; i < aff->n_div + 1; ++i) {
2035                 k = isl_basic_map_alloc_div(bmap);
2036                 if (k < 0)
2037                         goto error;
2038                 isl_int_set_si(bmap->div[k][0], 0);
2039         }
2040         for (i = 0; i < aff->n_eq; ++i) {
2041                 if (!is_eq_stride(aff, i))
2042                         continue;
2043                 k = isl_basic_map_alloc_equality(bmap);
2044                 if (k < 0)
2045                         goto error;
2046                 isl_seq_clr(bmap->eq[k], 1 + nparam);
2047                 isl_seq_cpy(bmap->eq[k] + 1 + nparam + d,
2048                                 aff->eq[i] + 1 + nparam, d);
2049                 isl_seq_neg(bmap->eq[k] + 1 + nparam,
2050                                 aff->eq[i] + 1 + nparam, d);
2051                 isl_seq_cpy(bmap->eq[k] + 1 + nparam + 2 * d,
2052                                 aff->eq[i] + 1 + nparam + d, aff->n_div);
2053                 isl_int_set_si(bmap->eq[k][1 + total + aff->n_div], 0);
2054         }
2055         obj = isl_vec_alloc(map->ctx, 1 + nparam + d);
2056         if (!obj)
2057                 goto error;
2058         isl_seq_clr(obj->el, 1 + nparam + d);
2059         for (i = 0; i < d; ++ i) {
2060                 enum isl_lp_result res;
2061
2062                 isl_int_set_si(obj->el[1 + nparam + i], 1);
2063
2064                 res = isl_set_solve_lp(delta, 0, obj->el, map->ctx->one, &opt,
2065                                         NULL, NULL);
2066                 if (res == isl_lp_error)
2067                         goto error;
2068                 if (res == isl_lp_ok) {
2069                         k = isl_basic_map_alloc_inequality(bmap);
2070                         if (k < 0)
2071                                 goto error;
2072                         isl_seq_clr(bmap->ineq[k],
2073                                         1 + nparam + 2 * d + bmap->n_div);
2074                         isl_int_set_si(bmap->ineq[k][1 + nparam + i], -1);
2075                         isl_int_set_si(bmap->ineq[k][1 + nparam + d + i], 1);
2076                         isl_int_neg(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], opt);
2077                 }
2078
2079                 res = isl_set_solve_lp(delta, 1, obj->el, map->ctx->one, &opt,
2080                                         NULL, NULL);
2081                 if (res == isl_lp_error)
2082                         goto error;
2083                 if (res == isl_lp_ok) {
2084                         k = isl_basic_map_alloc_inequality(bmap);
2085                         if (k < 0)
2086                                 goto error;
2087                         isl_seq_clr(bmap->ineq[k],
2088                                         1 + nparam + 2 * d + bmap->n_div);
2089                         isl_int_set_si(bmap->ineq[k][1 + nparam + i], 1);
2090                         isl_int_set_si(bmap->ineq[k][1 + nparam + d + i], -1);
2091                         isl_int_set(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], opt);
2092                 }
2093
2094                 isl_int_set_si(obj->el[1 + nparam + i], 0);
2095         }
2096         k = isl_basic_map_alloc_inequality(bmap);
2097         if (k < 0)
2098                 goto error;
2099         isl_seq_clr(bmap->ineq[k],
2100                         1 + nparam + 2 * d + bmap->n_div);
2101         if (!with_id)
2102                 isl_int_set_si(bmap->ineq[k][0], -1);
2103         isl_int_set_si(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], 1);
2104
2105         app = isl_map_from_domain_and_range(dom, ran);
2106
2107         isl_vec_free(obj);
2108         isl_basic_set_free(aff);
2109         isl_map_free(map);
2110         bmap = isl_basic_map_finalize(bmap);
2111         isl_set_free(delta);
2112         isl_int_clear(opt);
2113
2114         map = isl_map_from_basic_map(bmap);
2115         map = isl_map_intersect(map, app);
2116
2117         return map;
2118 error:
2119         isl_vec_free(obj);
2120         isl_basic_map_free(bmap);
2121         isl_basic_set_free(aff);
2122         isl_set_free(dom);
2123         isl_set_free(ran);
2124         isl_map_free(map);
2125         isl_set_free(delta);
2126         isl_int_clear(opt);
2127         return NULL;
2128 }
2129
2130 /* Given a map, compute the smallest superset of this map that is of the form
2131  *
2132  *      { i -> j : L <= j - i <= U and exists a_p: j_p - i_p = M_p a_p }
2133  *
2134  * (where p ranges over the (non-parametric) dimensions),
2135  * compute the transitive closure of this map, i.e.,
2136  *
2137  *      { i -> j : exists k > 0:
2138  *              k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p }
2139  *
2140  * and intersect domain and range of this transitive closure with
2141  * domain and range of the original map.
2142  */
2143 static __isl_give isl_map *box_closure(__isl_take isl_map *map)
2144 {
2145         isl_set *domain;
2146         isl_set *range;
2147
2148         domain = isl_map_domain(isl_map_copy(map));
2149         domain = isl_set_coalesce(domain);
2150         range = isl_map_range(isl_map_copy(map));
2151         range = isl_set_coalesce(range);
2152
2153         return box_closure_on_domain(map, domain, range, 0);
2154 }
2155
2156 /* Given a map, compute the smallest superset of this map that is of the form
2157  *
2158  *      { i -> j : L <= j - i <= U and exists a_p: j_p - i_p = M_p a_p }
2159  *
2160  * (where p ranges over the (non-parametric) dimensions),
2161  * compute the transitive and partially reflexive closure of this map, i.e.,
2162  *
2163  *      { i -> j : exists k >= 0:
2164  *              k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p }
2165  *
2166  * and intersect domain and range of this transitive closure with
2167  * the given domain.
2168  */
2169 static __isl_give isl_map *box_closure_with_identity(__isl_take isl_map *map,
2170         __isl_take isl_set *dom)
2171 {
2172         return box_closure_on_domain(map, dom, isl_set_copy(dom), 1);
2173 }
2174
2175 /* Check whether app is the transitive closure of map.
2176  * In particular, check that app is acyclic and, if so,
2177  * check that
2178  *
2179  *      app \subset (map \cup (map \circ app))
2180  */
2181 static int check_exactness_omega(__isl_keep isl_map *map,
2182         __isl_keep isl_map *app)
2183 {
2184         isl_set *delta;
2185         int i;
2186         int is_empty, is_exact;
2187         unsigned d;
2188         isl_map *test;
2189
2190         delta = isl_map_deltas(isl_map_copy(app));
2191         d = isl_set_dim(delta, isl_dim_set);
2192         for (i = 0; i < d; ++i)
2193                 delta = isl_set_fix_si(delta, isl_dim_set, i, 0);
2194         is_empty = isl_set_is_empty(delta);
2195         isl_set_free(delta);
2196         if (is_empty < 0)
2197                 return -1;
2198         if (!is_empty)
2199                 return 0;
2200
2201         test = isl_map_apply_range(isl_map_copy(app), isl_map_copy(map));
2202         test = isl_map_union(test, isl_map_copy(map));
2203         is_exact = isl_map_is_subset(app, test);
2204         isl_map_free(test);
2205
2206         return is_exact;
2207 }
2208
2209 /* Check if basic map M_i can be combined with all the other
2210  * basic maps such that
2211  *
2212  *      (\cup_j M_j)^+
2213  *
2214  * can be computed as
2215  *
2216  *      M_i \cup (\cup_{j \ne i} M_i^* \circ M_j \circ M_i^*)^+
2217  *
2218  * In particular, check if we can compute a compact representation
2219  * of
2220  *
2221  *              M_i^* \circ M_j \circ M_i^*
2222  *
2223  * for each j != i.
2224  * Let M_i^? be an extension of M_i^+ that allows paths
2225  * of length zero, i.e., the result of box_closure(., 1).
2226  * The criterion, as proposed by Kelly et al., is that
2227  * id = M_i^? - M_i^+ can be represented as a basic map
2228  * and that
2229  *
2230  *      id \circ M_j \circ id = M_j
2231  *
2232  * for each j != i.
2233  *
2234  * If this function returns 1, then tc and qc are set to
2235  * M_i^+ and M_i^?, respectively.
2236  */
2237 static int can_be_split_off(__isl_keep isl_map *map, int i,
2238         __isl_give isl_map **tc, __isl_give isl_map **qc)
2239 {
2240         isl_map *map_i, *id = NULL;
2241         int j = -1;
2242         isl_set *C;
2243
2244         *tc = NULL;
2245         *qc = NULL;
2246
2247         C = isl_set_union(isl_map_domain(isl_map_copy(map)),
2248                           isl_map_range(isl_map_copy(map)));
2249         C = isl_set_from_basic_set(isl_set_simple_hull(C));
2250         if (!C)
2251                 goto error;
2252
2253         map_i = isl_map_from_basic_map(isl_basic_map_copy(map->p[i]));
2254         *tc = box_closure(isl_map_copy(map_i));
2255         *qc = box_closure_with_identity(map_i, C);
2256         id = isl_map_subtract(isl_map_copy(*qc), isl_map_copy(*tc));
2257
2258         if (!id || !*qc)
2259                 goto error;
2260         if (id->n != 1 || (*qc)->n != 1)
2261                 goto done;
2262
2263         for (j = 0; j < map->n; ++j) {
2264                 isl_map *map_j, *test;
2265                 int is_ok;
2266
2267                 if (i == j)
2268                         continue;
2269                 map_j = isl_map_from_basic_map(
2270                                         isl_basic_map_copy(map->p[j]));
2271                 test = isl_map_apply_range(isl_map_copy(id),
2272                                                 isl_map_copy(map_j));
2273                 test = isl_map_apply_range(test, isl_map_copy(id));
2274                 is_ok = isl_map_is_equal(test, map_j);
2275                 isl_map_free(map_j);
2276                 isl_map_free(test);
2277                 if (is_ok < 0)
2278                         goto error;
2279                 if (!is_ok)
2280                         break;
2281         }
2282
2283 done:
2284         isl_map_free(id);
2285         if (j == map->n)
2286                 return 1;
2287
2288         isl_map_free(*qc);
2289         isl_map_free(*tc);
2290         *qc = NULL;
2291         *tc = NULL;
2292
2293         return 0;
2294 error:
2295         isl_map_free(id);
2296         isl_map_free(*qc);
2297         isl_map_free(*tc);
2298         *qc = NULL;
2299         *tc = NULL;
2300         return -1;
2301 }
2302
2303 static __isl_give isl_map *box_closure_with_check(__isl_take isl_map *map,
2304         int *exact)
2305 {
2306         isl_map *app;
2307
2308         app = box_closure(isl_map_copy(map));
2309         if (exact)
2310                 *exact = check_exactness_omega(map, app);
2311
2312         isl_map_free(map);
2313         return app;
2314 }
2315
2316 /* Compute an overapproximation of the transitive closure of "map"
2317  * using a variation of the algorithm from
2318  * "Transitive Closure of Infinite Graphs and its Applications"
2319  * by Kelly et al.
2320  *
2321  * We first check whether we can can split of any basic map M_i and
2322  * compute
2323  *
2324  *      (\cup_j M_j)^+
2325  *
2326  * as
2327  *
2328  *      M_i \cup (\cup_{j \ne i} M_i^* \circ M_j \circ M_i^*)^+
2329  *
2330  * using a recursive call on the remaining map.
2331  *
2332  * If not, we simply call box_closure on the whole map.
2333  */
2334 static __isl_give isl_map *transitive_closure_omega(__isl_take isl_map *map,
2335         int *exact)
2336 {
2337         int i, j;
2338         int exact_i;
2339         isl_map *app;
2340
2341         if (!map)
2342                 return NULL;
2343         if (map->n == 1)
2344                 return box_closure_with_check(map, exact);
2345
2346         for (i = 0; i < map->n; ++i) {
2347                 int ok;
2348                 isl_map *qc, *tc;
2349                 ok = can_be_split_off(map, i, &tc, &qc);
2350                 if (ok < 0)
2351                         goto error;
2352                 if (!ok)
2353                         continue;
2354
2355                 app = isl_map_alloc_dim(isl_map_get_dim(map), map->n - 1, 0);
2356
2357                 for (j = 0; j < map->n; ++j) {
2358                         if (j == i)
2359                                 continue;
2360                         app = isl_map_add_basic_map(app,
2361                                                 isl_basic_map_copy(map->p[j]));
2362                 }
2363
2364                 app = isl_map_apply_range(isl_map_copy(qc), app);
2365                 app = isl_map_apply_range(app, qc);
2366
2367                 app = isl_map_union(tc, transitive_closure_omega(app, NULL));
2368                 exact_i = check_exactness_omega(map, app);
2369                 if (exact_i == 1) {
2370                         if (exact)
2371                                 *exact = exact_i;
2372                         isl_map_free(map);
2373                         return app;
2374                 }
2375                 isl_map_free(app);
2376                 if (exact_i < 0)
2377                         goto error;
2378         }
2379
2380         return box_closure_with_check(map, exact);
2381 error:
2382         isl_map_free(map);
2383         return NULL;
2384 }
2385
2386 int isl_map_is_transitively_closed(__isl_keep isl_map *map)
2387 {
2388         isl_map *map2;
2389         int closed;
2390
2391         map2 = isl_map_apply_range(isl_map_copy(map), isl_map_copy(map));
2392         closed = isl_map_is_subset(map2, map);
2393         isl_map_free(map2);
2394
2395         return closed;
2396 }
2397
2398 /* Compute the transitive closure  of "map", or an overapproximation.
2399  * If the result is exact, then *exact is set to 1.
2400  * Simply use map_power to compute the powers of map, but tell
2401  * it to project out the lengths of the paths instead of equating
2402  * the length to a parameter.
2403  */
2404 __isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
2405         int *exact)
2406 {
2407         unsigned param;
2408         int closed;
2409
2410         if (!map)
2411                 goto error;
2412
2413         if (map->ctx->opt->closure == ISL_CLOSURE_OMEGA)
2414                 return transitive_closure_omega(map, exact);
2415
2416         map = isl_map_compute_divs(map);
2417         map = isl_map_coalesce(map);
2418         closed = isl_map_is_transitively_closed(map);
2419         if (closed < 0)
2420                 goto error;
2421         if (closed) {
2422                 if (exact)
2423                         *exact = 1;
2424                 return map;
2425         }
2426
2427         param = isl_map_dim(map, isl_dim_param);
2428         map = map_power(map, param, exact, 1);
2429
2430         return map;
2431 error:
2432         isl_map_free(map);
2433         return NULL;
2434 }