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