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