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
  14 /*
  15  * The transitive closure implementation is based on the paper
  16  * "Computing the Transitive Closure of a Union of Affine Integer
  17  * Tuple Relations" by Anna Beletska, Denis Barthou, Wlodzimierz Bielecki and
  18  * Albert Cohen.
  19  */
  20
  21 /* Given a set of n offsets v_i (the rows of "steps"), construct a relation
  22  * of the given dimension specification (Z^{n+1} -> Z^{n+1})
  23  * that maps an element x to any element that can be reached
  24  * by taking a non-negative number of steps along any of
  25  * the extended offsets v'_i = [v_i 1].
  26  * That is, construct
  27  *
  28  * { [x] -> [y] : exists k_i >= 0, y = x + \sum_i k_i v'_i }
  29  *
  30  * For any element in this relation, the number of steps taken
  31  * is equal to the difference in the final coordinates.
  32  */
  33 static __isl_give isl_map *path_along_steps(__isl_take isl_dim *dim,
  34         __isl_keep isl_mat *steps)
  35 {
  36         int i, j, k;
  37         struct isl_basic_map *path = NULL;
  38         unsigned d;
  39         unsigned n;
  40         unsigned nparam;
  41
  42         if (!dim || !steps)
  43                 goto error;
  44
  45         d = isl_dim_size(dim, isl_dim_in);
  46         n = steps->n_row;
  47         nparam = isl_dim_size(dim, isl_dim_param);
  48
  49         path = isl_basic_map_alloc_dim(isl_dim_copy(dim), n, d, n);
  50
  51         for (i = 0; i < n; ++i) {
  52                 k = isl_basic_map_alloc_div(path);
  53                 if (k < 0)
  54                         goto error;
  55                 isl_assert(steps->ctx, i == k, goto error);
  56                 isl_int_set_si(path->div[k][0], 0);
  57         }
  58
  59         for (i = 0; i < d; ++i) {
  60                 k = isl_basic_map_alloc_equality(path);
  61                 if (k < 0)
  62                         goto error;
  63                 isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
  64                 isl_int_set_si(path->eq[k][1 + nparam + i], 1);
  65                 isl_int_set_si(path->eq[k][1 + nparam + d + i], -1);
  66                 if (i == d - 1)
  67                         for (j = 0; j < n; ++j)
  68                                 isl_int_set_si(path->eq[k][1 + nparam + 2 * d + j], 1);
  69                 else
  70                         for (j = 0; j < n; ++j)
  71                                 isl_int_set(path->eq[k][1 + nparam + 2 * d + j],
  72                                             steps->row[j][i]);
  73         }
  74
  75         for (i = 0; i < n; ++i) {
  76                 k = isl_basic_map_alloc_inequality(path);
  77                 if (k < 0)
  78                         goto error;
  79                 isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
  80                 isl_int_set_si(path->ineq[k][1 + nparam + 2 * d + i], 1);
  81         }
  82
  83         isl_dim_free(dim);
  84
  85         path = isl_basic_map_simplify(path);
  86         path = isl_basic_map_finalize(path);
  87         return isl_map_from_basic_map(path);
  88 error:
  89         isl_dim_free(dim);
  90         isl_basic_map_free(path);
  91         return NULL;
  92 }
  93
  94 /* Given a dimenion specification Z^{n+1} -> Z^{n+1} and a parameter "param",
  95  * construct a map that equates the parameter to the difference
  96  * in the final coordinates and imposes that this difference is positive.
  97  * That is, construct
  98  *
  99  *      { [x,x_s] -> [y,y_s] : k = y_s - x_s > 0 }
 100  */
 101 static __isl_give isl_map *equate_parameter_to_length(__isl_take isl_dim *dim,
 102         unsigned param)
 103 {
 104         struct isl_basic_map *bmap;
 105         unsigned d;
 106         unsigned nparam;
 107         int k;
 108
 109         d = isl_dim_size(dim, isl_dim_in);
 110         nparam = isl_dim_size(dim, isl_dim_param);
 111         bmap = isl_basic_map_alloc_dim(dim, 0, 1, 1);
 112         k = isl_basic_map_alloc_equality(bmap);
 113         if (k < 0)
 114                 goto error;
 115         isl_seq_clr(bmap->eq[k], 1 + isl_basic_map_total_dim(bmap));
 116         isl_int_set_si(bmap->eq[k][1 + param], -1);
 117         isl_int_set_si(bmap->eq[k][1 + nparam + d - 1], -1);
 118         isl_int_set_si(bmap->eq[k][1 + nparam + d + d - 1], 1);
 119
 120         k = isl_basic_map_alloc_inequality(bmap);
 121         if (k < 0)
 122                 goto error;
 123         isl_seq_clr(bmap->ineq[k], 1 + isl_basic_map_total_dim(bmap));
 124         isl_int_set_si(bmap->ineq[k][1 + param], 1);
 125         isl_int_set_si(bmap->ineq[k][0], -1);
 126
 127         bmap = isl_basic_map_finalize(bmap);
 128         return isl_map_from_basic_map(bmap);
 129 error:
 130         isl_basic_map_free(bmap);
 131         return NULL;
 132 }
 133
 134 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D,
 135  * construct a map that is an overapproximation of the map
 136  * that takes an element from the space D to another
 137  * element from the same space, such that the difference between
 138  * them is a strictly positive sum of differences between images
 139  * and pre-images in one of the R_i.
 140  * The number of differences in the sum is equated to parameter "param".
 141  * That is, let
 142  *
 143  *      \Delta_i = { y - x | (x, y) in R_i }
 144  *
 145  * then the constructed map is an overapproximation of
 146  *
 147  *      { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
 148  *                              d = \sum_i k_i and k = \sum_i k_i > 0 }
 149  *
 150  * We first construct an extended mapping with an extra coordinate
 151  * that indicates the number of steps taken.  In particular,
 152  * the difference in the last coordinate is equal to the number
 153  * of steps taken to move from a domain element to the corresponding
 154  * image element(s).
 155  * In the final step, this difference is equated to the parameter "param"
 156  * and made positive.  The extra coordinates are subsequently projected out.
 157  *
 158  * The elements of the singleton \Delta_i's are collected as the
 159  * rows of the steps matrix.  For all these \Delta_i's together,
 160  * a single path is constructed.
 161  * For each of the other \Delta_i's
 162  * we currently simply construct a universal map { (x) -> (y) }.
 163  */
 164 static __isl_give isl_map *construct_path(__isl_keep isl_map *map,
 165         unsigned param)
 166 {
 167         struct isl_mat *steps = NULL;
 168         struct isl_map *path = NULL;
 169         struct isl_map *diff;
 170         struct isl_dim *dim = NULL;
 171         unsigned d;
 172         int i, j, n;
 173
 174         if (!map)
 175                 return NULL;
 176
 177         dim = isl_map_get_dim(map);
 178
 179         d = isl_dim_size(dim, isl_dim_in);
 180         dim = isl_dim_add(dim, isl_dim_in, 1);
 181         dim = isl_dim_add(dim, isl_dim_out, 1);
 182
 183         path = isl_map_identity(isl_dim_domain(isl_dim_copy(dim)));
 184
 185         steps = isl_mat_alloc(map->ctx, map->n, d);
 186         if (!steps)
 187                 goto error;
 188
 189         n = 0;
 190         for (i = 0; i < map->n; ++i) {
 191                 struct isl_basic_set *delta;
 192
 193                 delta = isl_basic_map_deltas(isl_basic_map_copy(map->p[i]));
 194
 195                 for (j = 0; j < d; ++j) {
 196                         int fixed;
 197
 198                         fixed = isl_basic_set_fast_dim_is_fixed(delta, j,
 199                                                             &steps->row[n][j]);
 200                         if (fixed < 0) {
 201                                 isl_basic_set_free(delta);
 202                                 goto error;
 203                         }
 204                         if (!fixed)
 205                                 break;
 206                 }
 207
 208                 isl_basic_set_free(delta);
 209
 210                 if (j < d) {
 211                         isl_map_free(path);
 212                         path = isl_map_universe(isl_dim_copy(dim));
 213                 } else
 214                         ++n;
 215         }
 216
 217         if (n > 0) {
 218                 steps->n_row = n;
 219                 path = isl_map_apply_range(path,
 220                                 path_along_steps(isl_dim_copy(dim), steps));
 221         }
 222
 223         diff = equate_parameter_to_length(dim, param);
 224         path = isl_map_intersect(path, diff);
 225         path = isl_map_project_out(path, isl_dim_in, d, 1);
 226         path = isl_map_project_out(path, isl_dim_out, d, 1);
 227
 228         isl_mat_free(steps);
 229         return path;
 230 error:
 231         isl_dim_free(dim);
 232         isl_map_free(path);
 233         return NULL;
 234 }
 235
 236 /* Check whether "path" is acyclic.
 237  * That is, check whether
 238  *
 239  *      { d | d = y - x and (x,y) in path }
 240  *
 241  * does not contain the origin.
 242  */
 243 static int is_acyclic(__isl_take isl_map *path)
 244 {
 245         int i;
 246         int acyclic;
 247         unsigned dim;
 248         struct isl_set *delta;
 249
 250         delta = isl_map_deltas(path);
 251         dim = isl_set_dim(delta, isl_dim_set);
 252         for (i = 0; i < dim; ++i)
 253                 delta = isl_set_fix_si(delta, isl_dim_set, i, 0);
 254
 255         acyclic = isl_set_is_empty(delta);
 256         isl_set_free(delta);
 257
 258         return acyclic;
 259 }
 260
 261 /* Shift variable at position "pos" up by one.
 262  * That is, replace the corresponding variable v by v - 1.
 263  */
 264 static __isl_give isl_basic_map *basic_map_shift_pos(
 265         __isl_take isl_basic_map *bmap, unsigned pos)
 266 {
 267         int i;
 268
 269         bmap = isl_basic_map_cow(bmap);
 270         if (!bmap)
 271                 return NULL;
 272
 273         for (i = 0; i < bmap->n_eq; ++i)
 274                 isl_int_sub(bmap->eq[i][0], bmap->eq[i][0], bmap->eq[i][pos]);
 275
 276         for (i = 0; i < bmap->n_ineq; ++i)
 277                 isl_int_sub(bmap->ineq[i][0],
 278                                 bmap->ineq[i][0], bmap->ineq[i][pos]);
 279
 280         for (i = 0; i < bmap->n_div; ++i) {
 281                 if (isl_int_is_zero(bmap->div[i][0]))
 282                         continue;
 283                 isl_int_sub(bmap->div[i][1],
 284                                 bmap->div[i][1], bmap->div[i][1 + pos]);
 285         }
 286
 287         return bmap;
 288 }
 289
 290 /* Shift variable at position "pos" up by one.
 291  * That is, replace the corresponding variable v by v - 1.
 292  */
 293 static __isl_give isl_map *map_shift_pos(__isl_take isl_map *map, unsigned pos)
 294 {
 295         int i;
 296
 297         map = isl_map_cow(map);
 298         if (!map)
 299                 return NULL;
 300
 301         for (i = 0; i < map->n; ++i) {
 302                 map->p[i] = basic_map_shift_pos(map->p[i], pos);
 303                 if (!map->p[i])
 304                         goto error;
 305         }
 306         ISL_F_CLR(map, ISL_MAP_NORMALIZED);
 307         return map;
 308 error:
 309         isl_map_free(map);
 310         return NULL;
 311 }
 312
 313 /* Check whether the overapproximation of the power of "map" is exactly
 314  * the power of "map".  Let R be "map" and A_k the overapproximation.
 315  * The approximation is exact if
 316  *
 317  *      A_1 = R
 318  *      A_k = A_{k-1} \circ R                   k >= 2
 319  *
 320  * Since A_k is known to be an overapproximation, we only need to check
 321  *
 322  *      A_1 \subset R
 323  *      A_k \subset A_{k-1} \circ R             k >= 2
 324  *
 325  */
 326 static int check_power_exactness(__isl_take isl_map *map,
 327         __isl_take isl_map *app, unsigned param)
 328 {
 329         int exact;
 330         isl_map *app_1;
 331         isl_map *app_2;
 332
 333         app_1 = isl_map_fix_si(isl_map_copy(app), isl_dim_param, param, 1);
 334
 335         exact = isl_map_is_subset(app_1, map);
 336         isl_map_free(app_1);
 337
 338         if (!exact || exact < 0) {
 339                 isl_map_free(app);
 340                 isl_map_free(map);
 341                 return exact;
 342         }
 343
 344         app_2 = isl_map_lower_bound_si(isl_map_copy(app),
 345                                         isl_dim_param, param, 2);
 346         app_1 = map_shift_pos(app, 1 + param);
 347         app_1 = isl_map_apply_range(map, app_1);
 348
 349         exact = isl_map_is_subset(app_2, app_1);
 350
 351         isl_map_free(app_1);
 352         isl_map_free(app_2);
 353
 354         return exact;
 355 }
 356
 357 /* Check whether the overapproximation of the power of "map" is exactly
 358  * the power of "map", possibly after projecting out the power (if "project"
 359  * is set).
 360  *
 361  * If "project" is set and if "steps" can only result in acyclic paths,
 362  * then we check
 363  *
 364  *      A = R \cup (A \circ R)
 365  *
 366  * where A is the overapproximation with the power projected out, i.e.,
 367  * an overapproximation of the transitive closure.
 368  * More specifically, since A is known to be an overapproximation, we check
 369  *
 370  *      A \subset R \cup (A \circ R)
 371  *
 372  * Otherwise, we check if the power is exact.
 373  */
 374 static int check_exactness(__isl_take isl_map *map, __isl_take isl_map *app,
 375         __isl_take isl_map *path, unsigned param, int project)
 376 {
 377         isl_map *test;
 378         int exact;
 379
 380         if (project) {
 381                 project = is_acyclic(path);
 382                 if (project < 0)
 383                         goto error;
 384         } else
 385                 isl_map_free(path);
 386
 387         if (!project)
 388                 return check_power_exactness(map, app, param);
 389
 390         map = isl_map_project_out(map, isl_dim_param, param, 1);
 391         app = isl_map_project_out(app, isl_dim_param, param, 1);
 392
 393         test = isl_map_apply_range(isl_map_copy(map), isl_map_copy(app));
 394         test = isl_map_union(test, isl_map_copy(map));
 395
 396         exact = isl_map_is_subset(app, test);
 397
 398         isl_map_free(app);
 399         isl_map_free(test);
 400
 401         isl_map_free(map);
 402
 403         return exact;
 404 error:
 405         isl_map_free(app);
 406         isl_map_free(map);
 407         return -1;
 408 }
 409
 410 /* Compute the positive powers of "map", or an overapproximation.
 411  * The power is given by parameter "param".  If the result is exact,
 412  * then *exact is set to 1.
 413  * If project is set, then we are actually interested in the transitive
 414  * closure, so we can use a more relaxed exactness check.
 415  */
 416 static __isl_give isl_map *map_power(__isl_take isl_map *map, unsigned param,
 417         int *exact, int project)
 418 {
 419         struct isl_set *domain = NULL;
 420         struct isl_set *range = NULL;
 421         struct isl_map *app = NULL;
 422         struct isl_map *path = NULL;
 423
 424         if (exact)
 425                 *exact = 1;
 426
 427         map = isl_map_remove_empty_parts(map);
 428         if (!map)
 429                 return NULL;
 430
 431         if (isl_map_fast_is_empty(map))
 432                 return map;
 433
 434         isl_assert(map->ctx, param < isl_map_dim(map, isl_dim_param), goto error);
 435         isl_assert(map->ctx,
 436                 isl_map_dim(map, isl_dim_in) == isl_map_dim(map, isl_dim_out),
 437                 goto error);
 438
 439         domain = isl_map_domain(isl_map_copy(map));
 440         domain = isl_set_coalesce(domain);
 441         range = isl_map_range(isl_map_copy(map));
 442         range = isl_set_coalesce(range);
 443         app = isl_map_from_domain_and_range(isl_set_copy(domain),
 444                                             isl_set_copy(range));
 445
 446         path = construct_path(map, param);
 447         app = isl_map_intersect(app, isl_map_copy(path));
 448
 449         if (exact &&
 450             (*exact = check_exactness(isl_map_copy(map), isl_map_copy(app),
 451                                   isl_map_copy(path), param, project)) < 0)
 452                 goto error;
 453
 454         isl_set_free(domain);
 455         isl_set_free(range);
 456         isl_map_free(path);
 457         isl_map_free(map);
 458         return app;
 459 error:
 460         isl_set_free(domain);
 461         isl_set_free(range);
 462         isl_map_free(path);
 463         isl_map_free(map);
 464         isl_map_free(app);
 465         return NULL;
 466 }
 467
 468 /* Compute the positive powers of "map", or an overapproximation.
 469  * The power is given by parameter "param".  If the result is exact,
 470  * then *exact is set to 1.
 471  */
 472 __isl_give isl_map *isl_map_power(__isl_take isl_map *map, unsigned param,
 473         int *exact)
 474 {
 475         return map_power(map, param, exact, 0);
 476 }
 477
 478 /* Compute the transitive closure  of "map", or an overapproximation.
 479  * If the result is exact, then *exact is set to 1.
 480  * Simply compute the powers of map and then project out the parameter
 481  * describing the power.
 482  */
 483 __isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
 484         int *exact)
 485 {
 486         unsigned param;
 487
 488         if (!map)
 489                 goto error;
 490
 491         param = isl_map_dim(map, isl_dim_param);
 492         map = isl_map_add(map, isl_dim_param, 1);
 493         map = map_power(map, param, exact, 1);
 494         map = isl_map_project_out(map, isl_dim_param, param, 1);
 495
 496         return map;
 497 error:
 498         isl_map_free(map);
 499         return NULL;
 500 }