isl_transitive_closure.c

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