isl_map_subtract.c

   1 /*
   2  * Copyright 2008-2009 Katholieke Universiteit Leuven
   3  *
   4  * Use of this software is governed by the GNU LGPLv2.1 license
   5  *
   6  * Written by Sven Verdoolaege, K.U.Leuven, Departement
   7  * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
   8  */
   9
  10 #include "isl_seq.h"
  11 #include "isl_set.h"
  12 #include "isl_map.h"
  13 #include "isl_map_private.h"
  14 #include "isl_tab.h"
  15
  16 static void expand_constraint(isl_vec *v, unsigned dim,
  17         isl_int *c, int *div_map, unsigned n_div)
  18 {
  19         int i;
  20
  21         isl_seq_cpy(v->el, c, 1 + dim);
  22         isl_seq_clr(v->el + 1 + dim, v->size - (1 + dim));
  23
  24         for (i = 0; i < n_div; ++i)
  25                 isl_int_set(v->el[1 + dim + div_map[i]], c[1 + dim + i]);
  26 }
  27
  28 /* Add all constraints of bmap to tab.  The equalities of bmap
  29  * are added as a pair of inequalities.
  30  */
  31 static int tab_add_constraints(struct isl_tab *tab,
  32         __isl_keep isl_basic_map *bmap, int *div_map)
  33 {
  34         int i;
  35         unsigned dim;
  36         unsigned tab_total;
  37         unsigned bmap_total;
  38         isl_vec *v;
  39
  40         if (!tab || !bmap)
  41                 return -1;
  42
  43         tab_total = isl_basic_map_total_dim(tab->bmap);
  44         bmap_total = isl_basic_map_total_dim(bmap);
  45         dim = isl_dim_total(tab->bmap->dim);
  46
  47         if (isl_tab_extend_cons(tab, 2 * bmap->n_eq + bmap->n_ineq) < 0)
  48                 return -1;
  49
  50         v = isl_vec_alloc(bmap->ctx, 1 + tab_total);
  51
  52         for (i = 0; i < bmap->n_eq; ++i) {
  53                 expand_constraint(v, dim, bmap->eq[i], div_map, bmap->n_div);
  54                 if (isl_tab_add_ineq(tab, v->el) < 0)
  55                         goto error;
  56                 isl_seq_neg(bmap->eq[i], bmap->eq[i], 1 + bmap_total);
  57                 expand_constraint(v, dim, bmap->eq[i], div_map, bmap->n_div);
  58                 if (isl_tab_add_ineq(tab, v->el) < 0)
  59                         goto error;
  60                 isl_seq_neg(bmap->eq[i], bmap->eq[i], 1 + bmap_total);
  61                 if (tab->empty)
  62                         break;
  63         }
  64
  65         for (i = 0; i < bmap->n_ineq; ++i) {
  66                 expand_constraint(v, dim, bmap->ineq[i], div_map, bmap->n_div);
  67                 if (isl_tab_add_ineq(tab, v->el) < 0)
  68                         goto error;
  69                 if (tab->empty)
  70                         break;
  71         }
  72
  73         isl_vec_free(v);
  74         return 0;
  75 error:
  76         isl_vec_free(v);
  77         return -1;
  78 }
  79
  80 /* Add a specific constraint of bmap (or its opposite) to tab.
  81  * The position of the constraint is specified by "c", where
  82  * the equalities of bmap are counted twice, once for the inequality
  83  * that is equal to the equality, and once for its negation.
  84  */
  85 static int tab_add_constraint(struct isl_tab *tab,
  86         __isl_keep isl_basic_map *bmap, int *div_map, int c, int oppose)
  87 {
  88         unsigned dim;
  89         unsigned tab_total;
  90         unsigned bmap_total;
  91         isl_vec *v;
  92         int r;
  93
  94         if (!tab || !bmap)
  95                 return -1;
  96
  97         tab_total = isl_basic_map_total_dim(tab->bmap);
  98         bmap_total = isl_basic_map_total_dim(bmap);
  99         dim = isl_dim_total(tab->bmap->dim);
 100
 101         v = isl_vec_alloc(bmap->ctx, 1 + tab_total);
 102         if (!v)
 103                 return -1;
 104
 105         if (c < 2 * bmap->n_eq) {
 106                 if ((c % 2) != oppose)
 107                         isl_seq_neg(bmap->eq[c/2], bmap->eq[c/2],
 108                                         1 + bmap_total);
 109                 if (oppose)
 110                         isl_int_sub_ui(bmap->eq[c/2][0], bmap->eq[c/2][0], 1);
 111                 expand_constraint(v, dim, bmap->eq[c/2], div_map, bmap->n_div);
 112                 r = isl_tab_add_ineq(tab, v->el);
 113                 if (oppose)
 114                         isl_int_add_ui(bmap->eq[c/2][0], bmap->eq[c/2][0], 1);
 115                 if ((c % 2) != oppose)
 116                         isl_seq_neg(bmap->eq[c/2], bmap->eq[c/2],
 117                                         1 + bmap_total);
 118         } else {
 119                 c -= 2 * bmap->n_eq;
 120                 if (oppose) {
 121                         isl_seq_neg(bmap->ineq[c], bmap->ineq[c],
 122                                         1 + bmap_total);
 123                         isl_int_sub_ui(bmap->ineq[c][0], bmap->ineq[c][0], 1);
 124                 }
 125                 expand_constraint(v, dim, bmap->ineq[c], div_map, bmap->n_div);
 126                 r = isl_tab_add_ineq(tab, v->el);
 127                 if (oppose) {
 128                         isl_int_add_ui(bmap->ineq[c][0], bmap->ineq[c][0], 1);
 129                         isl_seq_neg(bmap->ineq[c], bmap->ineq[c],
 130                                         1 + bmap_total);
 131                 }
 132         }
 133
 134         isl_vec_free(v);
 135         return r;
 136 }
 137
 138 static int tab_add_divs(struct isl_tab *tab, __isl_keep isl_basic_map *bmap,
 139         int **div_map)
 140 {
 141         int i, j;
 142         struct isl_vec *vec;
 143         unsigned total;
 144         unsigned dim;
 145
 146         if (!bmap)
 147                 return -1;
 148         if (!bmap->n_div)
 149                 return 0;
 150
 151         if (!*div_map)
 152                 *div_map = isl_alloc_array(bmap->ctx, int, bmap->n_div);
 153         if (!*div_map)
 154                 return -1;
 155
 156         total = isl_basic_map_total_dim(tab->bmap);
 157         dim = total - tab->bmap->n_div;
 158         vec = isl_vec_alloc(bmap->ctx, 2 + total + bmap->n_div);
 159         if (!vec)
 160                 return -1;
 161
 162         for (i = 0; i < bmap->n_div; ++i) {
 163                 isl_seq_cpy(vec->el, bmap->div[i], 2 + dim);
 164                 isl_seq_clr(vec->el + 2 + dim, tab->bmap->n_div);
 165                 for (j = 0; j < i; ++j)
 166                         isl_int_set(vec->el[2 + dim + (*div_map)[j]],
 167                                         bmap->div[i][2 + dim + j]);
 168                 for (j = 0; j < tab->bmap->n_div; ++j)
 169                         if (isl_seq_eq(tab->bmap->div[j],
 170                                         vec->el, 2 + dim + tab->bmap->n_div))
 171                                 break;
 172                 (*div_map)[i] = j;
 173                 if (j == tab->bmap->n_div) {
 174                         vec->size = 2 + dim + tab->bmap->n_div;
 175                         if (isl_tab_add_div(tab, vec, NULL, NULL) < 0)
 176                                 goto error;
 177                 }
 178         }
 179
 180         isl_vec_free(vec);
 181
 182         return 0;
 183 error:
 184         isl_vec_free(vec);
 185
 186         return -1;
 187 }
 188
 189 /* Freeze all constraints of tableau tab.
 190  */
 191 static int tab_freeze_constraints(struct isl_tab *tab)
 192 {
 193         int i;
 194
 195         for (i = 0; i < tab->n_con; ++i)
 196                 if (isl_tab_freeze_constraint(tab, i) < 0)
 197                         return -1;
 198
 199         return 0;
 200 }
 201
 202 /* Check for redundant constraints starting at offset.
 203  * Put the indices of the redundant constraints in index
 204  * and return the number of redundant constraints.
 205  */
 206 static int n_non_redundant(struct isl_tab *tab, int offset, int **index)
 207 {
 208         int i, n;
 209         int n_test = tab->n_con - offset;
 210
 211         if (isl_tab_detect_redundant(tab) < 0)
 212                 return -1;
 213
 214         if (!*index)
 215                 *index = isl_alloc_array(tab->mat->ctx, int, n_test);
 216         if (!*index)
 217                 return -1;
 218
 219         for (n = 0, i = 0; i < n_test; ++i) {
 220                 int r;
 221                 r = isl_tab_is_redundant(tab, offset + i);
 222                 if (r < 0)
 223                         return -1;
 224                 if (r)
 225                         continue;
 226                 (*index)[n++] = i;
 227         }
 228
 229         return n;
 230 }
 231
 232 /* basic_map_collect_diff calls add on each of the pieces of
 233  * the set difference between bmap and map until the add method
 234  * return a negative value.
 235  */
 236 struct isl_diff_collector {
 237         int (*add)(struct isl_diff_collector *dc,
 238                     __isl_take isl_basic_map *bmap);
 239 };
 240
 241 /* Compute the set difference between bmap and map and call
 242  * dc->add on each of the piece until this function returns
 243  * a negative value.
 244  * Return 0 on success and -1 on error.  dc->add returning
 245  * a negative value is treated as an error, but the calling
 246  * function can interpret the results based on the state of dc.
 247  *
 248  * Assumes that both bmap and map have known divs.
 249  *
 250  * The difference is computed by a backtracking algorithm.
 251  * Each level corresponds to a basic map in "map".
 252  * When a node in entered for the first time, we check
 253  * if the corresonding basic map intersects the current piece
 254  * of "bmap".  If not, we move to the next level.
 255  * Otherwise, we split the current piece into as many
 256  * pieces as there are non-redundant constraints of the current
 257  * basic map in the intersection.  Each of these pieces is
 258  * handled by a child of the current node.
 259  * In particular, if there are n non-redundant constraints,
 260  * then for each 0 <= i < n, a piece is cut off by adding
 261  * constraints 0 <= j < i and adding the opposite of constraint i.
 262  * If there are no non-redundant constraints, meaning that the current
 263  * piece is a subset of the current basic map, then we simply backtrack.
 264  *
 265  * In the leaves, we check if the remaining piece has any integer points
 266  * and if so, pass it along to dc->add.  As a special case, if nothing
 267  * has been removed when we end up in a leaf, we simply pass along
 268  * the original basic map.
 269  */
 270 static int basic_map_collect_diff(__isl_take isl_basic_map *bmap,
 271         __isl_take isl_map *map, struct isl_diff_collector *dc)
 272 {
 273         int i;
 274         int modified;
 275         int level;
 276         int init;
 277         int empty;
 278         struct isl_tab *tab = NULL;
 279         struct isl_tab_undo **snap = NULL;
 280         int *k = NULL;
 281         int *n = NULL;
 282         int **index = NULL;
 283         int **div_map = NULL;
 284
 285         empty = isl_basic_map_is_empty(bmap);
 286         if (empty) {
 287                 isl_basic_map_free(bmap);
 288                 isl_map_free(map);
 289                 return empty < 0 ? -1 : 0;
 290         }
 291
 292         bmap = isl_basic_map_cow(bmap);
 293         map = isl_map_cow(map);
 294
 295         if (!bmap || !map)
 296                 goto error;
 297
 298         snap = isl_alloc_array(map->ctx, struct isl_tab_undo *, map->n);
 299         k = isl_alloc_array(map->ctx, int, map->n);
 300         n = isl_alloc_array(map->ctx, int, map->n);
 301         index = isl_calloc_array(map->ctx, int *, map->n);
 302         div_map = isl_calloc_array(map->ctx, int *, map->n);
 303         if (!snap || !k || !n || !index || !div_map)
 304                 goto error;
 305
 306         bmap = isl_basic_map_order_divs(bmap);
 307         map = isl_map_order_divs(map);
 308
 309         tab = isl_tab_from_basic_map(bmap);
 310         if (isl_tab_track_bmap(tab, isl_basic_map_copy(bmap)) < 0)
 311                 goto error;
 312
 313         modified = 0;
 314         level = 0;
 315         init = 1;
 316
 317         while (level >= 0) {
 318                 if (level >= map->n) {
 319                         int empty;
 320                         struct isl_basic_map *bm;
 321                         if (!modified) {
 322                                 if (dc->add(dc, isl_basic_map_copy(bmap)) < 0)
 323                                         goto error;
 324                                 break;
 325                         }
 326                         bm = isl_basic_map_copy(tab->bmap);
 327                         bm = isl_basic_map_cow(bm);
 328                         bm = isl_basic_map_update_from_tab(bm, tab);
 329                         bm = isl_basic_map_simplify(bm);
 330                         bm = isl_basic_map_finalize(bm);
 331                         empty = isl_basic_map_is_empty(bm);
 332                         if (empty)
 333                                 isl_basic_map_free(bm);
 334                         else if (dc->add(dc, bm) < 0)
 335                                 goto error;
 336                         if (empty < 0)
 337                                 goto error;
 338                         level--;
 339                         init = 0;
 340                         continue;
 341                 }
 342                 if (init) {
 343                         int offset;
 344                         struct isl_tab_undo *snap2;
 345                         snap2 = isl_tab_snap(tab);
 346                         if (tab_add_divs(tab, map->p[level],
 347                                          &div_map[level]) < 0)
 348                                 goto error;
 349                         offset = tab->n_con;
 350                         snap[level] = isl_tab_snap(tab);
 351                         if (tab_freeze_constraints(tab) < 0)
 352                                 goto error;
 353                         if (tab_add_constraints(tab, map->p[level],
 354                                                 div_map[level]) < 0)
 355                                 goto error;
 356                         k[level] = 0;
 357                         n[level] = 0;
 358                         if (tab->empty) {
 359                                 if (isl_tab_rollback(tab, snap2) < 0)
 360                                         goto error;
 361                                 level++;
 362                                 continue;
 363                         }
 364                         modified = 1;
 365                         n[level] = n_non_redundant(tab, offset, &index[level]);
 366                         if (n[level] < 0)
 367                                 goto error;
 368                         if (n[level] == 0) {
 369                                 level--;
 370                                 init = 0;
 371                                 continue;
 372                         }
 373                         if (isl_tab_rollback(tab, snap[level]) < 0)
 374                                 goto error;
 375                         if (tab_add_constraint(tab, map->p[level],
 376                                         div_map[level], index[level][0], 1) < 0)
 377                                 goto error;
 378                         level++;
 379                         continue;
 380                 } else {
 381                         if (k[level] + 1 >= n[level]) {
 382                                 level--;
 383                                 continue;
 384                         }
 385                         if (isl_tab_rollback(tab, snap[level]) < 0)
 386                                 goto error;
 387                         if (tab_add_constraint(tab, map->p[level],
 388                                                 div_map[level],
 389                                                 index[level][k[level]], 0) < 0)
 390                                 goto error;
 391                         snap[level] = isl_tab_snap(tab);
 392                         k[level]++;
 393                         if (tab_add_constraint(tab, map->p[level],
 394                                                 div_map[level],
 395                                                 index[level][k[level]], 1) < 0)
 396                                 goto error;
 397                         level++;
 398                         init = 1;
 399                         continue;
 400                 }
 401         }
 402
 403         isl_tab_free(tab);
 404         free(snap);
 405         free(n);
 406         free(k);
 407         for (i = 0; index && i < map->n; ++i)
 408                 free(index[i]);
 409         free(index);
 410         for (i = 0; div_map && i < map->n; ++i)
 411                 free(div_map[i]);
 412         free(div_map);
 413
 414         isl_basic_map_free(bmap);
 415         isl_map_free(map);
 416
 417         return 0;
 418 error:
 419         isl_tab_free(tab);
 420         free(snap);
 421         free(n);
 422         free(k);
 423         for (i = 0; index && i < map->n; ++i)
 424                 free(index[i]);
 425         free(index);
 426         for (i = 0; div_map && i < map->n; ++i)
 427                 free(div_map[i]);
 428         free(div_map);
 429         isl_basic_map_free(bmap);
 430         isl_map_free(map);
 431         return -1;
 432 }
 433
 434 /* A diff collector that actually collects all parts of the
 435  * set difference in the field diff.
 436  */
 437 struct isl_subtract_diff_collector {
 438         struct isl_diff_collector dc;
 439         struct isl_map *diff;
 440 };
 441
 442 /* isl_subtract_diff_collector callback.
 443  */
 444 static int basic_map_subtract_add(struct isl_diff_collector *dc,
 445                             __isl_take isl_basic_map *bmap)
 446 {
 447         struct isl_subtract_diff_collector *sdc;
 448         sdc = (struct isl_subtract_diff_collector *)dc;
 449
 450         sdc->diff = isl_map_union_disjoint(sdc->diff,
 451                         isl_map_from_basic_map(bmap));
 452
 453         return sdc->diff ? 0 : -1;
 454 }
 455
 456 /* Return the set difference between bmap and map.
 457  */
 458 static __isl_give isl_map *basic_map_subtract(__isl_take isl_basic_map *bmap,
 459         __isl_take isl_map *map)
 460 {
 461         struct isl_subtract_diff_collector sdc;
 462         sdc.dc.add = &basic_map_subtract_add;
 463         sdc.diff = isl_map_empty_like_basic_map(bmap);
 464         if (basic_map_collect_diff(bmap, map, &sdc.dc) < 0) {
 465                 isl_map_free(sdc.diff);
 466                 sdc.diff = NULL;
 467         }
 468         return sdc.diff;
 469 }
 470
 471 /* Return the set difference between map1 and map2.
 472  * (U_i A_i) \ (U_j B_j) is computed as U_i (A_i \ (U_j B_j))
 473  */
 474 struct isl_map *isl_map_subtract(struct isl_map *map1, struct isl_map *map2)
 475 {
 476         int i;
 477         struct isl_map *diff;
 478
 479         if (!map1 || !map2)
 480                 goto error;
 481
 482         isl_assert(map1->ctx, isl_dim_equal(map1->dim, map2->dim), goto error);
 483
 484         if (isl_map_is_empty(map2)) {
 485                 isl_map_free(map2);
 486                 return map1;
 487         }
 488
 489         map1 = isl_map_compute_divs(map1);
 490         map2 = isl_map_compute_divs(map2);
 491         if (!map1 || !map2)
 492                 goto error;
 493
 494         map1 = isl_map_remove_empty_parts(map1);
 495         map2 = isl_map_remove_empty_parts(map2);
 496
 497         diff = isl_map_empty_like(map1);
 498         for (i = 0; i < map1->n; ++i) {
 499                 struct isl_map *d;
 500                 d = basic_map_subtract(isl_basic_map_copy(map1->p[i]),
 501                                        isl_map_copy(map2));
 502                 if (ISL_F_ISSET(map1, ISL_MAP_DISJOINT))
 503                         diff = isl_map_union_disjoint(diff, d);
 504                 else
 505                         diff = isl_map_union(diff, d);
 506         }
 507
 508         isl_map_free(map1);
 509         isl_map_free(map2);
 510
 511         return diff;
 512 error:
 513         isl_map_free(map1);
 514         isl_map_free(map2);
 515         return NULL;
 516 }
 517
 518 struct isl_set *isl_set_subtract(struct isl_set *set1, struct isl_set *set2)
 519 {
 520         return (struct isl_set *)
 521                 isl_map_subtract(
 522                         (struct isl_map *)set1, (struct isl_map *)set2);
 523 }
 524
 525 /* A diff collector that aborts as soon as its add function is called,
 526  * setting empty to 0.
 527  */
 528 struct isl_is_empty_diff_collector {
 529         struct isl_diff_collector dc;
 530         int empty;
 531 };
 532
 533 /* isl_is_empty_diff_collector callback.
 534  */
 535 static int basic_map_is_empty_add(struct isl_diff_collector *dc,
 536                             __isl_take isl_basic_map *bmap)
 537 {
 538         struct isl_is_empty_diff_collector *edc;
 539         edc = (struct isl_is_empty_diff_collector *)dc;
 540
 541         edc->empty = 0;
 542
 543         isl_basic_map_free(bmap);
 544         return -1;
 545 }
 546
 547 /* Check if bmap \ map is empty by computing this set difference
 548  * and breaking off as soon as the difference is known to be non-empty.
 549  */
 550 static int basic_map_diff_is_empty(__isl_keep isl_basic_map *bmap,
 551         __isl_keep isl_map *map)
 552 {
 553         int r;
 554         struct isl_is_empty_diff_collector edc;
 555
 556         r = isl_basic_map_fast_is_empty(bmap);
 557         if (r)
 558                 return r;
 559
 560         edc.dc.add = &basic_map_is_empty_add;
 561         edc.empty = 1;
 562         r = basic_map_collect_diff(isl_basic_map_copy(bmap),
 563                                    isl_map_copy(map), &edc.dc);
 564         if (!edc.empty)
 565                 return 0;
 566
 567         return r < 0 ? -1 : 1;
 568 }
 569
 570 /* Check if map1 \ map2 is empty by checking if the set difference is empty
 571  * for each of the basic maps in map1.
 572  */
 573 static int map_diff_is_empty(__isl_keep isl_map *map1, __isl_keep isl_map *map2)
 574 {
 575         int i;
 576         int is_empty = 1;
 577
 578         if (!map1 || !map2)
 579                 return -1;
 580
 581         for (i = 0; i < map1->n; ++i) {
 582                 is_empty = basic_map_diff_is_empty(map1->p[i], map2);
 583                 if (is_empty < 0 || !is_empty)
 584                          break;
 585         }
 586
 587         return is_empty;
 588 }
 589
 590 /* Return 1 if "bmap" contains a single element.
 591  */
 592 int isl_basic_map_is_singleton(__isl_keep isl_basic_map *bmap)
 593 {
 594         if (!bmap)
 595                 return -1;
 596         if (bmap->n_div)
 597                 return 0;
 598         if (bmap->n_ineq)
 599                 return 0;
 600         return bmap->n_eq == isl_basic_map_total_dim(bmap);
 601 }
 602
 603 /* Return 1 if "map" contains a single element.
 604  */
 605 int isl_map_is_singleton(__isl_keep isl_map *map)
 606 {
 607         if (!map)
 608                 return -1;
 609         if (map->n != 1)
 610                 return 0;
 611
 612         return isl_basic_map_is_singleton(map->p[0]);
 613 }
 614
 615 /* Given a singleton basic map, extract the single element
 616  * as an isl_vec.
 617  */
 618 static __isl_give isl_vec *singleton_extract_point(__isl_keep isl_basic_map *bmap)
 619 {
 620         int i, j;
 621         unsigned dim;
 622         struct isl_vec *point;
 623         isl_int m;
 624
 625         if (!bmap)
 626                 return NULL;
 627
 628         dim = isl_basic_map_total_dim(bmap);
 629         isl_assert(bmap->ctx, bmap->n_eq == dim, return NULL);
 630         point = isl_vec_alloc(bmap->ctx, 1 + dim);
 631         if (!point)
 632                 return NULL;
 633
 634         isl_int_init(m);
 635
 636         isl_int_set_si(point->el[0], 1);
 637         for (j = 0; j < bmap->n_eq; ++j) {
 638                 int s;
 639                 int i = dim - 1 - j;
 640                 isl_assert(bmap->ctx,
 641                     isl_seq_first_non_zero(bmap->eq[j] + 1, i) == -1,
 642                     goto error);
 643                 isl_assert(bmap->ctx,
 644                     isl_int_is_one(bmap->eq[j][1 + i]) ||
 645                     isl_int_is_negone(bmap->eq[j][1 + i]),
 646                     goto error);
 647                 isl_assert(bmap->ctx,
 648                     isl_seq_first_non_zero(bmap->eq[j]+1+i+1, dim-i-1) == -1,
 649                     goto error);
 650
 651                 isl_int_gcd(m, point->el[0], bmap->eq[j][1 + i]);
 652                 isl_int_divexact(m, bmap->eq[j][1 + i], m);
 653                 isl_int_abs(m, m);
 654                 isl_seq_scale(point->el, point->el, m, 1 + i);
 655                 isl_int_divexact(m, point->el[0], bmap->eq[j][1 + i]);
 656                 isl_int_neg(m, m);
 657                 isl_int_mul(point->el[1 + i], m, bmap->eq[j][0]);
 658         }
 659
 660         isl_int_clear(m);
 661         return point;
 662 error:
 663         isl_int_clear(m);
 664         isl_vec_free(point);
 665         return NULL;
 666 }
 667
 668 /* Return 1 if "bmap" contains the point "point".
 669  * "bmap" is assumed to have known divs.
 670  * The point is first extended with the divs and then passed
 671  * to basic_map_contains.
 672  */
 673 static int basic_map_contains_point(__isl_keep isl_basic_map *bmap,
 674         __isl_keep isl_vec *point)
 675 {
 676         int i;
 677         struct isl_vec *vec;
 678         unsigned dim;
 679         int contains;
 680
 681         if (!bmap || !point)
 682                 return -1;
 683         if (bmap->n_div == 0)
 684                 return isl_basic_map_contains(bmap, point);
 685
 686         dim = isl_basic_map_total_dim(bmap) - bmap->n_div;
 687         vec = isl_vec_alloc(bmap->ctx, 1 + dim + bmap->n_div);
 688         if (!vec)
 689                 return -1;
 690
 691         isl_seq_cpy(vec->el, point->el, point->size);
 692         for (i = 0; i < bmap->n_div; ++i) {
 693                 isl_seq_inner_product(bmap->div[i] + 1, vec->el,
 694                                         1 + dim + i, &vec->el[1+dim+i]);
 695                 isl_int_fdiv_q(vec->el[1+dim+i], vec->el[1+dim+i],
 696                                 bmap->div[i][0]);
 697         }
 698
 699         contains = isl_basic_map_contains(bmap, vec);
 700
 701         isl_vec_free(vec);
 702         return contains;
 703 }
 704
 705 /* Return 1 is the singleton map "map1" is a subset of "map2",
 706  * i.e., if the single element of "map1" is also an element of "map2".
 707  */
 708 static int map_is_singleton_subset(__isl_keep isl_map *map1,
 709         __isl_keep isl_map *map2)
 710 {
 711         int i;
 712         int is_subset = 0;
 713         struct isl_vec *point;
 714
 715         if (!map1 || !map2)
 716                 return -1;
 717         if (map1->n != 1)
 718                 return -1;
 719
 720         point = singleton_extract_point(map1->p[0]);
 721         if (!point)
 722                 return -1;
 723
 724         for (i = 0; i < map2->n; ++i) {
 725                 is_subset = basic_map_contains_point(map2->p[i], point);
 726                 if (is_subset)
 727                         break;
 728         }
 729
 730         isl_vec_free(point);
 731         return is_subset;
 732 }
 733
 734 int isl_map_is_subset(struct isl_map *map1, struct isl_map *map2)
 735 {
 736         int is_subset = 0;
 737         struct isl_map *diff;
 738
 739         if (!map1 || !map2)
 740                 return -1;
 741
 742         if (isl_map_is_empty(map1))
 743                 return 1;
 744
 745         if (isl_map_is_empty(map2))
 746                 return 0;
 747
 748         if (isl_map_fast_is_universe(map2))
 749                 return 1;
 750
 751         map1 = isl_map_compute_divs(isl_map_copy(map1));
 752         map2 = isl_map_compute_divs(isl_map_copy(map2));
 753         if (isl_map_is_singleton(map1)) {
 754                 is_subset = map_is_singleton_subset(map1, map2);
 755                 isl_map_free(map1);
 756                 isl_map_free(map2);
 757                 return is_subset;
 758         }
 759         is_subset = map_diff_is_empty(map1, map2);
 760         isl_map_free(map1);
 761         isl_map_free(map2);
 762
 763         return is_subset;
 764 }
 765
 766 int isl_set_is_subset(struct isl_set *set1, struct isl_set *set2)
 767 {
 768         return isl_map_is_subset(
 769                         (struct isl_map *)set1, (struct isl_map *)set2);
 770 }
 771
 772 __isl_give isl_map *isl_map_make_disjoint(__isl_take isl_map *map)
 773 {
 774         int i;
 775         struct isl_subtract_diff_collector sdc;
 776         sdc.dc.add = &basic_map_subtract_add;
 777
 778         if (!map)
 779                 return NULL;
 780         if (ISL_F_ISSET(map, ISL_MAP_DISJOINT))
 781                 return map;
 782         if (map->n <= 1)
 783                 return map;
 784
 785         map = isl_map_compute_divs(map);
 786         map = isl_map_remove_empty_parts(map);
 787
 788         if (!map || map->n <= 1)
 789                 return map;
 790
 791         sdc.diff = isl_map_from_basic_map(isl_basic_map_copy(map->p[0]));
 792
 793         for (i = 1; i < map->n; ++i) {
 794                 struct isl_basic_map *bmap = isl_basic_map_copy(map->p[i]);
 795                 struct isl_map *copy = isl_map_copy(sdc.diff);
 796                 if (basic_map_collect_diff(bmap, copy, &sdc.dc) < 0) {
 797                         isl_map_free(sdc.diff);
 798                         sdc.diff = NULL;
 799                         break;
 800                 }
 801         }
 802
 803         isl_map_free(map);
 804
 805         return sdc.diff;
 806 }
 807
 808 __isl_give isl_set *isl_set_make_disjoint(__isl_take isl_set *set)
 809 {
 810         return (struct isl_set *)isl_map_make_disjoint((struct isl_map *)set);
 811 }
 812
 813 __isl_give isl_set *isl_set_complement(__isl_take isl_set *set)
 814 {
 815         isl_set *universe;
 816
 817         if (!set)
 818                 return NULL;
 819
 820         universe = isl_set_universe(isl_set_get_dim(set));
 821
 822         return isl_set_subtract(universe, set);
 823 }