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