isl_transitive_closure.c

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