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