isl_vertices.c

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