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