add packaging
[platform/upstream/isl.git] / 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         if (!tab)
883                 goto error;
884         for (i = 0; i < bset->n_ineq; ++i)
885                 if (isl_tab_freeze_constraint(tab, i) < 0)
886                         goto error;
887         isl_basic_set_free(bset);
888
889         snap = isl_tab_snap(tab);
890
891         sample = isl_tab_get_sample_value(tab);
892
893         for (i = 0; i < vertices->n_vertices; ++i) {
894                 selection[i] = isl_basic_set_contains(vertices->v[i].dom, sample);
895                 if (selection[i] < 0)
896                         goto error;
897                 if (!selection[i])
898                         continue;
899                 selection[i] = can_intersect(tab, vertices->v[i].dom);
900                 if (selection[i] < 0)
901                         goto error;
902         }
903
904         if (isl_tab_detect_redundant(tab) < 0)
905                 goto error;
906
907         if (add_chamber(&list, vertices, tab, selection) < 0)
908                 goto error;
909         n_chambers++;
910
911         if (init_todo(&todo, tab) < 0)
912                 goto error;
913
914         while (todo) {
915                 struct isl_facet_todo *next;
916
917                 if (isl_tab_rollback(tab, snap) < 0)
918                         goto error;
919
920                 if (isl_tab_add_ineq(tab, todo->constraint->el) < 0)
921                         goto error;
922                 if (isl_tab_freeze_constraint(tab, tab->n_con - 1) < 0)
923                         goto error;
924
925                 for (i = 0; i < vertices->n_vertices; ++i) {
926                         selection[i] = bset_covers_tab(vertices->v[i].dom,
927                                                         todo->tab);
928                         if (selection[i] < 0)
929                                 goto error;
930                         if (!selection[i])
931                                 continue;
932                         selection[i] = can_intersect(tab, vertices->v[i].dom);
933                         if (selection[i] < 0)
934                                 goto error;
935                 }
936
937                 if (isl_tab_detect_redundant(tab) < 0)
938                         goto error;
939
940                 if (add_chamber(&list, vertices, tab, selection) < 0)
941                         goto error;
942                 n_chambers++;
943
944                 if (update_todo(todo, tab) < 0)
945                         goto error;
946
947                 next = todo->next;
948                 todo->next = NULL;
949                 free_todo(todo);
950                 todo = next;
951         }
952
953         isl_vec_free(sample);
954
955         isl_tab_free(tab);
956         free(selection);
957
958         vertices = vertices_add_chambers(vertices, n_chambers, list);
959
960         for (i = 0; vertices && i < vertices->n_vertices; ++i) {
961                 isl_basic_set_free(vertices->v[i].dom);
962                 vertices->v[i].dom = NULL;
963         }
964
965         return vertices;
966 error:
967         free_chamber_list(list);
968         free_todo(todo);
969         isl_vec_free(sample);
970         isl_tab_free(tab);
971         free(selection);
972         if (!tab)
973                 isl_basic_set_free(bset);
974         isl_vertices_free(vertices);
975         return NULL;
976 }
977
978 isl_ctx *isl_vertex_get_ctx(__isl_keep isl_vertex *vertex)
979 {
980         return vertex ? isl_vertices_get_ctx(vertex->vertices) : NULL;
981 }
982
983 int isl_vertex_get_id(__isl_keep isl_vertex *vertex)
984 {
985         return vertex ? vertex->id : -1;
986 }
987
988 __isl_give isl_basic_set *isl_vertex_get_domain(__isl_keep isl_vertex *vertex)
989 {
990         struct isl_vertex *v;
991
992         if (!vertex)
993                 return NULL;
994
995         v = &vertex->vertices->v[vertex->id];
996         if (!v->dom) {
997                 v->dom = isl_basic_set_copy(v->vertex);
998                 v->dom = isl_basic_set_params(v->dom);
999         }
1000
1001         return isl_basic_set_copy(v->dom);
1002 }
1003
1004 __isl_give isl_basic_set *isl_vertex_get_expr(__isl_keep isl_vertex *vertex)
1005 {
1006         struct isl_vertex *v;
1007
1008         if (!vertex)
1009                 return NULL;
1010
1011         v = &vertex->vertices->v[vertex->id];
1012
1013         return isl_basic_set_copy(v->vertex);
1014 }
1015
1016 static __isl_give isl_vertex *isl_vertex_alloc(__isl_take isl_vertices *vertices,
1017         int id)
1018 {
1019         isl_ctx *ctx;
1020         isl_vertex *vertex;
1021
1022         if (!vertices)
1023                 return NULL;
1024
1025         ctx = isl_vertices_get_ctx(vertices);
1026         vertex = isl_alloc_type(ctx, isl_vertex);
1027         if (!vertex)
1028                 goto error;
1029
1030         vertex->vertices = vertices;
1031         vertex->id = id;
1032
1033         return vertex;
1034 error:
1035         isl_vertices_free(vertices);
1036         return NULL;
1037 }
1038
1039 void isl_vertex_free(__isl_take isl_vertex *vertex)
1040 {
1041         if (!vertex)
1042                 return;
1043         isl_vertices_free(vertex->vertices);
1044         free(vertex);
1045 }
1046
1047 __isl_give isl_basic_set *isl_basic_set_set_integral(__isl_take isl_basic_set *bset)
1048 {
1049         if (!bset)
1050                 return NULL;
1051
1052         if (!ISL_F_ISSET(bset, ISL_BASIC_MAP_RATIONAL))
1053                 return bset;
1054
1055         bset = isl_basic_set_cow(bset);
1056         if (!bset)
1057                 return NULL;
1058
1059         ISL_F_CLR(bset, ISL_BASIC_MAP_RATIONAL);
1060
1061         return isl_basic_set_finalize(bset);
1062 }
1063
1064 isl_ctx *isl_cell_get_ctx(__isl_keep isl_cell *cell)
1065 {
1066         return cell ? cell->dom->ctx : NULL;
1067 }
1068
1069 __isl_give isl_basic_set *isl_cell_get_domain(__isl_keep isl_cell *cell)
1070 {
1071         return cell ? isl_basic_set_copy(cell->dom) : NULL;
1072 }
1073
1074 static __isl_give isl_cell *isl_cell_alloc(__isl_take isl_vertices *vertices,
1075         __isl_take isl_basic_set *dom, int id)
1076 {
1077         int i;
1078         isl_cell *cell = NULL;
1079
1080         if (!vertices || !dom)
1081                 goto error;
1082
1083         cell = isl_calloc_type(dom->ctx, isl_cell);
1084         if (!cell)
1085                 goto error;
1086
1087         cell->n_vertices = vertices->c[id].n_vertices;
1088         cell->ids = isl_alloc_array(dom->ctx, int, cell->n_vertices);
1089         if (!cell->ids)
1090                 goto error;
1091         for (i = 0; i < cell->n_vertices; ++i)
1092                 cell->ids[i] = vertices->c[id].vertices[i];
1093         cell->vertices = vertices;
1094         cell->dom = dom;
1095
1096         return cell;
1097 error:
1098         isl_cell_free(cell);
1099         isl_vertices_free(vertices);
1100         isl_basic_set_free(dom);
1101         return NULL;
1102 }
1103
1104 void isl_cell_free(__isl_take isl_cell *cell)
1105 {
1106         if (!cell)
1107                 return;
1108
1109         isl_vertices_free(cell->vertices);
1110         free(cell->ids);
1111         isl_basic_set_free(cell->dom);
1112         free(cell);
1113 }
1114
1115 /* Create a tableau of the cone obtained by first homogenizing the given
1116  * polytope and then making all inequalities strict by setting the
1117  * constant term to -1.
1118  */
1119 static struct isl_tab *tab_for_shifted_cone(__isl_keep isl_basic_set *bset)
1120 {
1121         int i;
1122         isl_vec *c = NULL;
1123         struct isl_tab *tab;
1124
1125         if (!bset)
1126                 return NULL;
1127         tab = isl_tab_alloc(bset->ctx, bset->n_ineq + 1,
1128                             1 + isl_basic_set_total_dim(bset), 0);
1129         if (!tab)
1130                 return NULL;
1131         tab->rational = ISL_F_ISSET(bset, ISL_BASIC_SET_RATIONAL);
1132         if (ISL_F_ISSET(bset, ISL_BASIC_MAP_EMPTY)) {
1133                 if (isl_tab_mark_empty(tab) < 0)
1134                         goto error;
1135                 return tab;
1136         }
1137
1138         c = isl_vec_alloc(bset->ctx, 1 + 1 + isl_basic_set_total_dim(bset));
1139         if (!c)
1140                 goto error;
1141
1142         isl_int_set_si(c->el[0], 0);
1143         for (i = 0; i < bset->n_eq; ++i) {
1144                 isl_seq_cpy(c->el + 1, bset->eq[i], c->size - 1);
1145                 if (isl_tab_add_eq(tab, c->el) < 0)
1146                         goto error;
1147         }
1148
1149         isl_int_set_si(c->el[0], -1);
1150         for (i = 0; i < bset->n_ineq; ++i) {
1151                 isl_seq_cpy(c->el + 1, bset->ineq[i], c->size - 1);
1152                 if (isl_tab_add_ineq(tab, c->el) < 0)
1153                         goto error;
1154                 if (tab->empty) {
1155                         isl_vec_free(c);
1156                         return tab;
1157                 }
1158         }
1159
1160         isl_seq_clr(c->el + 1, c->size - 1);
1161         isl_int_set_si(c->el[1], 1);
1162         if (isl_tab_add_ineq(tab, c->el) < 0)
1163                 goto error;
1164
1165         isl_vec_free(c);
1166         return tab;
1167 error:
1168         isl_vec_free(c);
1169         isl_tab_free(tab);
1170         return NULL;
1171 }
1172
1173 /* Compute an interior point of "bset" by selecting an interior
1174  * point in homogeneous space and projecting the point back down.
1175  */
1176 static __isl_give isl_vec *isl_basic_set_interior_point(
1177         __isl_keep isl_basic_set *bset)
1178 {
1179         isl_vec *vec;
1180         struct isl_tab *tab;
1181
1182         tab = tab_for_shifted_cone(bset);
1183         vec = isl_tab_get_sample_value(tab);
1184         isl_tab_free(tab);
1185         if (!vec)
1186                 return NULL;
1187
1188         isl_seq_cpy(vec->el, vec->el + 1, vec->size - 1);
1189         vec->size--;
1190
1191         return vec;
1192 }
1193
1194 /* Call "fn" on all chambers of the parametric polytope with the shared
1195  * facets of neighboring chambers only appearing in one of the chambers.
1196  *
1197  * We pick an interior point from one of the chambers and then make
1198  * all constraints that do not satisfy this point strict.
1199  */
1200 int isl_vertices_foreach_disjoint_cell(__isl_keep isl_vertices *vertices,
1201         int (*fn)(__isl_take isl_cell *cell, void *user), void *user)
1202 {
1203         int i, j;
1204         isl_vec *vec;
1205         isl_int v;
1206         isl_cell *cell;
1207
1208         if (!vertices)
1209                 return -1;
1210
1211         if (vertices->n_chambers == 0)
1212                 return 0;
1213
1214         if (vertices->n_chambers == 1) {
1215                 isl_basic_set *dom = isl_basic_set_copy(vertices->c[0].dom);
1216                 dom = isl_basic_set_set_integral(dom);
1217                 cell = isl_cell_alloc(isl_vertices_copy(vertices), dom, 0);
1218                 if (!cell)
1219                         return -1;
1220                 return fn(cell, user);
1221         }
1222
1223         vec = isl_basic_set_interior_point(vertices->c[0].dom);
1224         if (!vec)
1225                 return -1;
1226
1227         isl_int_init(v);
1228
1229         for (i = 0; i < vertices->n_chambers; ++i) {
1230                 int r;
1231                 isl_basic_set *dom = isl_basic_set_copy(vertices->c[i].dom);
1232                 dom = isl_basic_set_cow(dom);
1233                 if (!dom)
1234                         goto error;
1235                 for (j = 0; i && j < dom->n_ineq; ++j) {
1236                         isl_seq_inner_product(vec->el, dom->ineq[j], vec->size,
1237                                                 &v);
1238                         if (!isl_int_is_neg(v))
1239                                 continue;
1240                         isl_int_sub_ui(dom->ineq[j][0], dom->ineq[j][0], 1);
1241                 }
1242                 dom = isl_basic_set_set_integral(dom);
1243                 cell = isl_cell_alloc(isl_vertices_copy(vertices), dom, i);
1244                 if (!cell)
1245                         goto error;
1246                 r = fn(cell, user);
1247                 if (r < 0)
1248                         goto error;
1249         }
1250
1251         isl_int_clear(v);
1252         isl_vec_free(vec);
1253
1254         return 0;
1255 error:
1256         isl_int_clear(v);
1257         isl_vec_free(vec);
1258         return -1;
1259 }
1260
1261 int isl_vertices_foreach_cell(__isl_keep isl_vertices *vertices,
1262         int (*fn)(__isl_take isl_cell *cell, void *user), void *user)
1263 {
1264         int i;
1265         isl_cell *cell;
1266
1267         if (!vertices)
1268                 return -1;
1269
1270         if (vertices->n_chambers == 0)
1271                 return 0;
1272
1273         for (i = 0; i < vertices->n_chambers; ++i) {
1274                 int r;
1275                 isl_basic_set *dom = isl_basic_set_copy(vertices->c[i].dom);
1276
1277                 cell = isl_cell_alloc(isl_vertices_copy(vertices), dom, i);
1278                 if (!cell)
1279                         return -1;
1280
1281                 r = fn(cell, user);
1282                 if (r < 0)
1283                         return -1;
1284         }
1285
1286         return 0;
1287 }
1288
1289 int isl_vertices_foreach_vertex(__isl_keep isl_vertices *vertices,
1290         int (*fn)(__isl_take isl_vertex *vertex, void *user), void *user)
1291 {
1292         int i;
1293         isl_vertex *vertex;
1294
1295         if (!vertices)
1296                 return -1;
1297
1298         if (vertices->n_vertices == 0)
1299                 return 0;
1300
1301         for (i = 0; i < vertices->n_vertices; ++i) {
1302                 int r;
1303
1304                 vertex = isl_vertex_alloc(isl_vertices_copy(vertices), i);
1305                 if (!vertex)
1306                         return -1;
1307
1308                 r = fn(vertex, user);
1309                 if (r < 0)
1310                         return -1;
1311         }
1312
1313         return 0;
1314 }
1315
1316 int isl_cell_foreach_vertex(__isl_keep isl_cell *cell,
1317         int (*fn)(__isl_take isl_vertex *vertex, void *user), void *user)
1318 {
1319         int i;
1320         isl_vertex *vertex;
1321
1322         if (!cell)
1323                 return -1;
1324
1325         if (cell->n_vertices == 0)
1326                 return 0;
1327
1328         for (i = 0; i < cell->n_vertices; ++i) {
1329                 int r;
1330
1331                 vertex = isl_vertex_alloc(isl_vertices_copy(cell->vertices),
1332                                           cell->ids[i]);
1333                 if (!vertex)
1334                         return -1;
1335
1336                 r = fn(vertex, user);
1337                 if (r < 0)
1338                         return -1;
1339         }
1340
1341         return 0;
1342 }
1343
1344 isl_ctx *isl_vertices_get_ctx(__isl_keep isl_vertices *vertices)
1345 {
1346         return vertices ? vertices->bset->ctx : NULL;
1347 }
1348
1349 int isl_vertices_get_n_vertices(__isl_keep isl_vertices *vertices)
1350 {
1351         return vertices ? vertices->n_vertices : -1;
1352 }
1353
1354 __isl_give isl_vertices *isl_morph_vertices(__isl_take isl_morph *morph,
1355         __isl_take isl_vertices *vertices)
1356 {
1357         int i;
1358         isl_morph *param_morph = NULL;
1359
1360         if (!morph || !vertices)
1361                 goto error;
1362
1363         isl_assert(vertices->bset->ctx, vertices->ref == 1, goto error);
1364
1365         param_morph = isl_morph_copy(morph);
1366         param_morph = isl_morph_dom_params(param_morph);
1367         param_morph = isl_morph_ran_params(param_morph);
1368
1369         for (i = 0; i < vertices->n_vertices; ++i) {
1370                 vertices->v[i].dom = isl_morph_basic_set(
1371                         isl_morph_copy(param_morph), vertices->v[i].dom);
1372                 vertices->v[i].vertex = isl_morph_basic_set(
1373                         isl_morph_copy(morph), vertices->v[i].vertex);
1374                 if (!vertices->v[i].vertex)
1375                         goto error;
1376         }
1377
1378         for (i = 0; i < vertices->n_chambers; ++i) {
1379                 vertices->c[i].dom = isl_morph_basic_set(
1380                         isl_morph_copy(param_morph), vertices->c[i].dom);
1381                 if (!vertices->c[i].dom)
1382                         goto error;
1383         }
1384
1385         isl_morph_free(param_morph);
1386         isl_morph_free(morph);
1387         return vertices;
1388 error:
1389         isl_morph_free(param_morph);
1390         isl_morph_free(morph);
1391         isl_vertices_free(vertices);
1392         return NULL;
1393 }
1394
1395 /* Construct a simplex isl_cell spanned by the vertices with indices in
1396  * "simplex_ids" and "other_ids" and call "fn" on this isl_cell.
1397  */
1398 static int call_on_simplex(__isl_keep isl_cell *cell,
1399         int *simplex_ids, int n_simplex, int *other_ids, int n_other,
1400         int (*fn)(__isl_take isl_cell *simplex, void *user), void *user)
1401 {
1402         int i;
1403         isl_ctx *ctx;
1404         struct isl_cell *simplex;
1405
1406         ctx = isl_cell_get_ctx(cell);
1407
1408         simplex = isl_calloc_type(ctx, struct isl_cell);
1409         if (!simplex)
1410                 return -1;
1411         simplex->vertices = isl_vertices_copy(cell->vertices);
1412         if (!simplex->vertices)
1413                 goto error;
1414         simplex->dom = isl_basic_set_copy(cell->dom);
1415         if (!simplex->dom)
1416                 goto error;
1417         simplex->n_vertices = n_simplex + n_other;
1418         simplex->ids = isl_alloc_array(ctx, int, simplex->n_vertices);
1419         if (!simplex->ids)
1420                 goto error;
1421
1422         for (i = 0; i < n_simplex; ++i)
1423                 simplex->ids[i] = simplex_ids[i];
1424         for (i = 0; i < n_other; ++i)
1425                 simplex->ids[n_simplex + i] = other_ids[i];
1426
1427         return fn(simplex, user);
1428 error:
1429         isl_cell_free(simplex);
1430         return -1;
1431 }
1432
1433 /* Check whether the parametric vertex described by "vertex"
1434  * lies on the facet corresponding to constraint "facet" of "bset".
1435  * The isl_vec "v" is a temporary vector than can be used by this function.
1436  *
1437  * We eliminate the variables from the facet constraint using the
1438  * equalities defining the vertex and check if the result is identical
1439  * to zero.
1440  *
1441  * It would probably be better to keep track of the constraints defining
1442  * a vertex during the vertex construction so that we could simply look
1443  * it up here.
1444  */
1445 static int vertex_on_facet(__isl_keep isl_basic_set *vertex,
1446         __isl_keep isl_basic_set *bset, int facet, __isl_keep isl_vec *v)
1447 {
1448         int i;
1449         isl_int m;
1450
1451         isl_seq_cpy(v->el, bset->ineq[facet], v->size);
1452
1453         isl_int_init(m);
1454         for (i = 0; i < vertex->n_eq; ++i) {
1455                 int k = isl_seq_last_non_zero(vertex->eq[i], v->size);
1456                 isl_seq_elim(v->el, vertex->eq[i], k, v->size, &m);
1457         }
1458         isl_int_clear(m);
1459
1460         return isl_seq_first_non_zero(v->el, v->size) == -1;
1461 }
1462
1463 /* Triangulate the polytope spanned by the vertices with ids
1464  * in "simplex_ids" and "other_ids" and call "fn" on each of
1465  * the resulting simplices.
1466  * If the input polytope is already a simplex, we simply call "fn".
1467  * Otherwise, we pick a point from "other_ids" and add it to "simplex_ids".
1468  * Then we consider each facet of "bset" that does not contain the point
1469  * we just picked, but does contain some of the other points in "other_ids"
1470  * and call ourselves recursively on the polytope spanned by the new
1471  * "simplex_ids" and those points in "other_ids" that lie on the facet.
1472  */
1473 static int triangulate(__isl_keep isl_cell *cell, __isl_keep isl_vec *v,
1474         int *simplex_ids, int n_simplex, int *other_ids, int n_other,
1475         int (*fn)(__isl_take isl_cell *simplex, void *user), void *user)
1476 {
1477         int i, j, k;
1478         int d, nparam;
1479         int *ids;
1480         isl_ctx *ctx;
1481         isl_basic_set *vertex;
1482         isl_basic_set *bset;
1483
1484         ctx = isl_cell_get_ctx(cell);
1485         d = isl_basic_set_dim(cell->vertices->bset, isl_dim_set);
1486         nparam = isl_basic_set_dim(cell->vertices->bset, isl_dim_param);
1487
1488         if (n_simplex + n_other == d + 1)
1489                 return call_on_simplex(cell, simplex_ids, n_simplex,
1490                                        other_ids, n_other, fn, user);
1491
1492         simplex_ids[n_simplex] = other_ids[0];
1493         vertex = cell->vertices->v[other_ids[0]].vertex;
1494         bset = cell->vertices->bset;
1495
1496         ids = isl_alloc_array(ctx, int, n_other - 1);
1497         for (i = 0; i < bset->n_ineq; ++i) {
1498                 if (isl_seq_first_non_zero(bset->ineq[i] + 1 + nparam, d) == -1)
1499                         continue;
1500                 if (vertex_on_facet(vertex, bset, i, v))
1501                         continue;
1502
1503                 for (j = 1, k = 0; j < n_other; ++j) {
1504                         isl_basic_set *ov;
1505                         ov = cell->vertices->v[other_ids[j]].vertex;
1506                         if (vertex_on_facet(ov, bset, i, v))
1507                                 ids[k++] = other_ids[j];
1508                 }
1509                 if (k == 0)
1510                         continue;
1511
1512                 if (triangulate(cell, v, simplex_ids, n_simplex + 1,
1513                                 ids, k, fn, user) < 0)
1514                         goto error;
1515         }
1516         free(ids);
1517
1518         return 0;
1519 error:
1520         free(ids);
1521         return -1;
1522 }
1523
1524 /* Triangulate the given cell and call "fn" on each of the resulting
1525  * simplices.
1526  */
1527 int isl_cell_foreach_simplex(__isl_take isl_cell *cell,
1528         int (*fn)(__isl_take isl_cell *simplex, void *user), void *user)
1529 {
1530         int d, total;
1531         int r;
1532         isl_ctx *ctx;
1533         isl_vec *v = NULL;
1534         int *simplex_ids = NULL;
1535
1536         if (!cell)
1537                 return -1;
1538
1539         d = isl_basic_set_dim(cell->vertices->bset, isl_dim_set);
1540         total = isl_basic_set_total_dim(cell->vertices->bset);
1541
1542         if (cell->n_vertices == d + 1)
1543                 return fn(cell, user);
1544
1545         ctx = isl_cell_get_ctx(cell);
1546         simplex_ids = isl_alloc_array(ctx, int, d + 1);
1547         if (!simplex_ids)
1548                 goto error;
1549
1550         v = isl_vec_alloc(ctx, 1 + total);
1551         if (!v)
1552                 goto error;
1553
1554         r = triangulate(cell, v, simplex_ids, 0,
1555                         cell->ids, cell->n_vertices, fn, user);
1556
1557         isl_vec_free(v);
1558         free(simplex_ids);
1559
1560         isl_cell_free(cell);
1561
1562         return r;
1563 error:
1564         free(simplex_ids);
1565         isl_vec_free(v);
1566         isl_cell_free(cell);
1567         return -1;
1568 }