isl_vertices.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_private.h>
  12 #include <isl/set.h>
  13 #include <isl/seq.h>
  14 #include <isl_tab.h>
  15 #include <isl_dim_private.h>
  16 #include <isl_morph.h>
  17 #include <isl_vertices_private.h>
  18 #include <isl_mat_private.h>
  19
  20 #define SELECTED        1
  21 #define DESELECTED      -1
  22 #define UNSELECTED      0
  23
  24 static __isl_give isl_vertices *compute_chambers(__isl_take isl_basic_set *bset,
  25         __isl_take isl_vertices *vertices);
  26
  27 __isl_give isl_vertices *isl_vertices_copy(__isl_keep isl_vertices *vertices)
  28 {
  29         if (!vertices)
  30                 return NULL;
  31
  32         vertices->ref++;
  33         return vertices;
  34 }
  35
  36 void isl_vertices_free(__isl_take isl_vertices *vertices)
  37 {
  38         int i;
  39
  40         if (!vertices)
  41                 return;
  42
  43         if (--vertices->ref > 0)
  44                 return;
  45
  46         for (i = 0; i < vertices->n_vertices; ++i) {
  47                 isl_basic_set_free(vertices->v[i].vertex);
  48                 isl_basic_set_free(vertices->v[i].dom);
  49         }
  50         free(vertices->v);
  51
  52         for (i = 0; i < vertices->n_chambers; ++i) {
  53                 free(vertices->c[i].vertices);
  54                 isl_basic_set_free(vertices->c[i].dom);
  55         }
  56         free(vertices->c);
  57
  58         isl_basic_set_free(vertices->bset);
  59         free(vertices);
  60 }
  61
  62 struct isl_vertex_list {
  63         struct isl_vertex v;
  64         struct isl_vertex_list *next;
  65 };
  66
  67 static void free_vertex_list(struct isl_vertex_list *list)
  68 {
  69         struct isl_vertex_list *next;
  70
  71         for (; list; list = next) {
  72                 next = list->next;
  73                 isl_basic_set_free(list->v.vertex);
  74                 isl_basic_set_free(list->v.dom);
  75                 free(list);
  76         }
  77 }
  78
  79 static __isl_give isl_vertices *vertices_from_list(__isl_keep isl_basic_set *bset,
  80         int n_vertices, struct isl_vertex_list *list)
  81 {
  82         int i;
  83         struct isl_vertex_list *next;
  84         isl_vertices *vertices;
  85
  86         vertices = isl_calloc_type(bset->ctx, isl_vertices);
  87         if (!vertices)
  88                 goto error;
  89         vertices->ref = 1;
  90         vertices->bset = isl_basic_set_copy(bset);
  91         vertices->v = isl_alloc_array(bset->ctx, struct isl_vertex, n_vertices);
  92         if (!vertices->v)
  93                 goto error;
  94         vertices->n_vertices = n_vertices;
  95
  96         for (i = 0; list; list = next, i++) {
  97                 next = list->next;
  98                 vertices->v[i] = list->v;
  99                 free(list);
 100         }
 101
 102         return vertices;
 103 error:
 104         free(vertices);
 105         free_vertex_list(list);
 106         return NULL;
 107 }
 108
 109 /* Prepend a vertex to the linked list "list" based on the equalities in "tab".
 110  */
 111 static int add_vertex(struct isl_vertex_list **list,
 112         __isl_keep isl_basic_set *bset, struct isl_tab *tab)
 113 {
 114         unsigned nvar;
 115         unsigned nparam;
 116         struct isl_vertex_list *v = NULL;
 117
 118         if (isl_tab_detect_implicit_equalities(tab) < 0)
 119                 return -1;
 120
 121         nvar = isl_basic_set_dim(bset, isl_dim_set);
 122         nparam = isl_basic_set_dim(bset, isl_dim_param);
 123
 124         v = isl_calloc_type(tab->mat->ctx, struct isl_vertex_list);
 125         if (!v)
 126                 goto error;
 127
 128         v->v.vertex = isl_basic_set_copy(bset);
 129         v->v.vertex = isl_basic_set_cow(v->v.vertex);
 130         v->v.vertex = isl_basic_set_update_from_tab(v->v.vertex, tab);
 131         v->v.vertex = isl_basic_set_simplify(v->v.vertex);
 132         v->v.vertex = isl_basic_set_finalize(v->v.vertex);
 133         if (!v->v.vertex)
 134                 goto error;
 135         isl_assert(bset->ctx, v->v.vertex->n_eq >= nvar, goto error);
 136         v->v.dom = isl_basic_set_copy(v->v.vertex);
 137         v->v.dom = isl_basic_set_project_out(v->v.dom, isl_dim_set, 0, nvar);
 138         if (!v->v.dom)
 139                 goto error;
 140
 141         v->next = *list;
 142         *list = v;
 143
 144         return 0;
 145 error:
 146         free_vertex_list(v);
 147         return -1;
 148 }
 149
 150 /* Compute the parametric vertices and the chamber decomposition
 151  * of an empty parametric polytope.
 152  */
 153 static __isl_give isl_vertices *vertices_empty(__isl_keep isl_basic_set *bset)
 154 {
 155         isl_vertices *vertices;
 156         unsigned nparam;
 157
 158         if (!bset)
 159                 return NULL;
 160
 161         nparam = isl_basic_set_dim(bset, isl_dim_param);
 162
 163         vertices = isl_calloc_type(bset->ctx, isl_vertices);
 164         if (!vertices)
 165                 return NULL;
 166         vertices->bset = isl_basic_set_copy(bset);
 167         vertices->ref = 1;
 168
 169         vertices->n_vertices = 0;
 170         vertices->n_chambers = 0;
 171
 172         return vertices;
 173 }
 174
 175 /* Compute the parametric vertices and the chamber decomposition
 176  * of the parametric polytope defined using the same constraints
 177  * as "bset" in the 0D case.
 178  * There is exactly one 0D vertex and a single chamber containing
 179  * the vertex.
 180  */
 181 static __isl_give isl_vertices *vertices_0D(__isl_keep isl_basic_set *bset)
 182 {
 183         isl_vertices *vertices;
 184         unsigned nparam;
 185
 186         if (!bset)
 187                 return NULL;
 188
 189         nparam = isl_basic_set_dim(bset, isl_dim_param);
 190
 191         vertices = isl_calloc_type(bset->ctx, isl_vertices);
 192         if (!vertices)
 193                 return NULL;
 194         vertices->ref = 1;
 195         vertices->bset = isl_basic_set_copy(bset);
 196
 197         vertices->v = isl_calloc_array(bset->ctx, struct isl_vertex, 1);
 198         if (!vertices->v)
 199                 goto error;
 200         vertices->n_vertices = 1;
 201         vertices->v[0].vertex = isl_basic_set_copy(bset);
 202         if (!vertices->v[0].vertex)
 203                 goto error;
 204
 205         vertices->c = isl_calloc_array(bset->ctx, struct isl_chamber, 1);
 206         if (!vertices->c)
 207                 goto error;
 208         vertices->n_chambers = 1;
 209         vertices->c[0].n_vertices = 1;
 210         vertices->c[0].vertices = isl_calloc_array(bset->ctx, int, 1);
 211         if (!vertices->c[0].vertices)
 212                 goto error;
 213         vertices->c[0].dom = isl_basic_set_copy(bset);
 214         if (!vertices->c[0].dom)
 215                 goto error;
 216
 217         return vertices;
 218 error:
 219         isl_vertices_free(vertices);
 220         return NULL;
 221 }
 222
 223 static int isl_mat_rank(__isl_keep isl_mat *mat)
 224 {
 225         int row, col;
 226         isl_mat *H;
 227
 228         H = isl_mat_left_hermite(isl_mat_copy(mat), 0, NULL, NULL);
 229         if (!H)
 230                 return -1;
 231
 232         for (col = 0; col < H->n_col; ++col) {
 233                 for (row = 0; row < H->n_row; ++row)
 234                         if (!isl_int_is_zero(H->row[row][col]))
 235                                 break;
 236                 if (row == H->n_row)
 237                         break;
 238         }
 239
 240         isl_mat_free(H);
 241
 242         return col;
 243 }
 244
 245 /* Is the row pointed to by "f" linearly independent of the "n" first
 246  * rows in "facets"?
 247  */
 248 static int is_independent(__isl_keep isl_mat *facets, int n, isl_int *f)
 249 {
 250         int rank;
 251
 252         if (isl_seq_first_non_zero(f, facets->n_col) < 0)
 253                 return 0;
 254
 255         isl_seq_cpy(facets->row[n], f, facets->n_col);
 256         facets->n_row = n + 1;
 257         rank = isl_mat_rank(facets);
 258         if (rank < 0)
 259                 return -1;
 260
 261         return rank == n + 1;
 262 }
 263
 264 /* Check whether we can select constraint "level", given the current selection
 265  * reflected by facets in "tab", the rows of "facets" and the earlier
 266  * "selected" elements of "selection".
 267  *
 268  * If the constraint is (strictly) redundant in the tableau, selecting it would
 269  * result in an empty tableau, so it can't be selected.
 270  * If the set variable part of the constraint is not linearly indepedent
 271  * of the set variable parts of the already selected constraints,
 272  * the constraint cannot be selected.
 273  * If selecting the constraint results in an empty tableau, the constraint
 274  * cannot be selected.
 275  * Finally, if selecting the constraint results in some explicitly
 276  * deselected constraints turning into equalities, then the corresponding
 277  * vertices have already been generated, so the constraint cannot be selected.
 278  */
 279 static int can_select(__isl_keep isl_basic_set *bset, int level,
 280         struct isl_tab *tab, __isl_keep isl_mat *facets, int selected,
 281         int *selection)
 282 {
 283         int i;
 284         int indep;
 285         unsigned ovar;
 286         struct isl_tab_undo *snap;
 287
 288         if (isl_tab_is_redundant(tab, level))
 289                 return 0;
 290
 291         ovar = isl_dim_offset(bset->dim, isl_dim_set);
 292
 293         indep = is_independent(facets, selected, bset->ineq[level] + 1 + ovar);
 294         if (indep < 0)
 295                 return -1;
 296         if (!indep)
 297                 return 0;
 298
 299         snap = isl_tab_snap(tab);
 300         if (isl_tab_select_facet(tab, level) < 0)
 301                 return -1;
 302
 303         if (tab->empty) {
 304                 if (isl_tab_rollback(tab, snap) < 0)
 305                         return -1;
 306                 return 0;
 307         }
 308
 309         for (i = 0; i < level; ++i) {
 310                 int sgn;
 311
 312                 if (selection[i] != DESELECTED)
 313                         continue;
 314
 315                 if (isl_tab_is_equality(tab, i))
 316                         sgn = 0;
 317                 else if (isl_tab_is_redundant(tab, i))
 318                         sgn = 1;
 319                 else
 320                         sgn = isl_tab_sign_of_max(tab, i);
 321                 if (sgn < -1)
 322                         return -1;
 323                 if (sgn <= 0) {
 324                         if (isl_tab_rollback(tab, snap) < 0)
 325                                 return -1;
 326                         return 0;
 327                 }
 328         }
 329
 330         return 1;
 331 }
 332
 333 /* Compute the parametric vertices and the chamber decomposition
 334  * of a parametric polytope that is not full-dimensional.
 335  *
 336  * Simply map the parametric polytope to a lower dimensional space
 337  * and map the resulting vertices back.
 338  */
 339 static __isl_give isl_vertices *lower_dim_vertices(
 340         __isl_keep isl_basic_set *bset)
 341 {
 342         isl_morph *morph;
 343         isl_vertices *vertices;
 344
 345         bset = isl_basic_set_copy(bset);
 346         morph = isl_basic_set_full_compression(bset);
 347         bset = isl_morph_basic_set(isl_morph_copy(morph), bset);
 348
 349         vertices = isl_basic_set_compute_vertices(bset);
 350         isl_basic_set_free(bset);
 351
 352         morph = isl_morph_inverse(morph);
 353
 354         vertices = isl_morph_vertices(morph, vertices);
 355
 356         return vertices;
 357 }
 358
 359 /* Compute the parametric vertices and the chamber decomposition
 360  * of the parametric polytope defined using the same constraints
 361  * as "bset".  "bset" is assumed to have no existentially quantified
 362  * variables.
 363  *
 364  * The vertices themselves are computed in a fairly simplistic way.
 365  * We simply run through all combinations of d constraints,
 366  * with d the number of set variables, and check if those d constraints
 367  * define a vertex.  To avoid the generation of duplicate vertices,
 368  * which we may happen if a vertex is defined by more that d constraints,
 369  * we make sure we only generate the vertex for the d constraints with
 370  * smallest index.
 371  *
 372  * We set up a tableau and keep track of which facets have been
 373  * selected.  The tableau is marked strict_redundant so that we can be
 374  * sure that any constraint that is marked redundant (and that is not
 375  * also marked zero) is not an equality.
 376  * If a constraint is marked DESELECTED, it means the constraint was
 377  * SELECTED before (in combination with the same selection of earlier
 378  * constraints).  If such a deselected constraint turns out to be an
 379  * equality, then any vertex that may still be found with the current
 380  * selection has already been generated when the constraint was selected.
 381  * A constraint is marked UNSELECTED when there is no way selecting
 382  * the constraint could lead to a vertex (in combination with the current
 383  * selection of earlier constraints).
 384  *
 385  * The set variable coefficients of the selected constraints are stored
 386  * in the facets matrix.
 387  */
 388 __isl_give isl_vertices *isl_basic_set_compute_vertices(
 389         __isl_keep isl_basic_set *bset)
 390 {
 391         struct isl_tab *tab;
 392         int level;
 393         int init;
 394         unsigned nvar;
 395         int *selection;
 396         int selected;
 397         struct isl_tab_undo **snap;
 398         isl_mat *facets;
 399         struct isl_vertex_list *list = NULL;
 400         int n_vertices = 0;
 401         isl_vertices *vertices;
 402
 403         if (!bset)
 404                 return NULL;
 405
 406         if (isl_basic_set_fast_is_empty(bset))
 407                 return vertices_empty(bset);
 408
 409         if (bset->n_eq != 0)
 410                 return lower_dim_vertices(bset);
 411
 412         isl_assert(bset->ctx, isl_basic_set_dim(bset, isl_dim_div) == 0,
 413                 return NULL);
 414
 415         if (isl_basic_set_dim(bset, isl_dim_set) == 0)
 416                 return vertices_0D(bset);
 417
 418         nvar = isl_basic_set_dim(bset, isl_dim_set);
 419
 420         bset = isl_basic_set_copy(bset);
 421         bset = isl_basic_set_set_rational(bset);
 422         if (!bset)
 423                 return NULL;
 424
 425         tab = isl_tab_from_basic_set(bset);
 426         if (!tab)
 427                 goto error;
 428         tab->strict_redundant = 1;
 429
 430         if (tab->empty) {
 431                 vertices = vertices_empty(bset);
 432                 isl_basic_set_free(bset);
 433                 isl_tab_free(tab);
 434                 return vertices;
 435         }
 436
 437         selection = isl_alloc_array(bset->ctx, int, bset->n_ineq);
 438         snap = isl_alloc_array(bset->ctx, struct isl_tab_undo *, bset->n_ineq);
 439         facets = isl_mat_alloc(bset->ctx, nvar, nvar);
 440         if (!selection || !snap || !facets)
 441                 goto error;
 442
 443         level = 0;
 444         init = 1;
 445         selected = 0;
 446
 447         while (level >= 0) {
 448                 if (level >= bset->n_ineq ||
 449                     (!init && selection[level] != SELECTED)) {
 450                         --level;
 451                         init = 0;
 452                         continue;
 453                 }
 454                 if (init) {
 455                         int ok;
 456                         snap[level] = isl_tab_snap(tab);
 457                         ok = can_select(bset, level, tab, facets, selected,
 458                                         selection);
 459                         if (ok < 0)
 460                                 goto error;
 461                         if (ok) {
 462                                 selection[level] = SELECTED;
 463                                 selected++;
 464                         } else
 465                                 selection[level] = UNSELECTED;
 466                 } else {
 467                         selection[level] = DESELECTED;
 468                         selected--;
 469                         if (isl_tab_rollback(tab, snap[level]) < 0)
 470                                 goto error;
 471                 }
 472                 if (selected == nvar) {
 473                         if (tab->n_dead == nvar) {
 474                                 if (add_vertex(&list, bset, tab) < 0)
 475                                         goto error;
 476                                 n_vertices++;
 477                         }
 478                         init = 0;
 479                         continue;
 480                 }
 481                 ++level;
 482                 init = 1;
 483         }
 484
 485         isl_mat_free(facets);
 486         free(selection);
 487         free(snap);
 488
 489         isl_tab_free(tab);
 490
 491         vertices = vertices_from_list(bset, n_vertices, list);
 492
 493         vertices = compute_chambers(bset, vertices);
 494
 495         return vertices;
 496 error:
 497         isl_mat_free(facets);
 498         free(selection);
 499         free(snap);
 500         isl_tab_free(tab);
 501         isl_basic_set_free(bset);
 502         return NULL;
 503 }
 504
 505 struct isl_chamber_list {
 506         struct isl_chamber c;
 507         struct isl_chamber_list *next;
 508 };
 509
 510 static void free_chamber_list(struct isl_chamber_list *list)
 511 {
 512         struct isl_chamber_list *next;
 513
 514         for (; list; list = next) {
 515                 next = list->next;
 516                 isl_basic_set_free(list->c.dom);
 517                 free(list->c.vertices);
 518                 free(list);
 519         }
 520 }
 521
 522 /* Check whether the basic set "bset" is a superset of the basic set described
 523  * by "tab", i.e., check whether all constraints of "bset" are redundant.
 524  */
 525 static int bset_covers_tab(__isl_keep isl_basic_set *bset, struct isl_tab *tab)
 526 {
 527         int i;
 528
 529         if (!bset || !tab)
 530                 return -1;
 531
 532         for (i = 0; i < bset->n_ineq; ++i) {
 533                 enum isl_ineq_type type = isl_tab_ineq_type(tab, bset->ineq[i]);
 534                 switch (type) {
 535                 case isl_ineq_error:            return -1;
 536                 case isl_ineq_redundant:        continue;
 537                 default:                        return 0;
 538                 }
 539         }
 540
 541         return 1;
 542 }
 543
 544 static __isl_give isl_vertices *vertices_add_chambers(
 545         __isl_take isl_vertices *vertices, int n_chambers,
 546         struct isl_chamber_list *list)
 547 {
 548         int i;
 549         isl_ctx *ctx;
 550         struct isl_chamber_list *next;
 551
 552         ctx = isl_vertices_get_ctx(vertices);
 553         vertices->c = isl_alloc_array(ctx, struct isl_chamber, n_chambers);
 554         if (!vertices->c)
 555                 goto error;
 556         vertices->n_chambers = n_chambers;
 557
 558         for (i = 0; list; list = next, i++) {
 559                 next = list->next;
 560                 vertices->c[i] = list->c;
 561                 free(list);
 562         }
 563
 564         return vertices;
 565 error:
 566         isl_vertices_free(vertices);
 567         free_chamber_list(list);
 568         return NULL;
 569 }
 570
 571 /* Can "tab" be intersected with "bset" without resulting in
 572  * a lower-dimensional set.
 573  */
 574 static int can_intersect(struct isl_tab *tab, __isl_keep isl_basic_set *bset)
 575 {
 576         int i;
 577         struct isl_tab_undo *snap;
 578
 579         if (isl_tab_extend_cons(tab, bset->n_ineq) < 0)
 580                 return -1;
 581
 582         snap = isl_tab_snap(tab);
 583
 584         for (i = 0; i < bset->n_ineq; ++i) {
 585                 if (isl_tab_ineq_type(tab, bset->ineq[i]) == isl_ineq_redundant)
 586                         continue;
 587                 if (isl_tab_add_ineq(tab, bset->ineq[i]) < 0)
 588                         return -1;
 589         }
 590
 591         if (isl_tab_detect_implicit_equalities(tab) < 0)
 592                 return -1;
 593         if (tab->n_dead) {
 594                 if (isl_tab_rollback(tab, snap) < 0)
 595                         return -1;
 596                 return 0;
 597         }
 598
 599         return 1;
 600 }
 601
 602 static int add_chamber(struct isl_chamber_list **list,
 603         __isl_keep isl_vertices *vertices, struct isl_tab *tab, int *selection)
 604 {
 605         int n_frozen;
 606         int i, j;
 607         int n_vertices = 0;
 608         struct isl_tab_undo *snap;
 609         struct isl_chamber_list *c = NULL;
 610
 611         for (i = 0; i < vertices->n_vertices; ++i)
 612                 if (selection[i])
 613                         n_vertices++;
 614
 615         snap = isl_tab_snap(tab);
 616
 617         for (i = 0; i < tab->n_con && tab->con[i].frozen; ++i)
 618                 tab->con[i].frozen = 0;
 619         n_frozen = i;
 620
 621         if (isl_tab_detect_redundant(tab) < 0)
 622                 return -1;
 623
 624         c = isl_calloc_type(tab->mat->ctx, struct isl_chamber_list);
 625         if (!c)
 626                 goto error;
 627         c->c.vertices = isl_alloc_array(tab->mat->ctx, int, n_vertices);
 628         if (!c->c.vertices)
 629                 goto error;
 630         c->c.dom = isl_basic_set_from_basic_map(isl_basic_map_copy(tab->bmap));
 631         c->c.dom = isl_basic_set_set_rational(c->c.dom);
 632         c->c.dom = isl_basic_set_cow(c->c.dom);
 633         c->c.dom = isl_basic_set_update_from_tab(c->c.dom, tab);
 634         c->c.dom = isl_basic_set_simplify(c->c.dom);
 635         c->c.dom = isl_basic_set_finalize(c->c.dom);
 636         if (!c->c.dom)
 637                 goto error;
 638
 639         c->c.n_vertices = n_vertices;
 640
 641         for (i = 0, j = 0; i < vertices->n_vertices; ++i)
 642                 if (selection[i]) {
 643                         c->c.vertices[j] = i;
 644                         j++;
 645                 }
 646
 647         c->next = *list;
 648         *list = c;
 649
 650         for (i = 0; i < n_frozen; ++i)
 651                 tab->con[i].frozen = 1;
 652
 653         if (isl_tab_rollback(tab, snap) < 0)
 654                 return -1;
 655
 656         return 0;
 657 error:
 658         free_chamber_list(c);
 659         return -1;
 660 }
 661
 662 struct isl_facet_todo {
 663         struct isl_tab *tab;    /* A tableau representation of the facet */
 664         isl_basic_set *bset;    /* A normalized basic set representation */
 665         isl_vec *constraint;    /* Constraint pointing to the other side */
 666         struct isl_facet_todo *next;
 667 };
 668
 669 static void free_todo(struct isl_facet_todo *todo)
 670 {
 671         while (todo) {
 672                 struct isl_facet_todo *next = todo->next;
 673
 674                 isl_tab_free(todo->tab);
 675                 isl_basic_set_free(todo->bset);
 676                 isl_vec_free(todo->constraint);
 677                 free(todo);
 678
 679                 todo = next;
 680         }
 681 }
 682
 683 static struct isl_facet_todo *create_todo(struct isl_tab *tab, int con)
 684 {
 685         int i;
 686         int n_frozen;
 687         struct isl_tab_undo *snap;
 688         struct isl_facet_todo *todo;
 689
 690         snap = isl_tab_snap(tab);
 691
 692         for (i = 0; i < tab->n_con && tab->con[i].frozen; ++i)
 693                 tab->con[i].frozen = 0;
 694         n_frozen = i;
 695
 696         if (isl_tab_detect_redundant(tab) < 0)
 697                 return NULL;
 698
 699         todo = isl_calloc_type(tab->mat->ctx, struct isl_facet_todo);
 700         if (!todo)
 701                 return NULL;
 702
 703         todo->constraint = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
 704         if (!todo->constraint)
 705                 goto error;
 706         isl_seq_neg(todo->constraint->el, tab->bmap->ineq[con], 1 + tab->n_var);
 707         todo->bset = isl_basic_set_from_basic_map(isl_basic_map_copy(tab->bmap));
 708         todo->bset = isl_basic_set_set_rational(todo->bset);
 709         todo->bset = isl_basic_set_cow(todo->bset);
 710         todo->bset = isl_basic_set_update_from_tab(todo->bset, tab);
 711         todo->bset = isl_basic_set_simplify(todo->bset);
 712         todo->bset = isl_basic_set_sort_constraints(todo->bset);
 713         if (!todo->bset)
 714                 goto error;
 715         ISL_F_SET(todo->bset, ISL_BASIC_SET_NORMALIZED);
 716         todo->tab = isl_tab_dup(tab);
 717         if (!todo->tab)
 718                 goto error;
 719
 720         for (i = 0; i < n_frozen; ++i)
 721                 tab->con[i].frozen = 1;
 722
 723         if (isl_tab_rollback(tab, snap) < 0)
 724                 goto error;
 725
 726         return todo;
 727 error:
 728         free_todo(todo);
 729         return NULL;
 730 }
 731
 732 /* Create todo items for all interior facets of the chamber represented
 733  * by "tab" and collect them in "next".
 734  */
 735 static int init_todo(struct isl_facet_todo **next, struct isl_tab *tab)
 736 {
 737         int i;
 738         struct isl_tab_undo *snap;
 739         struct isl_facet_todo *todo;
 740
 741         snap = isl_tab_snap(tab);
 742
 743         for (i = 0; i < tab->n_con; ++i) {
 744                 if (tab->con[i].frozen)
 745                         continue;
 746                 if (tab->con[i].is_redundant)
 747                         continue;
 748
 749                 if (isl_tab_select_facet(tab, i) < 0)
 750                         return -1;
 751
 752                 todo = create_todo(tab, i);
 753                 if (!todo)
 754                         return -1;
 755
 756                 todo->next = *next;
 757                 *next = todo;
 758
 759                 if (isl_tab_rollback(tab, snap) < 0)
 760                         return -1;
 761         }
 762
 763         return 0;
 764 }
 765
 766 /* Does the linked list contain a todo item that is the opposite of "todo".
 767  * If so, return 1 and remove the opposite todo item.
 768  */
 769 static int has_opposite(struct isl_facet_todo *todo,
 770         struct isl_facet_todo **list)
 771 {
 772         for (; *list; list = &(*list)->next) {
 773                 int eq;
 774                 eq = isl_basic_set_fast_is_equal(todo->bset, (*list)->bset);
 775                 if (eq < 0)
 776                         return -1;
 777                 if (!eq)
 778                         continue;
 779                 todo = *list;
 780                 *list = todo->next;
 781                 todo->next = NULL;
 782                 free_todo(todo);
 783                 return 1;
 784         }
 785
 786         return 0;
 787 }
 788
 789 /* Create todo items for all interior facets of the chamber represented
 790  * by "tab" and collect them in first->next, taking care to cancel
 791  * opposite todo items.
 792  */
 793 static int update_todo(struct isl_facet_todo *first, struct isl_tab *tab)
 794 {
 795         int i;
 796         struct isl_tab_undo *snap;
 797         struct isl_facet_todo *todo;
 798
 799         snap = isl_tab_snap(tab);
 800
 801         for (i = 0; i < tab->n_con; ++i) {
 802                 int drop;
 803
 804                 if (tab->con[i].frozen)
 805                         continue;
 806                 if (tab->con[i].is_redundant)
 807                         continue;
 808
 809                 if (isl_tab_select_facet(tab, i) < 0)
 810                         return -1;
 811
 812                 todo = create_todo(tab, i);
 813                 if (!todo)
 814                         return -1;
 815
 816                 drop = has_opposite(todo, &first->next);
 817                 if (drop < 0)
 818                         return -1;
 819
 820                 if (drop)
 821                         free_todo(todo);
 822                 else {
 823                         todo->next = first->next;
 824                         first->next = todo;
 825                 }
 826
 827                 if (isl_tab_rollback(tab, snap) < 0)
 828                         return -1;
 829         }
 830
 831         return 0;
 832 }
 833
 834 /* Compute the chamber decomposition of the parametric polytope respresented
 835  * by "bset" given the parametric vertices and their activity domains.
 836  *
 837  * We are only interested in full-dimensional chambers.
 838  * Each of these chambers is the intersection of the activity domains of
 839  * one or more vertices and the union of all chambers is equal to the
 840  * projection of the entire parametric polytope onto the parameter space.
 841  *
 842  * We first create an initial chamber by intersecting as many activity
 843  * domains as possible without ending up with an empty or lower-dimensional
 844  * set.  As a minor optimization, we only consider those activity domains
 845  * that contain some arbitrary point.
 846  *
 847  * For each of interior facets of the chamber, we construct a todo item,
 848  * containing the facet and a constraint containing the other side of the facet,
 849  * for constructing the chamber on the other side.
 850  * While their are any todo items left, we pick a todo item and
 851  * create the required chamber by intersecting all activity domains
 852  * that contain the facet and have a full-dimensional intersection with
 853  * the other side of the facet.  For each of the interior facets, we
 854  * again create todo items, taking care to cancel opposite todo items.
 855  */
 856 static __isl_give isl_vertices *compute_chambers(__isl_take isl_basic_set *bset,
 857         __isl_take isl_vertices *vertices)
 858 {
 859         int i;
 860         isl_ctx *ctx;
 861         isl_vec *sample = NULL;
 862         struct isl_tab *tab = NULL;
 863         struct isl_tab_undo *snap;
 864         unsigned nvar;
 865         int *selection = NULL;
 866         int n_chambers = 0;
 867         struct isl_chamber_list *list = NULL;
 868         struct isl_facet_todo *todo = NULL;
 869
 870         if (!bset || !vertices)
 871                 goto error;
 872
 873         ctx = isl_vertices_get_ctx(vertices);
 874         selection = isl_alloc_array(ctx, int, vertices->n_vertices);
 875         if (!selection)
 876                 goto error;
 877
 878         nvar = isl_basic_set_dim(bset, isl_dim_set);
 879         bset = isl_basic_set_project_out(bset, isl_dim_set, 0, nvar);
 880
 881         tab = isl_tab_from_basic_set(bset);
 882         for (i = 0; i < bset->n_ineq; ++i)
 883                 if (isl_tab_freeze_constraint(tab, i) < 0)
 884                         goto error;
 885         if (isl_tab_track_bset(tab, bset) < 0)
 886                 goto error;
 887
 888         snap = isl_tab_snap(tab);
 889
 890         sample = isl_tab_get_sample_value(tab);
 891
 892         for (i = 0; i < vertices->n_vertices; ++i) {
 893                 selection[i] = isl_basic_set_contains(vertices->v[i].dom, sample);
 894                 if (selection[i] < 0)
 895                         goto error;
 896                 if (!selection[i])
 897                         continue;
 898                 selection[i] = can_intersect(tab, vertices->v[i].dom);
 899                 if (selection[i] < 0)
 900                         goto error;
 901         }
 902
 903         if (isl_tab_detect_redundant(tab) < 0)
 904                 goto error;
 905
 906         if (add_chamber(&list, vertices, tab, selection) < 0)
 907                 goto error;
 908         n_chambers++;
 909
 910         if (init_todo(&todo, tab) < 0)
 911                 goto error;
 912
 913         while (todo) {
 914                 struct isl_facet_todo *next;
 915
 916                 if (isl_tab_rollback(tab, snap) < 0)
 917                         goto error;
 918
 919                 if (isl_tab_add_ineq(tab, todo->constraint->el) < 0)
 920                         goto error;
 921                 if (isl_tab_freeze_constraint(tab, tab->n_con - 1) < 0)
 922                         goto error;
 923
 924                 for (i = 0; i < vertices->n_vertices; ++i) {
 925                         selection[i] = bset_covers_tab(vertices->v[i].dom,
 926                                                         todo->tab);
 927                         if (selection[i] < 0)
 928                                 goto error;
 929                         if (!selection[i])
 930                                 continue;
 931                         selection[i] = can_intersect(tab, vertices->v[i].dom);
 932                         if (selection[i] < 0)
 933                                 goto error;
 934                 }
 935
 936                 if (isl_tab_detect_redundant(tab) < 0)
 937                         goto error;
 938
 939                 if (add_chamber(&list, vertices, tab, selection) < 0)
 940                         goto error;
 941                 n_chambers++;
 942
 943                 if (update_todo(todo, tab) < 0)
 944                         goto error;
 945
 946                 next = todo->next;
 947                 todo->next = NULL;
 948                 free_todo(todo);
 949                 todo = next;
 950         }
 951
 952         isl_vec_free(sample);
 953
 954         isl_tab_free(tab);
 955         free(selection);
 956
 957         vertices = vertices_add_chambers(vertices, n_chambers, list);
 958
 959         for (i = 0; vertices && i < vertices->n_vertices; ++i) {
 960                 isl_basic_set_free(vertices->v[i].dom);
 961                 vertices->v[i].dom = NULL;
 962         }
 963
 964         return vertices;
 965 error:
 966         free_chamber_list(list);
 967         free_todo(todo);
 968         isl_vec_free(sample);
 969         isl_tab_free(tab);
 970         free(selection);
 971         if (!tab)
 972                 isl_basic_set_free(bset);
 973         isl_vertices_free(vertices);
 974         return NULL;
 975 }
 976
 977 isl_ctx *isl_vertex_get_ctx(__isl_keep isl_vertex *vertex)
 978 {
 979         return vertex ? isl_vertices_get_ctx(vertex->vertices) : NULL;
 980 }
 981
 982 int isl_vertex_get_id(__isl_keep isl_vertex *vertex)
 983 {
 984         return vertex ? vertex->id : -1;
 985 }
 986
 987 __isl_give isl_basic_set *isl_vertex_get_domain(__isl_keep isl_vertex *vertex)
 988 {
 989         struct isl_vertex *v;
 990
 991         if (!vertex)
 992                 return NULL;
 993
 994         v = &vertex->vertices->v[vertex->id];
 995         if (!v->dom) {
 996                 unsigned nvar;
 997                 nvar = isl_basic_set_dim(v->vertex, isl_dim_set);
 998                 v->dom = isl_basic_set_copy(v->vertex);
 999                 v->dom = isl_basic_set_project_out(v->dom, isl_dim_set, 0, nvar);
1000         }
1001
1002         return isl_basic_set_copy(v->dom);
1003 }
1004
1005 __isl_give isl_basic_set *isl_vertex_get_expr(__isl_keep isl_vertex *vertex)
1006 {
1007         struct isl_vertex *v;
1008
1009         if (!vertex)
1010                 return NULL;
1011
1012         v = &vertex->vertices->v[vertex->id];
1013
1014         return isl_basic_set_copy(v->vertex);
1015 }
1016
1017 static __isl_give isl_vertex *isl_vertex_alloc(__isl_take isl_vertices *vertices,
1018         int id)
1019 {
1020         isl_ctx *ctx;
1021         isl_vertex *vertex;
1022
1023         if (!vertices)
1024                 return NULL;
1025
1026         ctx = isl_vertices_get_ctx(vertices);
1027         vertex = isl_alloc_type(ctx, isl_vertex);
1028         if (!vertex)
1029                 goto error;
1030
1031         vertex->vertices = vertices;
1032         vertex->id = id;
1033
1034         return vertex;
1035 error:
1036         isl_vertices_free(vertices);
1037         return NULL;
1038 }
1039
1040 void isl_vertex_free(__isl_take isl_vertex *vertex)
1041 {
1042         if (!vertex)
1043                 return;
1044         isl_vertices_free(vertex->vertices);
1045         free(vertex);
1046 }
1047
1048 __isl_give isl_basic_set *isl_basic_set_set_integral(__isl_take isl_basic_set *bset)
1049 {
1050         if (!bset)
1051                 return NULL;
1052
1053         if (!ISL_F_ISSET(bset, ISL_BASIC_MAP_RATIONAL))
1054                 return bset;
1055
1056         bset = isl_basic_set_cow(bset);
1057         if (!bset)
1058                 return NULL;
1059
1060         ISL_F_CLR(bset, ISL_BASIC_MAP_RATIONAL);
1061
1062         return isl_basic_set_finalize(bset);
1063 }
1064
1065 isl_ctx *isl_cell_get_ctx(__isl_keep isl_cell *cell)
1066 {
1067         return cell ? cell->dom->ctx : NULL;
1068 }
1069
1070 __isl_give isl_basic_set *isl_cell_get_domain(__isl_keep isl_cell *cell)
1071 {
1072         return cell ? isl_basic_set_copy(cell->dom) : NULL;
1073 }
1074
1075 static __isl_give isl_cell *isl_cell_alloc(__isl_take isl_vertices *vertices,
1076         __isl_take isl_basic_set *dom, int id)
1077 {
1078         int i;
1079         isl_cell *cell = NULL;
1080
1081         if (!vertices || !dom)
1082                 goto error;
1083
1084         cell = isl_calloc_type(dom->ctx, isl_cell);
1085         if (!cell)
1086                 goto error;
1087
1088         cell->n_vertices = vertices->c[id].n_vertices;
1089         cell->ids = isl_alloc_array(dom->ctx, int, cell->n_vertices);
1090         if (!cell->ids)
1091                 goto error;
1092         for (i = 0; i < cell->n_vertices; ++i)
1093                 cell->ids[i] = vertices->c[id].vertices[i];
1094         cell->vertices = vertices;
1095         cell->dom = dom;
1096
1097         return cell;
1098 error:
1099         isl_cell_free(cell);
1100         isl_vertices_free(vertices);
1101         isl_basic_set_free(dom);
1102         return NULL;
1103 }
1104
1105 void isl_cell_free(__isl_take isl_cell *cell)
1106 {
1107         if (!cell)
1108                 return;
1109
1110         isl_vertices_free(cell->vertices);
1111         free(cell->ids);
1112         isl_basic_set_free(cell->dom);
1113         free(cell);
1114 }
1115
1116 /* Create a tableau of the cone obtained by first homogenizing the given
1117  * polytope and then making all inequalities strict by setting the
1118  * constant term to -1.
1119  */
1120 static struct isl_tab *tab_for_shifted_cone(__isl_keep isl_basic_set *bset)
1121 {
1122         int i;
1123         isl_vec *c = NULL;
1124         struct isl_tab *tab;
1125
1126         if (!bset)
1127                 return NULL;
1128         tab = isl_tab_alloc(bset->ctx, bset->n_ineq + 1,
1129                             1 + isl_basic_set_total_dim(bset), 0);
1130         if (!tab)
1131                 return NULL;
1132         tab->rational = ISL_F_ISSET(bset, ISL_BASIC_SET_RATIONAL);
1133         if (ISL_F_ISSET(bset, ISL_BASIC_MAP_EMPTY)) {
1134                 if (isl_tab_mark_empty(tab) < 0)
1135                         goto error;
1136                 return tab;
1137         }
1138
1139         c = isl_vec_alloc(bset->ctx, 1 + 1 + isl_basic_set_total_dim(bset));
1140         if (!c)
1141                 goto error;
1142
1143         isl_int_set_si(c->el[0], 0);
1144         for (i = 0; i < bset->n_eq; ++i) {
1145                 isl_seq_cpy(c->el + 1, bset->eq[i], c->size - 1);
1146                 if (isl_tab_add_eq(tab, c->el) < 0)
1147                         goto error;
1148         }
1149
1150         isl_int_set_si(c->el[0], -1);
1151         for (i = 0; i < bset->n_ineq; ++i) {
1152                 isl_seq_cpy(c->el + 1, bset->ineq[i], c->size - 1);
1153                 if (isl_tab_add_ineq(tab, c->el) < 0)
1154                         goto error;
1155                 if (tab->empty) {
1156                         isl_vec_free(c);
1157                         return tab;
1158                 }
1159         }
1160
1161         isl_seq_clr(c->el + 1, c->size - 1);
1162         isl_int_set_si(c->el[1], 1);
1163         if (isl_tab_add_ineq(tab, c->el) < 0)
1164                 goto error;
1165
1166         isl_vec_free(c);
1167         return tab;
1168 error:
1169         isl_vec_free(c);
1170         isl_tab_free(tab);
1171         return NULL;
1172 }
1173
1174 /* Compute an interior point of "bset" by selecting an interior
1175  * point in homogeneous space and projecting the point back down.
1176  */
1177 static __isl_give isl_vec *isl_basic_set_interior_point(
1178         __isl_keep isl_basic_set *bset)
1179 {
1180         isl_vec *vec;
1181         struct isl_tab *tab;
1182
1183         tab = tab_for_shifted_cone(bset);
1184         vec = isl_tab_get_sample_value(tab);
1185         isl_tab_free(tab);
1186         if (!vec)
1187                 return NULL;
1188
1189         isl_seq_cpy(vec->el, vec->el + 1, vec->size - 1);
1190         vec->size--;
1191
1192         return vec;
1193 }
1194
1195 /* Call "fn" on all chambers of the parametric polytope with the shared
1196  * facets of neighboring chambers only appearing in one of the chambers.
1197  *
1198  * We pick an interior point from one of the chambers and then make
1199  * all constraints that do not satisfy this point strict.
1200  */
1201 int isl_vertices_foreach_disjoint_cell(__isl_keep isl_vertices *vertices,
1202         int (*fn)(__isl_take isl_cell *cell, void *user), void *user)
1203 {
1204         int i, j;
1205         isl_vec *vec;
1206         isl_int v;
1207         isl_cell *cell;
1208
1209         if (!vertices)
1210                 return -1;
1211
1212         if (vertices->n_chambers == 0)
1213                 return 0;
1214
1215         if (vertices->n_chambers == 1) {
1216                 isl_basic_set *dom = isl_basic_set_copy(vertices->c[0].dom);
1217                 dom = isl_basic_set_set_integral(dom);
1218                 cell = isl_cell_alloc(isl_vertices_copy(vertices), dom, 0);
1219                 if (!cell)
1220                         return -1;
1221                 return fn(cell, user);
1222         }
1223
1224         vec = isl_basic_set_interior_point(vertices->c[0].dom);
1225         if (!vec)
1226                 return -1;
1227
1228         isl_int_init(v);
1229
1230         for (i = 0; i < vertices->n_chambers; ++i) {
1231                 int r;
1232                 isl_basic_set *dom = isl_basic_set_copy(vertices->c[i].dom);
1233                 dom = isl_basic_set_cow(dom);
1234                 if (!dom)
1235                         goto error;
1236                 for (j = 0; i && j < dom->n_ineq; ++j) {
1237                         isl_seq_inner_product(vec->el, dom->ineq[j], vec->size,
1238                                                 &v);
1239                         if (!isl_int_is_neg(v))
1240                                 continue;
1241                         isl_int_sub_ui(dom->ineq[j][0], dom->ineq[j][0], 1);
1242                 }
1243                 dom = isl_basic_set_set_integral(dom);
1244                 cell = isl_cell_alloc(isl_vertices_copy(vertices), dom, i);
1245                 if (!cell)
1246                         goto error;
1247                 r = fn(cell, user);
1248                 if (r < 0)
1249                         goto error;
1250         }
1251
1252         isl_int_clear(v);
1253         isl_vec_free(vec);
1254
1255         return 0;
1256 error:
1257         isl_int_clear(v);
1258         isl_vec_free(vec);
1259         return -1;
1260 }
1261
1262 int isl_vertices_foreach_cell(__isl_keep isl_vertices *vertices,
1263         int (*fn)(__isl_take isl_cell *cell, void *user), void *user)
1264 {
1265         int i;
1266         isl_cell *cell;
1267
1268         if (!vertices)
1269                 return -1;
1270
1271         if (vertices->n_chambers == 0)
1272                 return 0;
1273
1274         for (i = 0; i < vertices->n_chambers; ++i) {
1275                 int r;
1276                 isl_basic_set *dom = isl_basic_set_copy(vertices->c[i].dom);
1277
1278                 cell = isl_cell_alloc(isl_vertices_copy(vertices), dom, i);
1279                 if (!cell)
1280                         return -1;
1281
1282                 r = fn(cell, user);
1283                 if (r < 0)
1284                         return -1;
1285         }
1286
1287         return 0;
1288 }
1289
1290 int isl_vertices_foreach_vertex(__isl_keep isl_vertices *vertices,
1291         int (*fn)(__isl_take isl_vertex *vertex, void *user), void *user)
1292 {
1293         int i;
1294         isl_vertex *vertex;
1295
1296         if (!vertices)
1297                 return -1;
1298
1299         if (vertices->n_vertices == 0)
1300                 return 0;
1301
1302         for (i = 0; i < vertices->n_vertices; ++i) {
1303                 int r;
1304
1305                 vertex = isl_vertex_alloc(isl_vertices_copy(vertices), i);
1306                 if (!vertex)
1307                         return -1;
1308
1309                 r = fn(vertex, user);
1310                 if (r < 0)
1311                         return -1;
1312         }
1313
1314         return 0;
1315 }
1316
1317 int isl_cell_foreach_vertex(__isl_keep isl_cell *cell,
1318         int (*fn)(__isl_take isl_vertex *vertex, void *user), void *user)
1319 {
1320         int i;
1321         isl_vertex *vertex;
1322
1323         if (!cell)
1324                 return -1;
1325
1326         if (cell->n_vertices == 0)
1327                 return 0;
1328
1329         for (i = 0; i < cell->n_vertices; ++i) {
1330                 int r;
1331
1332                 vertex = isl_vertex_alloc(isl_vertices_copy(cell->vertices),
1333                                           cell->ids[i]);
1334                 if (!vertex)
1335                         return -1;
1336
1337                 r = fn(vertex, user);
1338                 if (r < 0)
1339                         return -1;
1340         }
1341
1342         return 0;
1343 }
1344
1345 isl_ctx *isl_vertices_get_ctx(__isl_keep isl_vertices *vertices)
1346 {
1347         return vertices ? vertices->bset->ctx : NULL;
1348 }
1349
1350 int isl_vertices_get_n_vertices(__isl_keep isl_vertices *vertices)
1351 {
1352         return vertices ? vertices->n_vertices : -1;
1353 }
1354
1355 __isl_give isl_vertices *isl_morph_vertices(__isl_take isl_morph *morph,
1356         __isl_take isl_vertices *vertices)
1357 {
1358         int i;
1359         isl_morph *param_morph = NULL;
1360
1361         if (!morph || !vertices)
1362                 goto error;
1363
1364         isl_assert(vertices->bset->ctx, vertices->ref == 1, goto error);
1365
1366         param_morph = isl_morph_copy(morph);
1367         param_morph = isl_morph_remove_dom_dims(param_morph, isl_dim_set,
1368                                     0, isl_morph_dom_dim(morph, isl_dim_set));
1369         param_morph = isl_morph_remove_ran_dims(param_morph, isl_dim_set,
1370                                     0, isl_morph_ran_dim(morph, isl_dim_set));
1371
1372         for (i = 0; i < vertices->n_vertices; ++i) {
1373                 vertices->v[i].dom = isl_morph_basic_set(
1374                         isl_morph_copy(param_morph), vertices->v[i].dom);
1375                 vertices->v[i].vertex = isl_morph_basic_set(
1376                         isl_morph_copy(morph), vertices->v[i].vertex);
1377                 if (!vertices->v[i].vertex)
1378                         goto error;
1379         }
1380
1381         for (i = 0; i < vertices->n_chambers; ++i) {
1382                 vertices->c[i].dom = isl_morph_basic_set(
1383                         isl_morph_copy(param_morph), vertices->c[i].dom);
1384                 if (!vertices->c[i].dom)
1385                         goto error;
1386         }
1387
1388         isl_morph_free(param_morph);
1389         isl_morph_free(morph);
1390         return vertices;
1391 error:
1392         isl_morph_free(param_morph);
1393         isl_morph_free(morph);
1394         isl_vertices_free(vertices);
1395         return NULL;
1396 }
1397
1398 /* Construct a simplex isl_cell spanned by the vertices with indices in
1399  * "simplex_ids" and "other_ids" and call "fn" on this isl_cell.
1400  */
1401 static int call_on_simplex(__isl_keep isl_cell *cell,
1402         int *simplex_ids, int n_simplex, int *other_ids, int n_other,
1403         int (*fn)(__isl_take isl_cell *simplex, void *user), void *user)
1404 {
1405         int i;
1406         isl_ctx *ctx;
1407         struct isl_cell *simplex;
1408
1409         ctx = isl_cell_get_ctx(cell);
1410
1411         simplex = isl_calloc_type(ctx, struct isl_cell);
1412         if (!simplex)
1413                 return -1;
1414         simplex->vertices = isl_vertices_copy(cell->vertices);
1415         if (!simplex->vertices)
1416                 goto error;
1417         simplex->dom = isl_basic_set_copy(cell->dom);
1418         if (!simplex->dom)
1419                 goto error;
1420         simplex->n_vertices = n_simplex + n_other;
1421         simplex->ids = isl_alloc_array(ctx, int, simplex->n_vertices);
1422         if (!simplex->ids)
1423                 goto error;
1424
1425         for (i = 0; i < n_simplex; ++i)
1426                 simplex->ids[i] = simplex_ids[i];
1427         for (i = 0; i < n_other; ++i)
1428                 simplex->ids[n_simplex + i] = other_ids[i];
1429
1430         return fn(simplex, user);
1431 error:
1432         isl_cell_free(simplex);
1433         return -1;
1434 }
1435
1436 /* Check whether the parametric vertex described by "vertex"
1437  * lies on the facet corresponding to constraint "facet" of "bset".
1438  * The isl_vec "v" is a temporary vector than can be used by this function.
1439  *
1440  * We eliminate the variables from the facet constraint using the
1441  * equalities defining the vertex and check if the result is identical
1442  * to zero.
1443  *
1444  * It would probably be better to keep track of the constraints defining
1445  * a vertex during the vertex construction so that we could simply look
1446  * it up here.
1447  */
1448 static int vertex_on_facet(__isl_keep isl_basic_set *vertex,
1449         __isl_keep isl_basic_set *bset, int facet, __isl_keep isl_vec *v)
1450 {
1451         int i;
1452         isl_int m;
1453
1454         isl_seq_cpy(v->el, bset->ineq[facet], v->size);
1455
1456         isl_int_init(m);
1457         for (i = 0; i < vertex->n_eq; ++i) {
1458                 int k = isl_seq_last_non_zero(vertex->eq[i], v->size);
1459                 isl_seq_elim(v->el, vertex->eq[i], k, v->size, &m);
1460         }
1461         isl_int_clear(m);
1462
1463         return isl_seq_first_non_zero(v->el, v->size) == -1;
1464 }
1465
1466 /* Triangulate the polytope spanned by the vertices with ids
1467  * in "simplex_ids" and "other_ids" and call "fn" on each of
1468  * the resulting simplices.
1469  * If the input polytope is already a simplex, we simply call "fn".
1470  * Otherwise, we pick a point from "other_ids" and add it to "simplex_ids".
1471  * Then we consider each facet of "bset" that does not contain the point
1472  * we just picked, but does contain some of the other points in "other_ids"
1473  * and call ourselves recursively on the polytope spanned by the new
1474  * "simplex_ids" and those points in "other_ids" that lie on the facet.
1475  */
1476 static int triangulate(__isl_keep isl_cell *cell, __isl_keep isl_vec *v,
1477         int *simplex_ids, int n_simplex, int *other_ids, int n_other,
1478         int (*fn)(__isl_take isl_cell *simplex, void *user), void *user)
1479 {
1480         int i, j, k;
1481         int d, nparam;
1482         int *ids;
1483         isl_ctx *ctx;
1484         isl_basic_set *vertex;
1485         isl_basic_set *bset;
1486
1487         ctx = isl_cell_get_ctx(cell);
1488         d = isl_basic_set_dim(cell->vertices->bset, isl_dim_set);
1489         nparam = isl_basic_set_dim(cell->vertices->bset, isl_dim_param);
1490
1491         if (n_simplex + n_other == d + 1)
1492                 return call_on_simplex(cell, simplex_ids, n_simplex,
1493                                        other_ids, n_other, fn, user);
1494
1495         simplex_ids[n_simplex] = other_ids[0];
1496         vertex = cell->vertices->v[other_ids[0]].vertex;
1497         bset = cell->vertices->bset;
1498
1499         ids = isl_alloc_array(ctx, int, n_other - 1);
1500         for (i = 0; i < bset->n_ineq; ++i) {
1501                 if (isl_seq_first_non_zero(bset->ineq[i] + 1 + nparam, d) == -1)
1502                         continue;
1503                 if (vertex_on_facet(vertex, bset, i, v))
1504                         continue;
1505
1506                 for (j = 1, k = 0; j < n_other; ++j) {
1507                         isl_basic_set *ov;
1508                         ov = cell->vertices->v[other_ids[j]].vertex;
1509                         if (vertex_on_facet(ov, bset, i, v))
1510                                 ids[k++] = other_ids[j];
1511                 }
1512                 if (k == 0)
1513                         continue;
1514
1515                 if (triangulate(cell, v, simplex_ids, n_simplex + 1,
1516                                 ids, k, fn, user) < 0)
1517                         goto error;
1518         }
1519         free(ids);
1520
1521         return 0;
1522 error:
1523         free(ids);
1524         return -1;
1525 }
1526
1527 /* Triangulate the given cell and call "fn" on each of the resulting
1528  * simplices.
1529  */
1530 int isl_cell_foreach_simplex(__isl_take isl_cell *cell,
1531         int (*fn)(__isl_take isl_cell *simplex, void *user), void *user)
1532 {
1533         int d, total;
1534         int r;
1535         isl_ctx *ctx;
1536         isl_vec *v = NULL;
1537         int *simplex_ids = NULL;
1538
1539         if (!cell)
1540                 return -1;
1541
1542         d = isl_basic_set_dim(cell->vertices->bset, isl_dim_set);
1543         total = isl_basic_set_total_dim(cell->vertices->bset);
1544
1545         if (cell->n_vertices == d + 1)
1546                 return fn(cell, user);
1547
1548         ctx = isl_cell_get_ctx(cell);
1549         simplex_ids = isl_alloc_array(ctx, int, d + 1);
1550         if (!simplex_ids)
1551                 goto error;
1552
1553         v = isl_vec_alloc(ctx, 1 + total);
1554         if (!v)
1555                 goto error;
1556
1557         r = triangulate(cell, v, simplex_ids, 0,
1558                         cell->ids, cell->n_vertices, fn, user);
1559
1560         isl_vec_free(v);
1561         free(simplex_ids);
1562
1563         isl_cell_free(cell);
1564
1565         return r;
1566 error:
1567         free(simplex_ids);
1568         isl_vec_free(v);
1569         isl_cell_free(cell);
1570         return -1;
1571 }