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