1	/*
2	* Copyright 2008-2009 Katholieke Universiteit Leuven
3	* Copyright 2014 INRIA Rocquencourt
4	*
5	* Use of this software is governed by the MIT license
6	*
7	* Written by Sven Verdoolaege, K.U.Leuven, Departement
8	* Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
9	* and Inria Paris - Rocquencourt, Domaine de Voluceau - Rocquencourt,
10	* B.P. 105 - 78153 Le Chesnay, France
11	*/
12
13	#include <isl_ctx_private.h>
14	#include <isl_map_private.h>
15	#include <isl_lp_private.h>
16	#include <isl/map.h>
17	#include <isl_mat_private.h>
18	#include <isl_vec_private.h>
19	#include <isl/set.h>
20	#include <isl_seq.h>
21	#include <isl_options_private.h>
22	#include "isl_equalities.h"
23	#include "isl_tab.h"
24	#include <isl_sort.h>
25
26	#include <bset_to_bmap.c>
27	#include <bset_from_bmap.c>
28	#include <set_to_map.c>
29
30	static __isl_give isl_basic_set *uset_convex_hull_wrap_bounded(
31	__isl_take isl_set *set);
32
33	/* Remove redundant
34	* constraints. If the minimal value along the normal of a constraint
35	* is the same if the constraint is removed, then the constraint is redundant.
36	*
37	* Since some constraints may be mutually redundant, sort the constraints
38	* first such that constraints that involve existentially quantified
39	* variables are considered for removal before those that do not.
40	* The sorting is also needed for the use in map_simple_hull.
41	*
42	* Note that isl_tab_detect_implicit_equalities may also end up
43	* marking some constraints as redundant. Make sure the constraints
44	* are preserved and undo those marking such that isl_tab_detect_redundant
45	* can consider the constraints in the sorted order.
46	*
47	* Alternatively, we could have intersected the basic map with the
48	* corresponding equality and then checked if the dimension was that
49	* of a facet.
50	*/
51	__isl_give isl_basic_map *isl_basic_map_remove_redundancies(
52	__isl_take isl_basic_map *bmap)
53	{
54	struct isl_tab *tab;
55
56	if (!bmap)
57	return NULL;
58
59	bmap = isl_basic_map_gauss(bmap, NULL);
60	if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
61	return bmap;
62	if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_NO_REDUNDANT))
63	return bmap;
64	if (bmap->n_ineq <= 1)
65	return bmap;
66
67	bmap = isl_basic_map_sort_constraints(bmap);
68	tab = isl_tab_from_basic_map(bmap, track: 0);
69	if (!tab)
70	goto error;
71	tab->preserve = 1;
72	if (isl_tab_detect_implicit_equalities(tab) < 0)
73	goto error;
74	if (isl_tab_restore_redundant(tab) < 0)
75	goto error;
76	tab->preserve = 0;
77	if (isl_tab_detect_redundant(tab) < 0)
78	goto error;
79	bmap = isl_basic_map_update_from_tab(bmap, tab);
80	isl_tab_free(tab);
81	if (!bmap)
82	return NULL;
83	ISL_F_SET(bmap, ISL_BASIC_MAP_NO_IMPLICIT);
84	ISL_F_SET(bmap, ISL_BASIC_MAP_NO_REDUNDANT);
85	return bmap;
86	error:
87	isl_tab_free(tab);
88	isl_basic_map_free(bmap);
89	return NULL;
90	}
91
92	__isl_give isl_basic_set *isl_basic_set_remove_redundancies(
93	__isl_take isl_basic_set *bset)
94	{
95	return bset_from_bmap(
96	bmap: isl_basic_map_remove_redundancies(bmap: bset_to_bmap(bset)));
97	}
98
99	/* Remove redundant constraints in each of the basic maps.
100	*/
101	__isl_give isl_map *isl_map_remove_redundancies(__isl_take isl_map *map)
102	{
103	return isl_map_inline_foreach_basic_map(map,
104	fn: &isl_basic_map_remove_redundancies);
105	}
106
107	__isl_give isl_set *isl_set_remove_redundancies(__isl_take isl_set *set)
108	{
109	return isl_map_remove_redundancies(map: set);
110	}
111
112	/* Check if the set set is bound in the direction of the affine
113	* constraint c and if so, set the constant term such that the
114	* resulting constraint is a bounding constraint for the set.
115	*/
116	static isl_bool uset_is_bound(__isl_keep isl_set *set, isl_int *c, unsigned len)
117	{
118	int first;
119	int j;
120	isl_int opt;
121	isl_int opt_denom;
122
123	isl_int_init(opt);
124	isl_int_init(opt_denom);
125	first = 1;
126	for (j = 0; j < set->n; ++j) {
127	enum isl_lp_result res;
128
129	if (ISL_F_ISSET(set->p[j], ISL_BASIC_SET_EMPTY))
130	continue;
131
132	res = isl_basic_set_solve_lp(bset: set->p[j],
133	max: 0, f: c, denom: set->ctx->one, opt: &opt, opt_denom: &opt_denom, NULL);
134	if (res == isl_lp_unbounded)
135	break;
136	if (res == isl_lp_error)
137	goto error;
138	if (res == isl_lp_empty) {
139	set->p[j] = isl_basic_set_set_to_empty(bset: set->p[j]);
140	if (!set->p[j])
141	goto error;
142	continue;
143	}
144	if (first || isl_int_is_neg(opt)) {
145	if (!isl_int_is_one(opt_denom))
146	isl_seq_scale(dst: c, src: c, f: opt_denom, len);
147	isl_int_sub(c[0], c[0], opt);
148	}
149	first = 0;
150	}
151	isl_int_clear(opt);
152	isl_int_clear(opt_denom);
153	return isl_bool_ok(b: j >= set->n);
154	error:
155	isl_int_clear(opt);
156	isl_int_clear(opt_denom);
157	return isl_bool_error;
158	}
159
160	static __isl_give isl_set *isl_set_add_basic_set_equality(
161	__isl_take isl_set *set, isl_int *c)
162	{
163	int i;
164
165	set = isl_set_cow(set);
166	if (!set)
167	return NULL;
168	for (i = 0; i < set->n; ++i) {
169	set->p[i] = isl_basic_set_add_eq(bset: set->p[i], eq: c);
170	if (!set->p[i])
171	goto error;
172	}
173	return set;
174	error:
175	isl_set_free(set);
176	return NULL;
177	}
178
179	/* Given a union of basic sets, construct the constraints for wrapping
180	* a facet around one of its ridges.
181	* In particular, if each of n the d-dimensional basic sets i in "set"
182	* contains the origin, satisfies the constraints x_1 >= 0 and x_2 >= 0
183	* and is defined by the constraints
184	* [ 1 ]
185	* A_i [ x ] >= 0
186	*
187	* then the resulting set is of dimension n*(1+d) and has as constraints
188	*
189	* [ a_i ]
190	* A_i [ x_i ] >= 0
191	*
192	* a_i >= 0
193	*
194	* \sum_i x_{i,1} = 1
195	*/
196	static __isl_give isl_basic_set *wrap_constraints(__isl_keep isl_set *set)
197	{
198	struct isl_basic_set *lp;
199	unsigned n_eq;
200	unsigned n_ineq;
201	int i, j, k;
202	isl_size dim, lp_dim;
203
204	dim = isl_set_dim(set, type: isl_dim_set);
205	if (dim < 0)
206	return NULL;
207
208	dim += 1;
209	n_eq = 1;
210	n_ineq = set->n;
211	for (i = 0; i < set->n; ++i) {
212	n_eq += set->p[i]->n_eq;
213	n_ineq += set->p[i]->n_ineq;
214	}
215	lp = isl_basic_set_alloc(ctx: set->ctx, nparam: 0, dim: dim * set->n, extra: 0, n_eq, n_ineq);
216	lp = isl_basic_set_set_rational(bset: lp);
217	if (!lp)
218	return NULL;
219	lp_dim = isl_basic_set_dim(bset: lp, type: isl_dim_set);
220	if (lp_dim < 0)
221	return isl_basic_set_free(bset: lp);
222	k = isl_basic_set_alloc_equality(bset: lp);
223	isl_int_set_si(lp->eq[k][0], -1);
224	for (i = 0; i < set->n; ++i) {
225	isl_int_set_si(lp->eq[k][1+dim*i], 0);
226	isl_int_set_si(lp->eq[k][1+dim*i+1], 1);
227	isl_seq_clr(p: lp->eq[k]+1+dim*i+2, len: dim-2);
228	}
229	for (i = 0; i < set->n; ++i) {
230	k = isl_basic_set_alloc_inequality(bset: lp);
231	isl_seq_clr(p: lp->ineq[k], len: 1+lp_dim);
232	isl_int_set_si(lp->ineq[k][1+dim*i], 1);
233
234	for (j = 0; j < set->p[i]->n_eq; ++j) {
235	k = isl_basic_set_alloc_equality(bset: lp);
236	isl_seq_clr(p: lp->eq[k], len: 1+dim*i);
237	isl_seq_cpy(dst: lp->eq[k]+1+dim*i, src: set->p[i]->eq[j], len: dim);
238	isl_seq_clr(p: lp->eq[k]+1+dim*(i+1), len: dim*(set->n-i-1));
239	}
240
241	for (j = 0; j < set->p[i]->n_ineq; ++j) {
242	k = isl_basic_set_alloc_inequality(bset: lp);
243	isl_seq_clr(p: lp->ineq[k], len: 1+dim*i);
244	isl_seq_cpy(dst: lp->ineq[k]+1+dim*i, src: set->p[i]->ineq[j], len: dim);
245	isl_seq_clr(p: lp->ineq[k]+1+dim*(i+1), len: dim*(set->n-i-1));
246	}
247	}
248	return lp;
249	}
250
251	/* Given a facet "facet" of the convex hull of "set" and a facet "ridge"
252	* of that facet, compute the other facet of the convex hull that contains
253	* the ridge.
254	*
255	* We first transform the set such that the facet constraint becomes
256	*
257	* x_1 >= 0
258	*
259	* I.e., the facet lies in
260	*
261	* x_1 = 0
262	*
263	* and on that facet, the constraint that defines the ridge is
264	*
265	* x_2 >= 0
266	*
267	* (This transformation is not strictly needed, all that is needed is
268	* that the ridge contains the origin.)
269	*
270	* Since the ridge contains the origin, the cone of the convex hull
271	* will be of the form
272	*
273	* x_1 >= 0
274	* x_2 >= a x_1
275	*
276	* with this second constraint defining the new facet.
277	* The constant a is obtained by settting x_1 in the cone of the
278	* convex hull to 1 and minimizing x_2.
279	* Now, each element in the cone of the convex hull is the sum
280	* of elements in the cones of the basic sets.
281	* If a_i is the dilation factor of basic set i, then the problem
282	* we need to solve is
283	*
284	* min \sum_i x_{i,2}
285	* st
286	* \sum_i x_{i,1} = 1
287	* a_i >= 0
288	* [ a_i ]
289	* A [ x_i ] >= 0
290	*
291	* with
292	* [ 1 ]
293	* A_i [ x_i ] >= 0
294	*
295	* the constraints of each (transformed) basic set.
296	* If a = n/d, then the constraint defining the new facet (in the transformed
297	* space) is
298	*
299	* -n x_1 + d x_2 >= 0
300	*
301	* In the original space, we need to take the same combination of the
302	* corresponding constraints "facet" and "ridge".
303	*
304	* If a = -infty = "-1/0", then we just return the original facet constraint.
305	* This means that the facet is unbounded, but has a bounded intersection
306	* with the union of sets.
307	*/
308	isl_int *isl_set_wrap_facet(__isl_keep isl_set *set,
309	isl_int *facet, isl_int *ridge)
310	{
311	int i;
312	isl_ctx *ctx;
313	struct isl_mat *T = NULL;
314	struct isl_basic_set *lp = NULL;
315	struct isl_vec *obj;
316	enum isl_lp_result res;
317	isl_int num, den;
318	isl_size dim;
319
320	dim = isl_set_dim(set, type: isl_dim_set);
321	if (dim < 0)
322	return NULL;
323	ctx = set->ctx;
324	set = isl_set_copy(set);
325	set = isl_set_set_rational(set);
326
327	dim += 1;
328	T = isl_mat_alloc(ctx, n_row: 3, n_col: dim);
329	if (!T)
330	goto error;
331	isl_int_set_si(T->row[0][0], 1);
332	isl_seq_clr(p: T->row[0]+1, len: dim - 1);
333	isl_seq_cpy(dst: T->row[1], src: facet, len: dim);
334	isl_seq_cpy(dst: T->row[2], src: ridge, len: dim);
335	T = isl_mat_right_inverse(mat: T);
336	set = isl_set_preimage(set, mat: T);
337	T = NULL;
338	if (!set)
339	goto error;
340	lp = wrap_constraints(set);
341	obj = isl_vec_alloc(ctx, size: 1 + dim*set->n);
342	if (!obj)
343	goto error;
344	isl_int_set_si(obj->block.data[0], 0);
345	for (i = 0; i < set->n; ++i) {
346	isl_seq_clr(p: obj->block.data + 1 + dim*i, len: 2);
347	isl_int_set_si(obj->block.data[1 + dim*i+2], 1);
348	isl_seq_clr(p: obj->block.data + 1 + dim*i+3, len: dim-3);
349	}
350	isl_int_init(num);
351	isl_int_init(den);
352	res = isl_basic_set_solve_lp(bset: lp, max: 0,
353	f: obj->block.data, denom: ctx->one, opt: &num, opt_denom: &den, NULL);
354	if (res == isl_lp_ok) {
355	isl_int_neg(num, num);
356	isl_seq_combine(dst: facet, m1: num, src1: facet, m2: den, src2: ridge, len: dim);
357	isl_seq_normalize(ctx, p: facet, len: dim);
358	}
359	isl_int_clear(num);
360	isl_int_clear(den);
361	isl_vec_free(vec: obj);
362	isl_basic_set_free(bset: lp);
363	isl_set_free(set);
364	if (res == isl_lp_error)
365	return NULL;
366	isl_assert(ctx, res == isl_lp_ok || res == isl_lp_unbounded,
367	return NULL);
368	return facet;
369	error:
370	isl_basic_set_free(bset: lp);
371	isl_mat_free(mat: T);
372	isl_set_free(set);
373	return NULL;
374	}
375
376	/* Compute the constraint of a facet of "set".
377	*
378	* We first compute the intersection with a bounding constraint
379	* that is orthogonal to one of the coordinate axes.
380	* If the affine hull of this intersection has only one equality,
381	* we have found a facet.
382	* Otherwise, we wrap the current bounding constraint around
383	* one of the equalities of the face (one that is not equal to
384	* the current bounding constraint).
385	* This process continues until we have found a facet.
386	* The dimension of the intersection increases by at least
387	* one on each iteration, so termination is guaranteed.
388	*/
389	static __isl_give isl_mat *initial_facet_constraint(__isl_keep isl_set *set)
390	{
391	struct isl_set *slice = NULL;
392	struct isl_basic_set *face = NULL;
393	int i;
394	isl_size dim = isl_set_dim(set, type: isl_dim_set);
395	isl_bool is_bound;
396	isl_mat *bounds = NULL;
397
398	if (dim < 0)
399	return NULL;
400	isl_assert(set->ctx, set->n > 0, goto error);
401	bounds = isl_mat_alloc(ctx: set->ctx, n_row: 1, n_col: 1 + dim);
402	if (!bounds)
403	return NULL;
404
405	isl_seq_clr(p: bounds->row[0], len: dim);
406	isl_int_set_si(bounds->row[0][1 + dim - 1], 1);
407	is_bound = uset_is_bound(set, c: bounds->row[0], len: 1 + dim);
408	if (is_bound < 0)
409	goto error;
410	isl_assert(set->ctx, is_bound, goto error);
411	isl_seq_normalize(ctx: set->ctx, p: bounds->row[0], len: 1 + dim);
412	bounds->n_row = 1;
413
414	for (;;) {
415	slice = isl_set_copy(set);
416	slice = isl_set_add_basic_set_equality(set: slice, c: bounds->row[0]);
417	face = isl_set_affine_hull(set: slice);
418	if (!face)
419	goto error;
420	if (face->n_eq == 1) {
421	isl_basic_set_free(bset: face);
422	break;
423	}
424	for (i = 0; i < face->n_eq; ++i)
425	if (!isl_seq_eq(p1: bounds->row[0], p2: face->eq[i], len: 1 + dim) &&
426	!isl_seq_is_neg(p1: bounds->row[0],
427	p2: face->eq[i], len: 1 + dim))
428	break;
429	isl_assert(set->ctx, i < face->n_eq, goto error);
430	if (!isl_set_wrap_facet(set, facet: bounds->row[0], ridge: face->eq[i]))
431	goto error;
432	isl_seq_normalize(ctx: set->ctx, p: bounds->row[0], len: bounds->n_col);
433	isl_basic_set_free(bset: face);
434	}
435
436	return bounds;
437	error:
438	isl_basic_set_free(bset: face);
439	isl_mat_free(mat: bounds);
440	return NULL;
441	}
442
443	/* Given the bounding constraint "c" of a facet of the convex hull of "set",
444	* compute a hyperplane description of the facet, i.e., compute the facets
445	* of the facet.
446	*
447	* We compute an affine transformation that transforms the constraint
448	*
449	* [ 1 ]
450	* c [ x ] = 0
451	*
452	* to the constraint
453	*
454	* z_1 = 0
455	*
456	* by computing the right inverse U of a matrix that starts with the rows
457	*
458	* [ 1 0 ]
459	* [ c ]
460	*
461	* Then
462	* [ 1 ] [ 1 ]
463	* [ x ] = U [ z ]
464	* and
465	* [ 1 ] [ 1 ]
466	* [ z ] = Q [ x ]
467	*
468	* with Q = U^{-1}
469	* Since z_1 is zero, we can drop this variable as well as the corresponding
470	* column of U to obtain
471	*
472	* [ 1 ] [ 1 ]
473	* [ x ] = U' [ z' ]
474	* and
475	* [ 1 ] [ 1 ]
476	* [ z' ] = Q' [ x ]
477	*
478	* with Q' equal to Q, but without the corresponding row.
479	* After computing the facets of the facet in the z' space,
480	* we convert them back to the x space through Q.
481	*/
482	static __isl_give isl_basic_set *compute_facet(__isl_keep isl_set *set,
483	isl_int *c)
484	{
485	struct isl_mat *m, *U, *Q;
486	struct isl_basic_set *facet = NULL;
487	struct isl_ctx *ctx;
488	isl_size dim;
489
490	dim = isl_set_dim(set, type: isl_dim_set);
491	if (dim < 0)
492	return NULL;
493	ctx = set->ctx;
494	set = isl_set_copy(set);
495	m = isl_mat_alloc(ctx: set->ctx, n_row: 2, n_col: 1 + dim);
496	if (!m)
497	goto error;
498	isl_int_set_si(m->row[0][0], 1);
499	isl_seq_clr(p: m->row[0]+1, len: dim);
500	isl_seq_cpy(dst: m->row[1], src: c, len: 1+dim);
501	U = isl_mat_right_inverse(mat: m);
502	Q = isl_mat_right_inverse(mat: isl_mat_copy(mat: U));
503	U = isl_mat_drop_cols(mat: U, col: 1, n: 1);
504	Q = isl_mat_drop_rows(mat: Q, row: 1, n: 1);
505	set = isl_set_preimage(set, mat: U);
506	facet = uset_convex_hull_wrap_bounded(set);
507	facet = isl_basic_set_preimage(bset: facet, mat: Q);
508	if (facet && facet->n_eq != 0)
509	isl_die(ctx, isl_error_internal, "unexpected equality",
510	return isl_basic_set_free(facet));
511	return facet;
512	error:
513	isl_basic_set_free(bset: facet);
514	isl_set_free(set);
515	return NULL;
516	}
517
518	/* Given an initial facet constraint, compute the remaining facets.
519	* We do this by running through all facets found so far and computing
520	* the adjacent facets through wrapping, adding those facets that we
521	* hadn't already found before.
522	*
523	* For each facet we have found so far, we first compute its facets
524	* in the resulting convex hull. That is, we compute the ridges
525	* of the resulting convex hull contained in the facet.
526	* We also compute the corresponding facet in the current approximation
527	* of the convex hull. There is no need to wrap around the ridges
528	* in this facet since that would result in a facet that is already
529	* present in the current approximation.
530	*
531	* This function can still be significantly optimized by checking which of
532	* the facets of the basic sets are also facets of the convex hull and
533	* using all the facets so far to help in constructing the facets of the
534	* facets
535	* and/or
536	* using the technique in section "3.1 Ridge Generation" of
537	* "Extended Convex Hull" by Fukuda et al.
538	*/
539	static __isl_give isl_basic_set *extend(__isl_take isl_basic_set *hull,
540	__isl_keep isl_set *set)
541	{
542	int i, j, f;
543	int k;
544	struct isl_basic_set *facet = NULL;
545	struct isl_basic_set *hull_facet = NULL;
546	isl_size dim;
547
548	dim = isl_set_dim(set, type: isl_dim_set);
549	if (dim < 0 || !hull)
550	return isl_basic_set_free(bset: hull);
551
552	isl_assert(set->ctx, set->n > 0, goto error);
553
554	for (i = 0; i < hull->n_ineq; ++i) {
555	facet = compute_facet(set, c: hull->ineq[i]);
556	facet = isl_basic_set_add_eq(bset: facet, eq: hull->ineq[i]);
557	facet = isl_basic_set_gauss(bset: facet, NULL);
558	facet = isl_basic_set_normalize_constraints(bset: facet);
559	hull_facet = isl_basic_set_copy(bset: hull);
560	hull_facet = isl_basic_set_add_eq(bset: hull_facet, eq: hull->ineq[i]);
561	hull_facet = isl_basic_set_gauss(bset: hull_facet, NULL);
562	hull_facet = isl_basic_set_normalize_constraints(bset: hull_facet);
563	if (!facet || !hull_facet)
564	goto error;
565	hull = isl_basic_set_cow(bset: hull);
566	hull = isl_basic_set_extend(base: hull, extra: 0, n_eq: 0, n_ineq: facet->n_ineq);
567	if (!hull)
568	goto error;
569	for (j = 0; j < facet->n_ineq; ++j) {
570	for (f = 0; f < hull_facet->n_ineq; ++f)
571	if (isl_seq_eq(p1: facet->ineq[j],
572	p2: hull_facet->ineq[f], len: 1 + dim))
573	break;
574	if (f < hull_facet->n_ineq)
575	continue;
576	k = isl_basic_set_alloc_inequality(bset: hull);
577	if (k < 0)
578	goto error;
579	isl_seq_cpy(dst: hull->ineq[k], src: hull->ineq[i], len: 1+dim);
580	if (!isl_set_wrap_facet(set, facet: hull->ineq[k], ridge: facet->ineq[j]))
581	goto error;
582	}
583	isl_basic_set_free(bset: hull_facet);
584	isl_basic_set_free(bset: facet);
585	}
586	hull = isl_basic_set_simplify(bset: hull);
587	hull = isl_basic_set_finalize(bset: hull);
588	return hull;
589	error:
590	isl_basic_set_free(bset: hull_facet);
591	isl_basic_set_free(bset: facet);
592	isl_basic_set_free(bset: hull);
593	return NULL;
594	}
595
596	/* Special case for computing the convex hull of a one dimensional set.
597	* We simply collect the lower and upper bounds of each basic set
598	* and the biggest of those.
599	*/
600	static __isl_give isl_basic_set *convex_hull_1d(__isl_take isl_set *set)
601	{
602	struct isl_mat *c = NULL;
603	isl_int *lower = NULL;
604	isl_int *upper = NULL;
605	int i, j, k;
606	isl_int a, b;
607	struct isl_basic_set *hull;
608
609	for (i = 0; i < set->n; ++i) {
610	set->p[i] = isl_basic_set_simplify(bset: set->p[i]);
611	if (!set->p[i])
612	goto error;
613	}
614	set = isl_set_remove_empty_parts(set);
615	if (!set)
616	goto error;
617	isl_assert(set->ctx, set->n > 0, goto error);
618	c = isl_mat_alloc(ctx: set->ctx, n_row: 2, n_col: 2);
619	if (!c)
620	goto error;
621
622	if (set->p[0]->n_eq > 0) {
623	isl_assert(set->ctx, set->p[0]->n_eq == 1, goto error);
624	lower = c->row[0];
625	upper = c->row[1];
626	if (isl_int_is_pos(set->p[0]->eq[0][1])) {
627	isl_seq_cpy(dst: lower, src: set->p[0]->eq[0], len: 2);
628	isl_seq_neg(dst: upper, src: set->p[0]->eq[0], len: 2);
629	} else {
630	isl_seq_neg(dst: lower, src: set->p[0]->eq[0], len: 2);
631	isl_seq_cpy(dst: upper, src: set->p[0]->eq[0], len: 2);
632	}
633	} else {
634	for (j = 0; j < set->p[0]->n_ineq; ++j) {
635	if (isl_int_is_pos(set->p[0]->ineq[j][1])) {
636	lower = c->row[0];
637	isl_seq_cpy(dst: lower, src: set->p[0]->ineq[j], len: 2);
638	} else {
639	upper = c->row[1];
640	isl_seq_cpy(dst: upper, src: set->p[0]->ineq[j], len: 2);
641	}
642	}
643	}
644
645	isl_int_init(a);
646	isl_int_init(b);
647	for (i = 0; i < set->n; ++i) {
648	struct isl_basic_set *bset = set->p[i];
649	int has_lower = 0;
650	int has_upper = 0;
651
652	for (j = 0; j < bset->n_eq; ++j) {
653	has_lower = 1;
654	has_upper = 1;
655	if (lower) {
656	isl_int_mul(a, lower[0], bset->eq[j][1]);
657	isl_int_mul(b, lower[1], bset->eq[j][0]);
658	if (isl_int_lt(a, b) && isl_int_is_pos(bset->eq[j][1]))
659	isl_seq_cpy(dst: lower, src: bset->eq[j], len: 2);
660	if (isl_int_gt(a, b) && isl_int_is_neg(bset->eq[j][1]))
661	isl_seq_neg(dst: lower, src: bset->eq[j], len: 2);
662	}
663	if (upper) {
664	isl_int_mul(a, upper[0], bset->eq[j][1]);
665	isl_int_mul(b, upper[1], bset->eq[j][0]);
666	if (isl_int_lt(a, b) && isl_int_is_pos(bset->eq[j][1]))
667	isl_seq_neg(dst: upper, src: bset->eq[j], len: 2);
668	if (isl_int_gt(a, b) && isl_int_is_neg(bset->eq[j][1]))
669	isl_seq_cpy(dst: upper, src: bset->eq[j], len: 2);
670	}
671	}
672	for (j = 0; j < bset->n_ineq; ++j) {
673	if (isl_int_is_pos(bset->ineq[j][1]))
674	has_lower = 1;
675	if (isl_int_is_neg(bset->ineq[j][1]))
676	has_upper = 1;
677	if (lower && isl_int_is_pos(bset->ineq[j][1])) {
678	isl_int_mul(a, lower[0], bset->ineq[j][1]);
679	isl_int_mul(b, lower[1], bset->ineq[j][0]);
680	if (isl_int_lt(a, b))
681	isl_seq_cpy(dst: lower, src: bset->ineq[j], len: 2);
682	}
683	if (upper && isl_int_is_neg(bset->ineq[j][1])) {
684	isl_int_mul(a, upper[0], bset->ineq[j][1]);
685	isl_int_mul(b, upper[1], bset->ineq[j][0]);
686	if (isl_int_gt(a, b))
687	isl_seq_cpy(dst: upper, src: bset->ineq[j], len: 2);
688	}
689	}
690	if (!has_lower)
691	lower = NULL;
692	if (!has_upper)
693	upper = NULL;
694	}
695	isl_int_clear(a);
696	isl_int_clear(b);
697
698	hull = isl_basic_set_alloc(ctx: set->ctx, nparam: 0, dim: 1, extra: 0, n_eq: 0, n_ineq: 2);
699	hull = isl_basic_set_set_rational(bset: hull);
700	if (!hull)
701	goto error;
702	if (lower) {
703	k = isl_basic_set_alloc_inequality(bset: hull);
704	isl_seq_cpy(dst: hull->ineq[k], src: lower, len: 2);
705	}
706	if (upper) {
707	k = isl_basic_set_alloc_inequality(bset: hull);
708	isl_seq_cpy(dst: hull->ineq[k], src: upper, len: 2);
709	}
710	hull = isl_basic_set_finalize(bset: hull);
711	isl_set_free(set);
712	isl_mat_free(mat: c);
713	return hull;
714	error:
715	isl_set_free(set);
716	isl_mat_free(mat: c);
717	return NULL;
718	}
719
720	static __isl_give isl_basic_set *convex_hull_0d(__isl_take isl_set *set)
721	{
722	struct isl_basic_set *convex_hull;
723
724	if (!set)
725	return NULL;
726
727	if (isl_set_is_empty(set))
728	convex_hull = isl_basic_set_empty(space: isl_space_copy(space: set->dim));
729	else
730	convex_hull = isl_basic_set_universe(space: isl_space_copy(space: set->dim));
731	isl_set_free(set);
732	return convex_hull;
733	}
734
735	/* Compute the convex hull of a pair of basic sets without any parameters or
736	* integer divisions using Fourier-Motzkin elimination.
737	* The convex hull is the set of all points that can be written as
738	* the sum of points from both basic sets (in homogeneous coordinates).
739	* We set up the constraints in a space with dimensions for each of
740	* the three sets and then project out the dimensions corresponding
741	* to the two original basic sets, retaining only those corresponding
742	* to the convex hull.
743	*/
744	static __isl_give isl_basic_set *convex_hull_pair_elim(
745	__isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)
746	{
747	int i, j, k;
748	struct isl_basic_set *bset[2];
749	struct isl_basic_set *hull = NULL;
750	isl_size dim;
751
752	dim = isl_basic_set_dim(bset: bset1, type: isl_dim_set);
753	if (dim < 0 || !bset2)
754	goto error;
755
756	hull = isl_basic_set_alloc(ctx: bset1->ctx, nparam: 0, dim: 2 + 3 * dim, extra: 0,
757	n_eq: 1 + dim + bset1->n_eq + bset2->n_eq,
758	n_ineq: 2 + bset1->n_ineq + bset2->n_ineq);
759	bset[0] = bset1;
760	bset[1] = bset2;
761	for (i = 0; i < 2; ++i) {
762	for (j = 0; j < bset[i]->n_eq; ++j) {
763	k = isl_basic_set_alloc_equality(bset: hull);
764	if (k < 0)
765	goto error;
766	isl_seq_clr(p: hull->eq[k], len: (i+1) * (1+dim));
767	isl_seq_clr(p: hull->eq[k]+(i+2)*(1+dim), len: (1-i)*(1+dim));
768	isl_seq_cpy(dst: hull->eq[k]+(i+1)*(1+dim), src: bset[i]->eq[j],
769	len: 1+dim);
770	}
771	for (j = 0; j < bset[i]->n_ineq; ++j) {
772	k = isl_basic_set_alloc_inequality(bset: hull);
773	if (k < 0)
774	goto error;
775	isl_seq_clr(p: hull->ineq[k], len: (i+1) * (1+dim));
776	isl_seq_clr(p: hull->ineq[k]+(i+2)*(1+dim), len: (1-i)*(1+dim));
777	isl_seq_cpy(dst: hull->ineq[k]+(i+1)*(1+dim),
778	src: bset[i]->ineq[j], len: 1+dim);
779	}
780	k = isl_basic_set_alloc_inequality(bset: hull);
781	if (k < 0)
782	goto error;
783	isl_seq_clr(p: hull->ineq[k], len: 1+2+3*dim);
784	isl_int_set_si(hull->ineq[k][(i+1)*(1+dim)], 1);
785	}
786	for (j = 0; j < 1+dim; ++j) {
787	k = isl_basic_set_alloc_equality(bset: hull);
788	if (k < 0)
789	goto error;
790	isl_seq_clr(p: hull->eq[k], len: 1+2+3*dim);
791	isl_int_set_si(hull->eq[k][j], -1);
792	isl_int_set_si(hull->eq[k][1+dim+j], 1);
793	isl_int_set_si(hull->eq[k][2*(1+dim)+j], 1);
794	}
795	hull = isl_basic_set_set_rational(bset: hull);
796	hull = isl_basic_set_remove_dims(bset: hull, type: isl_dim_set, first: dim, n: 2*(1+dim));
797	hull = isl_basic_set_remove_redundancies(bset: hull);
798	isl_basic_set_free(bset: bset1);
799	isl_basic_set_free(bset: bset2);
800	return hull;
801	error:
802	isl_basic_set_free(bset: bset1);
803	isl_basic_set_free(bset: bset2);
804	isl_basic_set_free(bset: hull);
805	return NULL;
806	}
807
808	/* Is the set bounded for each value of the parameters?
809	*/
810	isl_bool isl_basic_set_is_bounded(__isl_keep isl_basic_set *bset)
811	{
812	struct isl_tab *tab;
813	isl_bool bounded;
814
815	if (!bset)
816	return isl_bool_error;
817	if (isl_basic_set_plain_is_empty(bset))
818	return isl_bool_true;
819
820	tab = isl_tab_from_recession_cone(bset, parametric: 1);
821	bounded = isl_tab_cone_is_bounded(tab);
822	isl_tab_free(tab);
823	return bounded;
824	}
825
826	/* Is the image bounded for each value of the parameters and
827	* the domain variables?
828	*/
829	isl_bool isl_basic_map_image_is_bounded(__isl_keep isl_basic_map *bmap)
830	{
831	isl_size nparam = isl_basic_map_dim(bmap, type: isl_dim_param);
832	isl_size n_in = isl_basic_map_dim(bmap, type: isl_dim_in);
833	isl_bool bounded;
834
835	if (nparam < 0 || n_in < 0)
836	return isl_bool_error;
837
838	bmap = isl_basic_map_copy(bmap);
839	bmap = isl_basic_map_cow(bmap);
840	bmap = isl_basic_map_move_dims(bmap, dst_type: isl_dim_param, dst_pos: nparam,
841	src_type: isl_dim_in, src_pos: 0, n: n_in);
842	bounded = isl_basic_set_is_bounded(bset: bset_from_bmap(bmap));
843	isl_basic_map_free(bmap);
844
845	return bounded;
846	}
847
848	/* Is the set bounded for each value of the parameters?
849	*/
850	isl_bool isl_set_is_bounded(__isl_keep isl_set *set)
851	{
852	int i;
853
854	if (!set)
855	return isl_bool_error;
856
857	for (i = 0; i < set->n; ++i) {
858	isl_bool bounded = isl_basic_set_is_bounded(bset: set->p[i]);
859	if (!bounded || bounded < 0)
860	return bounded;
861	}
862	return isl_bool_true;
863	}
864
865	/* Compute the lineality space of the convex hull of bset1 and bset2.
866	*
867	* We first compute the intersection of the recession cone of bset1
868	* with the negative of the recession cone of bset2 and then compute
869	* the linear hull of the resulting cone.
870	*/
871	static __isl_give isl_basic_set *induced_lineality_space(
872	__isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)
873	{
874	int i, k;
875	struct isl_basic_set *lin = NULL;
876	isl_size dim;
877
878	dim = isl_basic_set_dim(bset: bset1, type: isl_dim_all);
879	if (dim < 0 || !bset2)
880	goto error;
881
882	lin = isl_basic_set_alloc_space(space: isl_basic_set_get_space(bset: bset1), extra: 0,
883	n_eq: bset1->n_eq + bset2->n_eq,
884	n_ineq: bset1->n_ineq + bset2->n_ineq);
885	lin = isl_basic_set_set_rational(bset: lin);
886	if (!lin)
887	goto error;
888	for (i = 0; i < bset1->n_eq; ++i) {
889	k = isl_basic_set_alloc_equality(bset: lin);
890	if (k < 0)
891	goto error;
892	isl_int_set_si(lin->eq[k][0], 0);
893	isl_seq_cpy(dst: lin->eq[k] + 1, src: bset1->eq[i] + 1, len: dim);
894	}
895	for (i = 0; i < bset1->n_ineq; ++i) {
896	k = isl_basic_set_alloc_inequality(bset: lin);
897	if (k < 0)
898	goto error;
899	isl_int_set_si(lin->ineq[k][0], 0);
900	isl_seq_cpy(dst: lin->ineq[k] + 1, src: bset1->ineq[i] + 1, len: dim);
901	}
902	for (i = 0; i < bset2->n_eq; ++i) {
903	k = isl_basic_set_alloc_equality(bset: lin);
904	if (k < 0)
905	goto error;
906	isl_int_set_si(lin->eq[k][0], 0);
907	isl_seq_neg(dst: lin->eq[k] + 1, src: bset2->eq[i] + 1, len: dim);
908	}
909	for (i = 0; i < bset2->n_ineq; ++i) {
910	k = isl_basic_set_alloc_inequality(bset: lin);
911	if (k < 0)
912	goto error;
913	isl_int_set_si(lin->ineq[k][0], 0);
914	isl_seq_neg(dst: lin->ineq[k] + 1, src: bset2->ineq[i] + 1, len: dim);
915	}
916
917	isl_basic_set_free(bset: bset1);
918	isl_basic_set_free(bset: bset2);
919	return isl_basic_set_affine_hull(bset: lin);
920	error:
921	isl_basic_set_free(bset: lin);
922	isl_basic_set_free(bset: bset1);
923	isl_basic_set_free(bset: bset2);
924	return NULL;
925	}
926
927	static __isl_give isl_basic_set *uset_convex_hull(__isl_take isl_set *set);
928
929	/* Given a set and a linear space "lin" of dimension n > 0,
930	* project the linear space from the set, compute the convex hull
931	* and then map the set back to the original space.
932	*
933	* Let
934	*
935	* M x = 0
936	*
937	* describe the linear space. We first compute the Hermite normal
938	* form H = M U of M = H Q, to obtain
939	*
940	* H Q x = 0
941	*
942	* The last n rows of H will be zero, so the last n variables of x' = Q x
943	* are the one we want to project out. We do this by transforming each
944	* basic set A x >= b to A U x' >= b and then removing the last n dimensions.
945	* After computing the convex hull in x'_1, i.e., A' x'_1 >= b',
946	* we transform the hull back to the original space as A' Q_1 x >= b',
947	* with Q_1 all but the last n rows of Q.
948	*/
949	static __isl_give isl_basic_set *modulo_lineality(__isl_take isl_set *set,
950	__isl_take isl_basic_set *lin)
951	{
952	isl_size total = isl_basic_set_dim(bset: lin, type: isl_dim_all);
953	unsigned lin_dim;
954	struct isl_basic_set *hull;
955	struct isl_mat *M, *U, *Q;
956
957	if (!set || total < 0)
958	goto error;
959	lin_dim = total - lin->n_eq;
960	M = isl_mat_sub_alloc6(ctx: set->ctx, row: lin->eq, first_row: 0, n_row: lin->n_eq, first_col: 1, n_col: total);
961	M = isl_mat_left_hermite(M, neg: 0, U: &U, Q: &Q);
962	if (!M)
963	goto error;
964	isl_mat_free(mat: M);
965	isl_basic_set_free(bset: lin);
966
967	Q = isl_mat_drop_rows(mat: Q, row: Q->n_row - lin_dim, n: lin_dim);
968
969	U = isl_mat_lin_to_aff(mat: U);
970	Q = isl_mat_lin_to_aff(mat: Q);
971
972	set = isl_set_preimage(set, mat: U);
973	set = isl_set_remove_dims(bset: set, type: isl_dim_set, first: total - lin_dim, n: lin_dim);
974	hull = uset_convex_hull(set);
975	hull = isl_basic_set_preimage(bset: hull, mat: Q);
976
977	return hull;
978	error:
979	isl_basic_set_free(bset: lin);
980	isl_set_free(set);
981	return NULL;
982	}
983
984	/* Given two polyhedra with as constraints h_{ij} x >= 0 in homegeneous space,
985	* set up an LP for solving
986	*
987	* \sum_j \alpha_{1j} h_{1j} = \sum_j \alpha_{2j} h_{2j}
988	*
989	* \alpha{i0} corresponds to the (implicit) positivity constraint 1 >= 0
990	* The next \alpha{ij} correspond to the equalities and come in pairs.
991	* The final \alpha{ij} correspond to the inequalities.
992	*/
993	static __isl_give isl_basic_set *valid_direction_lp(
994	__isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)
995	{
996	isl_space *space;
997	struct isl_basic_set *lp;
998	unsigned d;
999	int n;
1000	int i, j, k;
1001	isl_size total;
1002
1003	total = isl_basic_set_dim(bset: bset1, type: isl_dim_all);
1004	if (total < 0 || !bset2)
1005	goto error;
1006	d = 1 + total;
1007	n = 2 +
1008	2 * bset1->n_eq + bset1->n_ineq + 2 * bset2->n_eq + bset2->n_ineq;
1009	space = isl_space_set_alloc(ctx: bset1->ctx, nparam: 0, dim: n);
1010	lp = isl_basic_set_alloc_space(space, extra: 0, n_eq: d, n_ineq: n);
1011	if (!lp)
1012	goto error;
1013	for (i = 0; i < n; ++i) {
1014	k = isl_basic_set_alloc_inequality(bset: lp);
1015	if (k < 0)
1016	goto error;
1017	isl_seq_clr(p: lp->ineq[k] + 1, len: n);
1018	isl_int_set_si(lp->ineq[k][0], -1);
1019	isl_int_set_si(lp->ineq[k][1 + i], 1);
1020	}
1021	for (i = 0; i < d; ++i) {
1022	k = isl_basic_set_alloc_equality(bset: lp);
1023	if (k < 0)
1024	goto error;
1025	n = 0;
1026	isl_int_set_si(lp->eq[k][n], 0); n++;
1027	/* positivity constraint 1 >= 0 */
1028	isl_int_set_si(lp->eq[k][n], i == 0); n++;
1029	for (j = 0; j < bset1->n_eq; ++j) {
1030	isl_int_set(lp->eq[k][n], bset1->eq[j][i]); n++;
1031	isl_int_neg(lp->eq[k][n], bset1->eq[j][i]); n++;
1032	}
1033	for (j = 0; j < bset1->n_ineq; ++j) {
1034	isl_int_set(lp->eq[k][n], bset1->ineq[j][i]); n++;
1035	}
1036	/* positivity constraint 1 >= 0 */
1037	isl_int_set_si(lp->eq[k][n], -(i == 0)); n++;
1038	for (j = 0; j < bset2->n_eq; ++j) {
1039	isl_int_neg(lp->eq[k][n], bset2->eq[j][i]); n++;
1040	isl_int_set(lp->eq[k][n], bset2->eq[j][i]); n++;
1041	}
1042	for (j = 0; j < bset2->n_ineq; ++j) {
1043	isl_int_neg(lp->eq[k][n], bset2->ineq[j][i]); n++;
1044	}
1045	}
1046	lp = isl_basic_set_gauss(bset: lp, NULL);
1047	isl_basic_set_free(bset: bset1);
1048	isl_basic_set_free(bset: bset2);
1049	return lp;
1050	error:
1051	isl_basic_set_free(bset: bset1);
1052	isl_basic_set_free(bset: bset2);
1053	return NULL;
1054	}
1055
1056	/* Compute a vector s in the homogeneous space such that <s, r> > 0
1057	* for all rays in the homogeneous space of the two cones that correspond
1058	* to the input polyhedra bset1 and bset2.
1059	*
1060	* We compute s as a vector that satisfies
1061	*
1062	* s = \sum_j \alpha_{ij} h_{ij} for i = 1,2 (*)
1063	*
1064	* with h_{ij} the normals of the facets of polyhedron i
1065	* (including the "positivity constraint" 1 >= 0) and \alpha_{ij}
1066	* strictly positive numbers. For simplicity we impose \alpha_{ij} >= 1.
1067	* We first set up an LP with as variables the \alpha{ij}.
1068	* In this formulation, for each polyhedron i,
1069	* the first constraint is the positivity constraint, followed by pairs
1070	* of variables for the equalities, followed by variables for the inequalities.
1071	* We then simply pick a feasible solution and compute s using (*).
1072	*
1073	* Note that we simply pick any valid direction and make no attempt
1074	* to pick a "good" or even the "best" valid direction.
1075	*/
1076	static __isl_give isl_vec *valid_direction(
1077	__isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)
1078	{
1079	struct isl_basic_set *lp;
1080	struct isl_tab *tab;
1081	struct isl_vec *sample = NULL;
1082	struct isl_vec *dir;
1083	isl_size d;
1084	int i;
1085	int n;
1086
1087	if (!bset1 || !bset2)
1088	goto error;
1089	lp = valid_direction_lp(bset1: isl_basic_set_copy(bset: bset1),
1090	bset2: isl_basic_set_copy(bset: bset2));
1091	tab = isl_tab_from_basic_set(bset: lp, track: 0);
1092	sample = isl_tab_get_sample_value(tab);
1093	isl_tab_free(tab);
1094	isl_basic_set_free(bset: lp);
1095	if (!sample)
1096	goto error;
1097	d = isl_basic_set_dim(bset: bset1, type: isl_dim_all);
1098	if (d < 0)
1099	goto error;
1100	dir = isl_vec_alloc(ctx: bset1->ctx, size: 1 + d);
1101	if (!dir)
1102	goto error;
1103	isl_seq_clr(p: dir->block.data + 1, len: dir->size - 1);
1104	n = 1;
1105	/* positivity constraint 1 >= 0 */
1106	isl_int_set(dir->block.data[0], sample->block.data[n]); n++;
1107	for (i = 0; i < bset1->n_eq; ++i) {
1108	isl_int_sub(sample->block.data[n],
1109	sample->block.data[n], sample->block.data[n+1]);
1110	isl_seq_combine(dst: dir->block.data,
1111	m1: bset1->ctx->one, src1: dir->block.data,
1112	m2: sample->block.data[n], src2: bset1->eq[i], len: 1 + d);
1113
1114	n += 2;
1115	}
1116	for (i = 0; i < bset1->n_ineq; ++i)
1117	isl_seq_combine(dst: dir->block.data,
1118	m1: bset1->ctx->one, src1: dir->block.data,
1119	m2: sample->block.data[n++], src2: bset1->ineq[i], len: 1 + d);
1120	isl_vec_free(vec: sample);
1121	isl_seq_normalize(ctx: bset1->ctx, p: dir->el, len: dir->size);
1122	isl_basic_set_free(bset: bset1);
1123	isl_basic_set_free(bset: bset2);
1124	return dir;
1125	error:
1126	isl_vec_free(vec: sample);
1127	isl_basic_set_free(bset: bset1);
1128	isl_basic_set_free(bset: bset2);
1129	return NULL;
1130	}
1131
1132	/* Given a polyhedron b_i + A_i x >= 0 and a map T = S^{-1},
1133	* compute b_i' + A_i' x' >= 0, with
1134	*
1135	* [ b_i A_i ] [ y' ] [ y' ]
1136	* [ 1 0 ] S^{-1} [ x' ] >= 0 or [ b_i' A_i' ] [ x' ] >= 0
1137	*
1138	* In particular, add the "positivity constraint" and then perform
1139	* the mapping.
1140	*/
1141	static __isl_give isl_basic_set *homogeneous_map(__isl_take isl_basic_set *bset,
1142	__isl_take isl_mat *T)
1143	{
1144	int k;
1145	isl_size total;
1146
1147	total = isl_basic_set_dim(bset, type: isl_dim_all);
1148	if (total < 0)
1149	goto error;
1150	bset = isl_basic_set_extend_constraints(base: bset, n_eq: 0, n_ineq: 1);
1151	k = isl_basic_set_alloc_inequality(bset);
1152	if (k < 0)
1153	goto error;
1154	isl_seq_clr(p: bset->ineq[k] + 1, len: total);
1155	isl_int_set_si(bset->ineq[k][0], 1);
1156	bset = isl_basic_set_preimage(bset, mat: T);
1157	return bset;
1158	error:
1159	isl_mat_free(mat: T);
1160	isl_basic_set_free(bset);
1161	return NULL;
1162	}
1163
1164	/* Compute the convex hull of a pair of basic sets without any parameters or
1165	* integer divisions, where the convex hull is known to be pointed,
1166	* but the basic sets may be unbounded.
1167	*
1168	* We turn this problem into the computation of a convex hull of a pair
1169	* _bounded_ polyhedra by "changing the direction of the homogeneous
1170	* dimension". This idea is due to Matthias Koeppe.
1171	*
1172	* Consider the cones in homogeneous space that correspond to the
1173	* input polyhedra. The rays of these cones are also rays of the
1174	* polyhedra if the coordinate that corresponds to the homogeneous
1175	* dimension is zero. That is, if the inner product of the rays
1176	* with the homogeneous direction is zero.
1177	* The cones in the homogeneous space can also be considered to
1178	* correspond to other pairs of polyhedra by chosing a different
1179	* homogeneous direction. To ensure that both of these polyhedra
1180	* are bounded, we need to make sure that all rays of the cones
1181	* correspond to vertices and not to rays.
1182	* Let s be a direction such that <s, r> > 0 for all rays r of both cones.
1183	* Then using s as a homogeneous direction, we obtain a pair of polytopes.
1184	* The vector s is computed in valid_direction.
1185	*
1186	* Note that we need to consider _all_ rays of the cones and not just
1187	* the rays that correspond to rays in the polyhedra. If we were to
1188	* only consider those rays and turn them into vertices, then we
1189	* may inadvertently turn some vertices into rays.
1190	*
1191	* The standard homogeneous direction is the unit vector in the 0th coordinate.
1192	* We therefore transform the two polyhedra such that the selected
1193	* direction is mapped onto this standard direction and then proceed
1194	* with the normal computation.
1195	* Let S be a non-singular square matrix with s as its first row,
1196	* then we want to map the polyhedra to the space
1197	*
1198	* [ y' ] [ y ] [ y ] [ y' ]
1199	* [ x' ] = S [ x ] i.e., [ x ] = S^{-1} [ x' ]
1200	*
1201	* We take S to be the unimodular completion of s to limit the growth
1202	* of the coefficients in the following computations.
1203	*
1204	* Let b_i + A_i x >= 0 be the constraints of polyhedron i.
1205	* We first move to the homogeneous dimension
1206	*
1207	* b_i y + A_i x >= 0 [ b_i A_i ] [ y ] [ 0 ]
1208	* y >= 0 or [ 1 0 ] [ x ] >= [ 0 ]
1209	*
1210	* Then we change directoin
1211	*
1212	* [ b_i A_i ] [ y' ] [ y' ]
1213	* [ 1 0 ] S^{-1} [ x' ] >= 0 or [ b_i' A_i' ] [ x' ] >= 0
1214	*
1215	* Then we compute the convex hull of the polytopes b_i' + A_i' x' >= 0
1216	* resulting in b' + A' x' >= 0, which we then convert back
1217	*
1218	* [ y ] [ y ]
1219	* [ b' A' ] S [ x ] >= 0 or [ b A ] [ x ] >= 0
1220	*
1221	* The polyhedron b + A x >= 0 is then the convex hull of the input polyhedra.
1222	*/
1223	static __isl_give isl_basic_set *convex_hull_pair_pointed(
1224	__isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)
1225	{
1226	struct isl_ctx *ctx = NULL;
1227	struct isl_vec *dir = NULL;
1228	struct isl_mat *T = NULL;
1229	struct isl_mat *T2 = NULL;
1230	struct isl_basic_set *hull;
1231	struct isl_set *set;
1232
1233	if (!bset1 || !bset2)
1234	goto error;
1235	ctx = isl_basic_set_get_ctx(bset: bset1);
1236	dir = valid_direction(bset1: isl_basic_set_copy(bset: bset1),
1237	bset2: isl_basic_set_copy(bset: bset2));
1238	if (!dir)
1239	goto error;
1240	T = isl_mat_alloc(ctx, n_row: dir->size, n_col: dir->size);
1241	if (!T)
1242	goto error;
1243	isl_seq_cpy(dst: T->row[0], src: dir->block.data, len: dir->size);
1244	T = isl_mat_unimodular_complete(M: T, row: 1);
1245	T2 = isl_mat_right_inverse(mat: isl_mat_copy(mat: T));
1246
1247	bset1 = homogeneous_map(bset: bset1, T: isl_mat_copy(mat: T2));
1248	bset2 = homogeneous_map(bset: bset2, T: T2);
1249	set = isl_set_alloc_space(space: isl_basic_set_get_space(bset: bset1), n: 2, flags: 0);
1250	set = isl_set_add_basic_set(set, bset: bset1);
1251	set = isl_set_add_basic_set(set, bset: bset2);
1252	hull = uset_convex_hull(set);
1253	hull = isl_basic_set_preimage(bset: hull, mat: T);
1254
1255	isl_vec_free(vec: dir);
1256
1257	return hull;
1258	error:
1259	isl_vec_free(vec: dir);
1260	isl_basic_set_free(bset: bset1);
1261	isl_basic_set_free(bset: bset2);
1262	return NULL;
1263	}
1264
1265	static __isl_give isl_basic_set *uset_convex_hull_wrap(__isl_take isl_set *set);
1266	static __isl_give isl_basic_set *modulo_affine_hull(
1267	__isl_take isl_set *set, __isl_take isl_basic_set *affine_hull);
1268
1269	/* Compute the convex hull of a pair of basic sets without any parameters or
1270	* integer divisions.
1271	*
1272	* This function is called from uset_convex_hull_unbounded, which
1273	* means that the complete convex hull is unbounded. Some pairs
1274	* of basic sets may still be bounded, though.
1275	* They may even lie inside a lower dimensional space, in which
1276	* case they need to be handled inside their affine hull since
1277	* the main algorithm assumes that the result is full-dimensional.
1278	*
1279	* If the convex hull of the two basic sets would have a non-trivial
1280	* lineality space, we first project out this lineality space.
1281	*/
1282	static __isl_give isl_basic_set *convex_hull_pair(
1283	__isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)
1284	{
1285	isl_basic_set *lin, *aff;
1286	isl_bool bounded1, bounded2;
1287	isl_size total;
1288
1289	if (bset1->ctx->opt->convex == ISL_CONVEX_HULL_FM)
1290	return convex_hull_pair_elim(bset1, bset2);
1291
1292	aff = isl_set_affine_hull(set: isl_basic_set_union(bset1: isl_basic_set_copy(bset: bset1),
1293	bset2: isl_basic_set_copy(bset: bset2)));
1294	if (!aff)
1295	goto error;
1296	if (aff->n_eq != 0)
1297	return modulo_affine_hull(set: isl_basic_set_union(bset1, bset2), affine_hull: aff);
1298	isl_basic_set_free(bset: aff);
1299
1300	bounded1 = isl_basic_set_is_bounded(bset: bset1);
1301	bounded2 = isl_basic_set_is_bounded(bset: bset2);
1302
1303	if (bounded1 < 0 || bounded2 < 0)
1304	goto error;
1305
1306	if (bounded1 && bounded2)
1307	return uset_convex_hull_wrap(set: isl_basic_set_union(bset1, bset2));
1308
1309	if (bounded1 || bounded2)
1310	return convex_hull_pair_pointed(bset1, bset2);
1311
1312	lin = induced_lineality_space(bset1: isl_basic_set_copy(bset: bset1),
1313	bset2: isl_basic_set_copy(bset: bset2));
1314	if (!lin)
1315	goto error;
1316	if (isl_basic_set_plain_is_universe(bset: lin)) {
1317	isl_basic_set_free(bset: bset1);
1318	isl_basic_set_free(bset: bset2);
1319	return lin;
1320	}
1321	total = isl_basic_set_dim(bset: lin, type: isl_dim_all);
1322	if (lin->n_eq < total) {
1323	struct isl_set *set;
1324	set = isl_set_alloc_space(space: isl_basic_set_get_space(bset: bset1), n: 2, flags: 0);
1325	set = isl_set_add_basic_set(set, bset: bset1);
1326	set = isl_set_add_basic_set(set, bset: bset2);
1327	return modulo_lineality(set, lin);
1328	}
1329	isl_basic_set_free(bset: lin);
1330	if (total < 0)
1331	goto error;
1332
1333	return convex_hull_pair_pointed(bset1, bset2);
1334	error:
1335	isl_basic_set_free(bset: bset1);
1336	isl_basic_set_free(bset: bset2);
1337	return NULL;
1338	}
1339
1340	/* Compute the lineality space of a basic set.
1341	* We basically just drop the constants and turn every inequality
1342	* into an equality.
1343	* Any explicit representations of local variables are removed
1344	* because they may no longer be valid representations
1345	* in the lineality space.
1346	*/
1347	__isl_give isl_basic_set *isl_basic_set_lineality_space(
1348	__isl_take isl_basic_set *bset)
1349	{
1350	int i, k;
1351	struct isl_basic_set *lin = NULL;
1352	isl_size n_div, dim;
1353
1354	n_div = isl_basic_set_dim(bset, type: isl_dim_div);
1355	dim = isl_basic_set_dim(bset, type: isl_dim_all);
1356	if (n_div < 0 || dim < 0)
1357	return isl_basic_set_free(bset);
1358
1359	lin = isl_basic_set_alloc_space(space: isl_basic_set_get_space(bset),
1360	extra: n_div, n_eq: dim, n_ineq: 0);
1361	for (i = 0; i < n_div; ++i)
1362	if (isl_basic_set_alloc_div(bset: lin) < 0)
1363	goto error;
1364	if (!lin)
1365	goto error;
1366	for (i = 0; i < bset->n_eq; ++i) {
1367	k = isl_basic_set_alloc_equality(bset: lin);
1368	if (k < 0)
1369	goto error;
1370	isl_int_set_si(lin->eq[k][0], 0);
1371	isl_seq_cpy(dst: lin->eq[k] + 1, src: bset->eq[i] + 1, len: dim);
1372	}
1373	lin = isl_basic_set_gauss(bset: lin, NULL);
1374	if (!lin)
1375	goto error;
1376	for (i = 0; i < bset->n_ineq && lin->n_eq < dim; ++i) {
1377	k = isl_basic_set_alloc_equality(bset: lin);
1378	if (k < 0)
1379	goto error;
1380	isl_int_set_si(lin->eq[k][0], 0);
1381	isl_seq_cpy(dst: lin->eq[k] + 1, src: bset->ineq[i] + 1, len: dim);
1382	lin = isl_basic_set_gauss(bset: lin, NULL);
1383	if (!lin)
1384	goto error;
1385	}
1386	isl_basic_set_free(bset);
1387	return lin;
1388	error:
1389	isl_basic_set_free(bset: lin);
1390	isl_basic_set_free(bset);
1391	return NULL;
1392	}
1393
1394	/* Compute the (linear) hull of the lineality spaces of the basic sets in the
1395	* set "set".
1396	*/
1397	__isl_give isl_basic_set *isl_set_combined_lineality_space(
1398	__isl_take isl_set *set)
1399	{
1400	int i;
1401	struct isl_set *lin = NULL;
1402
1403	if (!set)
1404	return NULL;
1405	if (set->n == 0) {
1406	isl_space *space = isl_set_get_space(set);
1407	isl_set_free(set);
1408	return isl_basic_set_empty(space);
1409	}
1410
1411	lin = isl_set_alloc_space(space: isl_set_get_space(set), n: set->n, flags: 0);
1412	for (i = 0; i < set->n; ++i)
1413	lin = isl_set_add_basic_set(set: lin,
1414	bset: isl_basic_set_lineality_space(bset: isl_basic_set_copy(bset: set->p[i])));
1415	isl_set_free(set);
1416	return isl_set_affine_hull(set: lin);
1417	}
1418
1419	/* Compute the convex hull of a set without any parameters or
1420	* integer divisions.
1421	* In each step, we combined two basic sets until only one
1422	* basic set is left.
1423	* The input basic sets are assumed not to have a non-trivial
1424	* lineality space. If any of the intermediate results has
1425	* a non-trivial lineality space, it is projected out.
1426	*/
1427	static __isl_give isl_basic_set *uset_convex_hull_unbounded(
1428	__isl_take isl_set *set)
1429	{
1430	isl_basic_set_list *list;
1431
1432	list = isl_set_get_basic_set_list(set);
1433	isl_set_free(set);
1434
1435	while (list) {
1436	isl_size n, total;
1437	struct isl_basic_set *t;
1438	isl_basic_set *bset1, *bset2;
1439
1440	n = isl_basic_set_list_n_basic_set(list);
1441	if (n < 0)
1442	goto error;
1443	if (n < 2)
1444	isl_die(isl_basic_set_list_get_ctx(list),
1445	isl_error_internal,
1446	"expecting at least two elements", goto error);
1447	bset1 = isl_basic_set_list_get_basic_set(list, index: n - 1);
1448	bset2 = isl_basic_set_list_get_basic_set(list, index: n - 2);
1449	bset1 = convex_hull_pair(bset1, bset2);
1450	if (n == 2) {
1451	isl_basic_set_list_free(list);
1452	return bset1;
1453	}
1454	bset1 = isl_basic_set_underlying_set(bset: bset1);
1455	list = isl_basic_set_list_drop(list, first: n - 2, n: 2);
1456	list = isl_basic_set_list_add(list, el: bset1);
1457
1458	t = isl_basic_set_list_get_basic_set(list, index: n - 2);
1459	t = isl_basic_set_lineality_space(bset: t);
1460	if (!t)
1461	goto error;
1462	if (isl_basic_set_plain_is_universe(bset: t)) {
1463	isl_basic_set_list_free(list);
1464	return t;
1465	}
1466	total = isl_basic_set_dim(bset: t, type: isl_dim_all);
1467	if (t->n_eq < total) {
1468	set = isl_basic_set_list_union(list);
1469	return modulo_lineality(set, lin: t);
1470	}
1471	isl_basic_set_free(bset: t);
1472	if (total < 0)
1473	goto error;
1474	}
1475
1476	return NULL;
1477	error:
1478	isl_basic_set_list_free(list);
1479	return NULL;
1480	}
1481
1482	/* Compute an initial hull for wrapping containing a single initial
1483	* facet.
1484	* This function assumes that the given set is bounded.
1485	*/
1486	static __isl_give isl_basic_set *initial_hull(__isl_take isl_basic_set *hull,
1487	__isl_keep isl_set *set)
1488	{
1489	struct isl_mat *bounds = NULL;
1490	isl_size dim;
1491	int k;
1492
1493	if (!hull)
1494	goto error;
1495	bounds = initial_facet_constraint(set);
1496	if (!bounds)
1497	goto error;
1498	k = isl_basic_set_alloc_inequality(bset: hull);
1499	if (k < 0)
1500	goto error;
1501	dim = isl_set_dim(set, type: isl_dim_set);
1502	if (dim < 0)
1503	goto error;
1504	isl_assert(set->ctx, 1 + dim == bounds->n_col, goto error);
1505	isl_seq_cpy(dst: hull->ineq[k], src: bounds->row[0], len: bounds->n_col);
1506	isl_mat_free(mat: bounds);
1507
1508	return hull;
1509	error:
1510	isl_basic_set_free(bset: hull);
1511	isl_mat_free(mat: bounds);
1512	return NULL;
1513	}
1514
1515	struct max_constraint {
1516	struct isl_mat *c;
1517	int count;
1518	int ineq;
1519	};
1520
1521	static isl_bool max_constraint_equal(const void *entry, const void *val)
1522	{
1523	struct max_constraint *a = (struct max_constraint *)entry;
1524	isl_int *b = (isl_int *)val;
1525
1526	return isl_bool_ok(b: isl_seq_eq(p1: a->c->row[0] + 1, p2: b, len: a->c->n_col - 1));
1527	}
1528
1529	static isl_stat update_constraint(struct isl_ctx *ctx,
1530	struct isl_hash_table *table,
1531	isl_int *con, unsigned len, int n, int ineq)
1532	{
1533	struct isl_hash_table_entry *entry;
1534	struct max_constraint *c;
1535	uint32_t c_hash;
1536
1537	c_hash = isl_seq_get_hash(p: con + 1, len);
1538	entry = isl_hash_table_find(ctx, table, key_hash: c_hash, eq: max_constraint_equal,
1539	val: con + 1, reserve: 0);
1540	if (!entry)
1541	return isl_stat_error;
1542	if (entry == isl_hash_table_entry_none)
1543	return isl_stat_ok;
1544	c = entry->data;
1545	if (c->count < n) {
1546	isl_hash_table_remove(ctx, table, entry);
1547	return isl_stat_ok;
1548	}
1549	c->count++;
1550	if (isl_int_gt(c->c->row[0][0], con[0]))
1551	return isl_stat_ok;
1552	if (isl_int_eq(c->c->row[0][0], con[0])) {
1553	if (ineq)
1554	c->ineq = ineq;
1555	return isl_stat_ok;
1556	}
1557	c->c = isl_mat_cow(mat: c->c);
1558	isl_int_set(c->c->row[0][0], con[0]);
1559	c->ineq = ineq;
1560
1561	return isl_stat_ok;
1562	}
1563
1564	/* Check whether the constraint hash table "table" contains the constraint
1565	* "con".
1566	*/
1567	static isl_bool has_constraint(struct isl_ctx *ctx,
1568	struct isl_hash_table *table, isl_int *con, unsigned len, int n)
1569	{
1570	struct isl_hash_table_entry *entry;
1571	struct max_constraint *c;
1572	uint32_t c_hash;
1573
1574	c_hash = isl_seq_get_hash(p: con + 1, len);
1575	entry = isl_hash_table_find(ctx, table, key_hash: c_hash, eq: max_constraint_equal,
1576	val: con + 1, reserve: 0);
1577	if (!entry)
1578	return isl_bool_error;
1579	if (entry == isl_hash_table_entry_none)
1580	return isl_bool_false;
1581	c = entry->data;
1582	if (c->count < n)
1583	return isl_bool_false;
1584	return isl_bool_ok(isl_int_eq(c->c->row[0][0], con[0]));
1585	}
1586
1587	/* Are the constraints of "bset" known to be facets?
1588	* If there are any equality constraints, then they are not.
1589	* If there may be redundant constraints, then those
1590	* redundant constraints are not facets.
1591	*/
1592	static isl_bool has_facets(__isl_keep isl_basic_set *bset)
1593	{
1594	isl_size n_eq;
1595
1596	n_eq = isl_basic_set_n_equality(bset);
1597	if (n_eq < 0)
1598	return isl_bool_error;
1599	if (n_eq != 0)
1600	return isl_bool_false;
1601	return ISL_F_ISSET(bset, ISL_BASIC_SET_NO_REDUNDANT);
1602	}
1603
1604	/* Check for inequality constraints of a basic set without equalities
1605	* or redundant constraints
1606	* such that the same or more stringent copies of the constraint appear
1607	* in all of the basic sets. Such constraints are necessarily facet
1608	* constraints of the convex hull.
1609	*
1610	* If the resulting basic set is by chance identical to one of
1611	* the basic sets in "set", then we know that this basic set contains
1612	* all other basic sets and is therefore the convex hull of set.
1613	* In this case we set *is_hull to 1.
1614	*/
1615	static __isl_give isl_basic_set *common_constraints(
1616	__isl_take isl_basic_set *hull, __isl_keep isl_set *set, int *is_hull)
1617	{
1618	int i, j, s, n;
1619	int min_constraints;
1620	int best;
1621	struct max_constraint *constraints = NULL;
1622	struct isl_hash_table *table = NULL;
1623	isl_size total;
1624
1625	*is_hull = 0;
1626
1627	for (i = 0; i < set->n; ++i) {
1628	isl_bool facets = has_facets(bset: set->p[i]);
1629	if (facets < 0)
1630	return isl_basic_set_free(bset: hull);
1631	if (facets)
1632	break;
1633	}
1634	if (i >= set->n)
1635	return hull;
1636	min_constraints = set->p[i]->n_ineq;
1637	best = i;
1638	for (i = best + 1; i < set->n; ++i) {
1639	isl_bool facets = has_facets(bset: set->p[i]);
1640	if (facets < 0)
1641	return isl_basic_set_free(bset: hull);
1642	if (!facets)
1643	continue;
1644	if (set->p[i]->n_ineq >= min_constraints)
1645	continue;
1646	min_constraints = set->p[i]->n_ineq;
1647	best = i;
1648	}
1649	constraints = isl_calloc_array(hull->ctx, struct max_constraint,
1650	min_constraints);
1651	if (!constraints)
1652	return hull;
1653	table = isl_alloc_type(hull->ctx, struct isl_hash_table);
1654	if (isl_hash_table_init(ctx: hull->ctx, table, min_size: min_constraints))
1655	goto error;
1656
1657	total = isl_set_dim(set, type: isl_dim_all);
1658	if (total < 0)
1659	goto error;
1660	for (i = 0; i < set->p[best]->n_ineq; ++i) {
1661	constraints[i].c = isl_mat_sub_alloc6(ctx: hull->ctx,
1662	row: set->p[best]->ineq + i, first_row: 0, n_row: 1, first_col: 0, n_col: 1 + total);
1663	if (!constraints[i].c)
1664	goto error;
1665	constraints[i].ineq = 1;
1666	}
1667	for (i = 0; i < min_constraints; ++i) {
1668	struct isl_hash_table_entry *entry;
1669	uint32_t c_hash;
1670	c_hash = isl_seq_get_hash(p: constraints[i].c->row[0] + 1, len: total);
1671	entry = isl_hash_table_find(ctx: hull->ctx, table, key_hash: c_hash,
1672	eq: max_constraint_equal, val: constraints[i].c->row[0] + 1, reserve: 1);
1673	if (!entry)
1674	goto error;
1675	isl_assert(hull->ctx, !entry->data, goto error);
1676	entry->data = &constraints[i];
1677	}
1678
1679	n = 0;
1680	for (s = 0; s < set->n; ++s) {
1681	if (s == best)
1682	continue;
1683
1684	for (i = 0; i < set->p[s]->n_eq; ++i) {
1685	isl_int *eq = set->p[s]->eq[i];
1686	for (j = 0; j < 2; ++j) {
1687	isl_seq_neg(dst: eq, src: eq, len: 1 + total);
1688	if (update_constraint(ctx: hull->ctx, table,
1689	con: eq, len: total, n, ineq: 0) < 0)
1690	goto error;
1691	}
1692	}
1693	for (i = 0; i < set->p[s]->n_ineq; ++i) {
1694	isl_int *ineq = set->p[s]->ineq[i];
1695	if (update_constraint(ctx: hull->ctx, table, con: ineq, len: total, n,
1696	ineq: set->p[s]->n_eq == 0) < 0)
1697	goto error;
1698	}
1699	++n;
1700	}
1701
1702	for (i = 0; i < min_constraints; ++i) {
1703	if (constraints[i].count < n)
1704	continue;
1705	if (!constraints[i].ineq)
1706	continue;
1707	j = isl_basic_set_alloc_inequality(bset: hull);
1708	if (j < 0)
1709	goto error;
1710	isl_seq_cpy(dst: hull->ineq[j], src: constraints[i].c->row[0], len: 1 + total);
1711	}
1712
1713	for (s = 0; s < set->n; ++s) {
1714	if (set->p[s]->n_eq)
1715	continue;
1716	if (set->p[s]->n_ineq != hull->n_ineq)
1717	continue;
1718	for (i = 0; i < set->p[s]->n_ineq; ++i) {
1719	isl_bool has;
1720	isl_int *ineq = set->p[s]->ineq[i];
1721	has = has_constraint(ctx: hull->ctx, table, con: ineq, len: total, n);
1722	if (has < 0)
1723	goto error;
1724	if (!has)
1725	break;
1726	}
1727	if (i == set->p[s]->n_ineq)
1728	*is_hull = 1;
1729	}
1730
1731	isl_hash_table_clear(table);
1732	for (i = 0; i < min_constraints; ++i)
1733	isl_mat_free(mat: constraints[i].c);
1734	free(ptr: constraints);
1735	free(ptr: table);
1736	return hull;
1737	error:
1738	isl_hash_table_clear(table);
1739	free(ptr: table);
1740	if (constraints)
1741	for (i = 0; i < min_constraints; ++i)
1742	isl_mat_free(mat: constraints[i].c);
1743	free(ptr: constraints);
1744	return hull;
1745	}
1746
1747	/* Create a template for the convex hull of "set" and fill it up
1748	* obvious facet constraints, if any. If the result happens to
1749	* be the convex hull of "set" then *is_hull is set to 1.
1750	*/
1751	static __isl_give isl_basic_set *proto_hull(__isl_keep isl_set *set,
1752	int *is_hull)
1753	{
1754	struct isl_basic_set *hull;
1755	unsigned n_ineq;
1756	int i;
1757
1758	n_ineq = 1;
1759	for (i = 0; i < set->n; ++i) {
1760	n_ineq += set->p[i]->n_eq;
1761	n_ineq += set->p[i]->n_ineq;
1762	}
1763	hull = isl_basic_set_alloc_space(space: isl_space_copy(space: set->dim), extra: 0, n_eq: 0, n_ineq);
1764	hull = isl_basic_set_set_rational(bset: hull);
1765	if (!hull)
1766	return NULL;
1767	return common_constraints(hull, set, is_hull);
1768	}
1769
1770	static __isl_give isl_basic_set *uset_convex_hull_wrap(__isl_take isl_set *set)
1771	{
1772	struct isl_basic_set *hull;
1773	int is_hull;
1774
1775	hull = proto_hull(set, is_hull: &is_hull);
1776	if (hull && !is_hull) {
1777	if (hull->n_ineq == 0)
1778	hull = initial_hull(hull, set);
1779	hull = extend(hull, set);
1780	}
1781	isl_set_free(set);
1782
1783	return hull;
1784	}
1785
1786	/* Compute the convex hull of a set without any parameters or
1787	* integer divisions. Depending on whether the set is bounded,
1788	* we pass control to the wrapping based convex hull or
1789	* the Fourier-Motzkin elimination based convex hull.
1790	* We also handle a few special cases before checking the boundedness.
1791	*/
1792	static __isl_give isl_basic_set *uset_convex_hull(__isl_take isl_set *set)
1793	{
1794	isl_bool bounded;
1795	isl_size dim;
1796	struct isl_basic_set *convex_hull = NULL;
1797	struct isl_basic_set *lin;
1798
1799	dim = isl_set_dim(set, type: isl_dim_all);
1800	if (dim < 0)
1801	goto error;
1802	if (dim == 0)
1803	return convex_hull_0d(set);
1804
1805	set = isl_set_coalesce(set);
1806	set = isl_set_set_rational(set);
1807
1808	if (!set)
1809	return NULL;
1810	if (set->n == 1) {
1811	convex_hull = isl_basic_set_copy(bset: set->p[0]);
1812	isl_set_free(set);
1813	return convex_hull;
1814	}
1815	if (dim == 1)
1816	return convex_hull_1d(set);
1817
1818	bounded = isl_set_is_bounded(set);
1819	if (bounded < 0)
1820	goto error;
1821	if (bounded && set->ctx->opt->convex == ISL_CONVEX_HULL_WRAP)
1822	return uset_convex_hull_wrap(set);
1823
1824	lin = isl_set_combined_lineality_space(set: isl_set_copy(set));
1825	if (!lin)
1826	goto error;
1827	if (isl_basic_set_plain_is_universe(bset: lin)) {
1828	isl_set_free(set);
1829	return lin;
1830	}
1831	if (lin->n_eq < dim)
1832	return modulo_lineality(set, lin);
1833	isl_basic_set_free(bset: lin);
1834
1835	return uset_convex_hull_unbounded(set);
1836	error:
1837	isl_set_free(set);
1838	isl_basic_set_free(bset: convex_hull);
1839	return NULL;
1840	}
1841
1842	/* This is the core procedure, where "set" is a "pure" set, i.e.,
1843	* without parameters or divs and where the convex hull of set is
1844	* known to be full-dimensional.
1845	*/
1846	static __isl_give isl_basic_set *uset_convex_hull_wrap_bounded(
1847	__isl_take isl_set *set)
1848	{
1849	struct isl_basic_set *convex_hull = NULL;
1850	isl_size dim;
1851
1852	dim = isl_set_dim(set, type: isl_dim_all);
1853	if (dim < 0)
1854	goto error;
1855
1856	if (dim == 0) {
1857	convex_hull = isl_basic_set_universe(space: isl_space_copy(space: set->dim));
1858	isl_set_free(set);
1859	convex_hull = isl_basic_set_set_rational(bset: convex_hull);
1860	return convex_hull;
1861	}
1862
1863	set = isl_set_set_rational(set);
1864	set = isl_set_coalesce(set);
1865	if (!set)
1866	goto error;
1867	if (set->n == 1) {
1868	convex_hull = isl_basic_set_copy(bset: set->p[0]);
1869	isl_set_free(set);
1870	convex_hull = isl_basic_map_remove_redundancies(bmap: convex_hull);
1871	return convex_hull;
1872	}
1873	if (dim == 1)
1874	return convex_hull_1d(set);
1875
1876	return uset_convex_hull_wrap(set);
1877	error:
1878	isl_set_free(set);
1879	return NULL;
1880	}
1881
1882	/* Compute the convex hull of set "set" with affine hull "affine_hull",
1883	* We first remove the equalities (transforming the set), compute the
1884	* convex hull of the transformed set and then add the equalities back
1885	* (after performing the inverse transformation.
1886	*/
1887	static __isl_give isl_basic_set *modulo_affine_hull(
1888	__isl_take isl_set *set, __isl_take isl_basic_set *affine_hull)
1889	{
1890	struct isl_mat *T;
1891	struct isl_mat *T2;
1892	struct isl_basic_set *dummy;
1893	struct isl_basic_set *convex_hull;
1894
1895	dummy = isl_basic_set_remove_equalities(
1896	bset: isl_basic_set_copy(bset: affine_hull), T: &T, T2: &T2);
1897	if (!dummy)
1898	goto error;
1899	isl_basic_set_free(bset: dummy);
1900	set = isl_set_preimage(set, mat: T);
1901	convex_hull = uset_convex_hull(set);
1902	convex_hull = isl_basic_set_preimage(bset: convex_hull, mat: T2);
1903	convex_hull = isl_basic_set_intersect(bset1: convex_hull, bset2: affine_hull);
1904	return convex_hull;
1905	error:
1906	isl_mat_free(mat: T);
1907	isl_mat_free(mat: T2);
1908	isl_basic_set_free(bset: affine_hull);
1909	isl_set_free(set);
1910	return NULL;
1911	}
1912
1913	/* Return an empty basic map living in the same space as "map".
1914	*/
1915	static __isl_give isl_basic_map *replace_map_by_empty_basic_map(
1916	__isl_take isl_map *map)
1917	{
1918	isl_space *space;
1919
1920	space = isl_map_get_space(map);
1921	isl_map_free(map);
1922	return isl_basic_map_empty(space);
1923	}
1924
1925	/* Compute the convex hull of a map.
1926	*
1927	* The implementation was inspired by "Extended Convex Hull" by Fukuda et al.,
1928	* specifically, the wrapping of facets to obtain new facets.
1929	*/
1930	__isl_give isl_basic_map *isl_map_convex_hull(__isl_take isl_map *map)
1931	{
1932	struct isl_basic_set *bset;
1933	struct isl_basic_map *model = NULL;
1934	struct isl_basic_set *affine_hull = NULL;
1935	struct isl_basic_map *convex_hull = NULL;
1936	struct isl_set *set = NULL;
1937
1938	map = isl_map_detect_equalities(map);
1939	map = isl_map_align_divs_internal(map);
1940	if (!map)
1941	goto error;
1942
1943	if (map->n == 0)
1944	return replace_map_by_empty_basic_map(map);
1945
1946	model = isl_basic_map_copy(bmap: map->p[0]);
1947	set = isl_map_underlying_set(map);
1948	if (!set)
1949	goto error;
1950
1951	affine_hull = isl_set_affine_hull(set: isl_set_copy(set));
1952	if (!affine_hull)
1953	goto error;
1954	if (affine_hull->n_eq != 0)
1955	bset = modulo_affine_hull(set, affine_hull);
1956	else {
1957	isl_basic_set_free(bset: affine_hull);
1958	bset = uset_convex_hull(set);
1959	}
1960
1961	convex_hull = isl_basic_map_overlying_set(bset, like: model);
1962	if (!convex_hull)
1963	return NULL;
1964
1965	ISL_F_SET(convex_hull, ISL_BASIC_MAP_NO_IMPLICIT);
1966	ISL_F_SET(convex_hull, ISL_BASIC_MAP_ALL_EQUALITIES);
1967	ISL_F_CLR(convex_hull, ISL_BASIC_MAP_RATIONAL);
1968	return convex_hull;
1969	error:
1970	isl_set_free(set);
1971	isl_basic_map_free(bmap: model);
1972	return NULL;
1973	}
1974
1975	__isl_give isl_basic_set *isl_set_convex_hull(__isl_take isl_set *set)
1976	{
1977	return bset_from_bmap(bmap: isl_map_convex_hull(map: set_to_map(set)));
1978	}
1979
1980	__isl_give isl_basic_map *isl_map_polyhedral_hull(__isl_take isl_map *map)
1981	{
1982	isl_basic_map *hull;
1983
1984	hull = isl_map_convex_hull(map);
1985	return isl_basic_map_remove_divs(bmap: hull);
1986	}
1987
1988	__isl_give isl_basic_set *isl_set_polyhedral_hull(__isl_take isl_set *set)
1989	{
1990	return bset_from_bmap(bmap: isl_map_polyhedral_hull(map: set_to_map(set)));
1991	}
1992
1993	struct sh_data_entry {
1994	struct isl_hash_table *table;
1995	struct isl_tab *tab;
1996	};
1997
1998	/* Holds the data needed during the simple hull computation.
1999	* In particular,
2000	* n the number of basic sets in the original set
2001	* hull_table a hash table of already computed constraints
2002	* in the simple hull
2003	* p for each basic set,
2004	* table a hash table of the constraints
2005	* tab the tableau corresponding to the basic set
2006	*/
2007	struct sh_data {
2008	struct isl_ctx *ctx;
2009	unsigned n;
2010	struct isl_hash_table *hull_table;
2011	struct sh_data_entry p[1];
2012	};
2013
2014	static void sh_data_free(struct sh_data *data)
2015	{
2016	int i;
2017
2018	if (!data)
2019	return;
2020	isl_hash_table_free(ctx: data->ctx, table: data->hull_table);
2021	for (i = 0; i < data->n; ++i) {
2022	isl_hash_table_free(ctx: data->ctx, table: data->p[i].table);
2023	isl_tab_free(tab: data->p[i].tab);
2024	}
2025	free(ptr: data);
2026	}
2027
2028	struct ineq_cmp_data {
2029	unsigned len;
2030	isl_int *p;
2031	};
2032
2033	static isl_bool has_ineq(const void *entry, const void *val)
2034	{
2035	isl_int *row = (isl_int *)entry;
2036	struct ineq_cmp_data *v = (struct ineq_cmp_data *)val;
2037
2038	return isl_bool_ok(b: isl_seq_eq(p1: row + 1, p2: v->p + 1, len: v->len) ||
2039	isl_seq_is_neg(p1: row + 1, p2: v->p + 1, len: v->len));
2040	}
2041
2042	static int hash_ineq(struct isl_ctx *ctx, struct isl_hash_table *table,
2043	isl_int *ineq, unsigned len)
2044	{
2045	uint32_t c_hash;
2046	struct ineq_cmp_data v;
2047	struct isl_hash_table_entry *entry;
2048
2049	v.len = len;
2050	v.p = ineq;
2051	c_hash = isl_seq_get_hash(p: ineq + 1, len);
2052	entry = isl_hash_table_find(ctx, table, key_hash: c_hash, eq: has_ineq, val: &v, reserve: 1);
2053	if (!entry)
2054	return - 1;
2055	entry->data = ineq;
2056	return 0;
2057	}
2058
2059	/* Fill hash table "table" with the constraints of "bset".
2060	* Equalities are added as two inequalities.
2061	* The value in the hash table is a pointer to the (in)equality of "bset".
2062	*/
2063	static isl_stat hash_basic_set(struct isl_hash_table *table,
2064	__isl_keep isl_basic_set *bset)
2065	{
2066	int i, j;
2067	isl_size dim = isl_basic_set_dim(bset, type: isl_dim_all);
2068
2069	if (dim < 0)
2070	return isl_stat_error;
2071	for (i = 0; i < bset->n_eq; ++i) {
2072	for (j = 0; j < 2; ++j) {
2073	isl_seq_neg(dst: bset->eq[i], src: bset->eq[i], len: 1 + dim);
2074	if (hash_ineq(ctx: bset->ctx, table, ineq: bset->eq[i], len: dim) < 0)
2075	return isl_stat_error;
2076	}
2077	}
2078	for (i = 0; i < bset->n_ineq; ++i) {
2079	if (hash_ineq(ctx: bset->ctx, table, ineq: bset->ineq[i], len: dim) < 0)
2080	return isl_stat_error;
2081	}
2082	return isl_stat_ok;
2083	}
2084
2085	static struct sh_data *sh_data_alloc(__isl_keep isl_set *set, unsigned n_ineq)
2086	{
2087	struct sh_data *data;
2088	int i;
2089
2090	data = isl_calloc(set->ctx, struct sh_data,
2091	sizeof(struct sh_data) +
2092	(set->n - 1) * sizeof(struct sh_data_entry));
2093	if (!data)
2094	return NULL;
2095	data->ctx = set->ctx;
2096	data->n = set->n;
2097	data->hull_table = isl_hash_table_alloc(ctx: set->ctx, min_size: n_ineq);
2098	if (!data->hull_table)
2099	goto error;
2100	for (i = 0; i < set->n; ++i) {
2101	data->p[i].table = isl_hash_table_alloc(ctx: set->ctx,
2102	min_size: 2 * set->p[i]->n_eq + set->p[i]->n_ineq);
2103	if (!data->p[i].table)
2104	goto error;
2105	if (hash_basic_set(table: data->p[i].table, bset: set->p[i]) < 0)
2106	goto error;
2107	}
2108	return data;
2109	error:
2110	sh_data_free(data);
2111	return NULL;
2112	}
2113
2114	/* Check if inequality "ineq" is a bound for basic set "j" or if
2115	* it can be relaxed (by increasing the constant term) to become
2116	* a bound for that basic set. In the latter case, the constant
2117	* term is updated.
2118	* Relaxation of the constant term is only allowed if "shift" is set.
2119	*
2120	* Return 1 if "ineq" is a bound
2121	* 0 if "ineq" may attain arbitrarily small values on basic set "j"
2122	* -1 if some error occurred
2123	*/
2124	static int is_bound(struct sh_data *data, __isl_keep isl_set *set, int j,
2125	isl_int *ineq, int shift)
2126	{
2127	enum isl_lp_result res;
2128	isl_int opt;
2129
2130	if (!data->p[j].tab) {
2131	data->p[j].tab = isl_tab_from_basic_set(bset: set->p[j], track: 0);
2132	if (!data->p[j].tab)
2133	return -1;
2134	}
2135
2136	isl_int_init(opt);
2137
2138	res = isl_tab_min(tab: data->p[j].tab, f: ineq, denom: data->ctx->one,
2139	opt: &opt, NULL, flags: 0);
2140	if (res == isl_lp_ok && isl_int_is_neg(opt)) {
2141	if (shift)
2142	isl_int_sub(ineq[0], ineq[0], opt);
2143	else
2144	res = isl_lp_unbounded;
2145	}
2146
2147	isl_int_clear(opt);
2148
2149	return (res == isl_lp_ok || res == isl_lp_empty) ? 1 :
2150	res == isl_lp_unbounded ? 0 : -1;
2151	}
2152
2153	/* Set the constant term of "ineq" to the maximum of those of the constraints
2154	* in the basic sets of "set" following "i" that are parallel to "ineq".
2155	* That is, if any of the basic sets of "set" following "i" have a more
2156	* relaxed copy of "ineq", then replace "ineq" by the most relaxed copy.
2157	* "c_hash" is the hash value of the linear part of "ineq".
2158	* "v" has been set up for use by has_ineq.
2159	*
2160	* Note that the two inequality constraints corresponding to an equality are
2161	* represented by the same inequality constraint in data->p[j].table
2162	* (but with different hash values). This means the constraint (or at
2163	* least its constant term) may need to be temporarily negated to get
2164	* the actually hashed constraint.
2165	*/
2166	static isl_stat set_max_constant_term(struct sh_data *data,
2167	__isl_keep isl_set *set,
2168	int i, isl_int *ineq, uint32_t c_hash, struct ineq_cmp_data *v)
2169	{
2170	int j;
2171	isl_ctx *ctx;
2172	struct isl_hash_table_entry *entry;
2173
2174	ctx = isl_set_get_ctx(set);
2175	for (j = i + 1; j < set->n; ++j) {
2176	int neg;
2177	isl_int *ineq_j;
2178
2179	entry = isl_hash_table_find(ctx, table: data->p[j].table,
2180	key_hash: c_hash, eq: &has_ineq, val: v, reserve: 0);
2181	if (!entry)
2182	return isl_stat_error;
2183	if (entry == isl_hash_table_entry_none)
2184	continue;
2185
2186	ineq_j = entry->data;
2187	neg = isl_seq_is_neg(p1: ineq_j + 1, p2: ineq + 1, len: v->len);
2188	if (neg)
2189	isl_int_neg(ineq_j[0], ineq_j[0]);
2190	if (isl_int_gt(ineq_j[0], ineq[0]))
2191	isl_int_set(ineq[0], ineq_j[0]);
2192	if (neg)
2193	isl_int_neg(ineq_j[0], ineq_j[0]);
2194	}
2195
2196	return isl_stat_ok;
2197	}
2198
2199	/* Check if inequality "ineq" from basic set "i" is or can be relaxed to
2200	* become a bound on the whole set. If so, add the (relaxed) inequality
2201	* to "hull". Relaxation is only allowed if "shift" is set.
2202	*
2203	* We first check if "hull" already contains a translate of the inequality.
2204	* If so, we are done.
2205	* Then, we check if any of the previous basic sets contains a translate
2206	* of the inequality. If so, then we have already considered this
2207	* inequality and we are done.
2208	* Otherwise, for each basic set other than "i", we check if the inequality
2209	* is a bound on the basic set, but first replace the constant term
2210	* by the maximal value of any translate of the inequality in any
2211	* of the following basic sets.
2212	* For previous basic sets, we know that they do not contain a translate
2213	* of the inequality, so we directly call is_bound.
2214	* For following basic sets, we first check if a translate of the
2215	* inequality appears in its description. If so, the constant term
2216	* of the inequality has already been updated with respect to this
2217	* translate and the inequality is therefore known to be a bound
2218	* of this basic set.
2219	*/
2220	static __isl_give isl_basic_set *add_bound(__isl_take isl_basic_set *hull,
2221	struct sh_data *data, __isl_keep isl_set *set, int i, isl_int *ineq,
2222	int shift)
2223	{
2224	uint32_t c_hash;
2225	struct ineq_cmp_data v;
2226	struct isl_hash_table_entry *entry;
2227	int j, k;
2228	isl_size total;
2229
2230	total = isl_basic_set_dim(bset: hull, type: isl_dim_all);
2231	if (total < 0)
2232	return isl_basic_set_free(bset: hull);
2233
2234	v.len = total;
2235	v.p = ineq;
2236	c_hash = isl_seq_get_hash(p: ineq + 1, len: v.len);
2237
2238	entry = isl_hash_table_find(ctx: hull->ctx, table: data->hull_table, key_hash: c_hash,
2239	eq: has_ineq, val: &v, reserve: 0);
2240	if (!entry)
2241	return isl_basic_set_free(bset: hull);
2242	if (entry != isl_hash_table_entry_none)
2243	return hull;
2244
2245	for (j = 0; j < i; ++j) {
2246	entry = isl_hash_table_find(ctx: hull->ctx, table: data->p[j].table,
2247	key_hash: c_hash, eq: has_ineq, val: &v, reserve: 0);
2248	if (!entry)
2249	return isl_basic_set_free(bset: hull);
2250	if (entry != isl_hash_table_entry_none)
2251	break;
2252	}
2253	if (j < i)
2254	return hull;
2255
2256	k = isl_basic_set_alloc_inequality(bset: hull);
2257	if (k < 0)
2258	goto error;
2259	isl_seq_cpy(dst: hull->ineq[k], src: ineq, len: 1 + v.len);
2260
2261	if (set_max_constant_term(data, set, i, ineq: hull->ineq[k], c_hash, v: &v) < 0)
2262	goto error;
2263	for (j = 0; j < i; ++j) {
2264	int bound;
2265	bound = is_bound(data, set, j, ineq: hull->ineq[k], shift);
2266	if (bound < 0)
2267	goto error;
2268	if (!bound)
2269	break;
2270	}
2271	if (j < i)
2272	return isl_basic_set_free_inequality(bset: hull, n: 1);
2273
2274	for (j = i + 1; j < set->n; ++j) {
2275	int bound;
2276	entry = isl_hash_table_find(ctx: hull->ctx, table: data->p[j].table,
2277	key_hash: c_hash, eq: has_ineq, val: &v, reserve: 0);
2278	if (!entry)
2279	return isl_basic_set_free(bset: hull);
2280	if (entry != isl_hash_table_entry_none)
2281	continue;
2282	bound = is_bound(data, set, j, ineq: hull->ineq[k], shift);
2283	if (bound < 0)
2284	goto error;
2285	if (!bound)
2286	break;
2287	}
2288	if (j < set->n)
2289	return isl_basic_set_free_inequality(bset: hull, n: 1);
2290
2291	entry = isl_hash_table_find(ctx: hull->ctx, table: data->hull_table, key_hash: c_hash,
2292	eq: has_ineq, val: &v, reserve: 1);
2293	if (!entry)
2294	goto error;
2295	entry->data = hull->ineq[k];
2296
2297	return hull;
2298	error:
2299	isl_basic_set_free(bset: hull);
2300	return NULL;
2301	}
2302
2303	/* Check if any inequality from basic set "i" is or can be relaxed to
2304	* become a bound on the whole set. If so, add the (relaxed) inequality
2305	* to "hull". Relaxation is only allowed if "shift" is set.
2306	*/
2307	static __isl_give isl_basic_set *add_bounds(__isl_take isl_basic_set *bset,
2308	struct sh_data *data, __isl_keep isl_set *set, int i, int shift)
2309	{
2310	int j, k;
2311	isl_size dim = isl_basic_set_dim(bset, type: isl_dim_all);
2312
2313	if (dim < 0)
2314	return isl_basic_set_free(bset);
2315
2316	for (j = 0; j < set->p[i]->n_eq; ++j) {
2317	for (k = 0; k < 2; ++k) {
2318	isl_seq_neg(dst: set->p[i]->eq[j], src: set->p[i]->eq[j], len: 1+dim);
2319	bset = add_bound(hull: bset, data, set, i, ineq: set->p[i]->eq[j],
2320	shift);
2321	}
2322	}
2323	for (j = 0; j < set->p[i]->n_ineq; ++j)
2324	bset = add_bound(hull: bset, data, set, i, ineq: set->p[i]->ineq[j], shift);
2325	return bset;
2326	}
2327
2328	/* Compute a superset of the convex hull of set that is described
2329	* by only (translates of) the constraints in the constituents of set.
2330	* Translation is only allowed if "shift" is set.
2331	*/
2332	static __isl_give isl_basic_set *uset_simple_hull(__isl_take isl_set *set,
2333	int shift)
2334	{
2335	struct sh_data *data = NULL;
2336	struct isl_basic_set *hull = NULL;
2337	unsigned n_ineq;
2338	int i;
2339
2340	if (!set)
2341	return NULL;
2342
2343	n_ineq = 0;
2344	for (i = 0; i < set->n; ++i) {
2345	if (!set->p[i])
2346	goto error;
2347	n_ineq += 2 * set->p[i]->n_eq + set->p[i]->n_ineq;
2348	}
2349
2350	hull = isl_basic_set_alloc_space(space: isl_space_copy(space: set->dim), extra: 0, n_eq: 0, n_ineq);
2351	if (!hull)
2352	goto error;
2353
2354	data = sh_data_alloc(set, n_ineq);
2355	if (!data)
2356	goto error;
2357
2358	for (i = 0; i < set->n; ++i)
2359	hull = add_bounds(bset: hull, data, set, i, shift);
2360
2361	sh_data_free(data);
2362	isl_set_free(set);
2363
2364	return hull;
2365	error:
2366	sh_data_free(data);
2367	isl_basic_set_free(bset: hull);
2368	isl_set_free(set);
2369	return NULL;
2370	}
2371
2372	/* Compute a superset of the convex hull of map that is described
2373	* by only (translates of) the constraints in the constituents of map.
2374	* Handle trivial cases where map is NULL or contains at most one disjunct.
2375	*/
2376	static __isl_give isl_basic_map *map_simple_hull_trivial(
2377	__isl_take isl_map *map)
2378	{
2379	isl_basic_map *hull;
2380
2381	if (!map)
2382	return NULL;
2383	if (map->n == 0)
2384	return replace_map_by_empty_basic_map(map);
2385
2386	hull = isl_basic_map_copy(bmap: map->p[0]);
2387	isl_map_free(map);
2388	return hull;
2389	}
2390
2391	/* Return a copy of the simple hull cached inside "map".
2392	* "shift" determines whether to return the cached unshifted or shifted
2393	* simple hull.
2394	*/
2395	static __isl_give isl_basic_map *cached_simple_hull(__isl_take isl_map *map,
2396	int shift)
2397	{
2398	isl_basic_map *hull;
2399
2400	hull = isl_basic_map_copy(bmap: map->cached_simple_hull[shift]);
2401	isl_map_free(map);
2402
2403	return hull;
2404	}
2405
2406	/* Compute a superset of the convex hull of map that is described
2407	* by only (translates of) the constraints in the constituents of map.
2408	* Translation is only allowed if "shift" is set.
2409	*
2410	* The constraints are sorted while removing redundant constraints
2411	* in order to indicate a preference of which constraints should
2412	* be preserved. In particular, pairs of constraints that are
2413	* sorted together are preferred to either both be preserved
2414	* or both be removed. The sorting is performed inside
2415	* isl_basic_map_remove_redundancies.
2416	*
2417	* The result of the computation is stored in map->cached_simple_hull[shift]
2418	* such that it can be reused in subsequent calls. The cache is cleared
2419	* whenever the map is modified (in isl_map_cow).
2420	* Note that the results need to be stored in the input map for there
2421	* to be any chance that they may get reused. In particular, they
2422	* are stored in a copy of the input map that is saved before
2423	* the integer division alignment.
2424	*/
2425	static __isl_give isl_basic_map *map_simple_hull(__isl_take isl_map *map,
2426	int shift)
2427	{
2428	struct isl_set *set = NULL;
2429	struct isl_basic_map *model = NULL;
2430	struct isl_basic_map *hull;
2431	struct isl_basic_map *affine_hull;
2432	struct isl_basic_set *bset = NULL;
2433	isl_map *input;
2434
2435	if (!map || map->n <= 1)
2436	return map_simple_hull_trivial(map);
2437
2438	if (map->cached_simple_hull[shift])
2439	return cached_simple_hull(map, shift);
2440
2441	map = isl_map_detect_equalities(map);
2442	if (!map || map->n <= 1)
2443	return map_simple_hull_trivial(map);
2444	affine_hull = isl_map_affine_hull(map: isl_map_copy(map));
2445	input = isl_map_copy(map);
2446	map = isl_map_align_divs_internal(map);
2447	model = map ? isl_basic_map_copy(bmap: map->p[0]) : NULL;
2448
2449	set = isl_map_underlying_set(map);
2450
2451	bset = uset_simple_hull(set, shift);
2452
2453	hull = isl_basic_map_overlying_set(bset, like: model);
2454
2455	hull = isl_basic_map_intersect(bmap1: hull, bmap2: affine_hull);
2456	hull = isl_basic_map_remove_redundancies(bmap: hull);
2457
2458	if (hull) {
2459	ISL_F_SET(hull, ISL_BASIC_MAP_NO_IMPLICIT);
2460	ISL_F_SET(hull, ISL_BASIC_MAP_ALL_EQUALITIES);
2461	}
2462
2463	hull = isl_basic_map_finalize(bmap: hull);
2464	if (input)
2465	input->cached_simple_hull[shift] = isl_basic_map_copy(bmap: hull);
2466	isl_map_free(map: input);
2467
2468	return hull;
2469	}
2470
2471	/* Compute a superset of the convex hull of map that is described
2472	* by only translates of the constraints in the constituents of map.
2473	*/
2474	__isl_give isl_basic_map *isl_map_simple_hull(__isl_take isl_map *map)
2475	{
2476	return map_simple_hull(map, shift: 1);
2477	}
2478
2479	__isl_give isl_basic_set *isl_set_simple_hull(__isl_take isl_set *set)
2480	{
2481	return bset_from_bmap(bmap: isl_map_simple_hull(map: set_to_map(set)));
2482	}
2483
2484	/* Compute a superset of the convex hull of map that is described
2485	* by only the constraints in the constituents of map.
2486	*/
2487	__isl_give isl_basic_map *isl_map_unshifted_simple_hull(
2488	__isl_take isl_map *map)
2489	{
2490	return map_simple_hull(map, shift: 0);
2491	}
2492
2493	__isl_give isl_basic_set *isl_set_unshifted_simple_hull(
2494	__isl_take isl_set *set)
2495	{
2496	return isl_map_unshifted_simple_hull(map: set);
2497	}
2498
2499	/* Drop all inequalities from "bmap1" that do not also appear in "bmap2".
2500	* A constraint that appears with different constant terms
2501	* in "bmap1" and "bmap2" is also kept, with the least restrictive
2502	* (i.e., greatest) constant term.
2503	* "bmap1" and "bmap2" are assumed to have the same (known)
2504	* integer divisions.
2505	* The constraints of both "bmap1" and "bmap2" are assumed
2506	* to have been sorted using isl_basic_map_sort_constraints.
2507	*
2508	* Run through the inequality constraints of "bmap1" and "bmap2"
2509	* in sorted order.
2510	* Each constraint of "bmap1" without a matching constraint in "bmap2"
2511	* is removed.
2512	* If a match is found, the constraint is kept. If needed, the constant
2513	* term of the constraint is adjusted.
2514	*/
2515	static __isl_give isl_basic_map *select_shared_inequalities(
2516	__isl_take isl_basic_map *bmap1, __isl_keep isl_basic_map *bmap2)
2517	{
2518	int i1, i2;
2519
2520	bmap1 = isl_basic_map_cow(bmap: bmap1);
2521	if (!bmap1 || !bmap2)
2522	return isl_basic_map_free(bmap: bmap1);
2523
2524	i1 = bmap1->n_ineq - 1;
2525	i2 = bmap2->n_ineq - 1;
2526	while (bmap1 && i1 >= 0 && i2 >= 0) {
2527	int cmp;
2528
2529	cmp = isl_basic_map_constraint_cmp(bmap: bmap1, c1: bmap1->ineq[i1],
2530	c2: bmap2->ineq[i2]);
2531	if (cmp < 0) {
2532	--i2;
2533	continue;
2534	}
2535	if (cmp > 0) {
2536	if (isl_basic_map_drop_inequality(bmap: bmap1, pos: i1) < 0)
2537	bmap1 = isl_basic_map_free(bmap: bmap1);
2538	--i1;
2539	continue;
2540	}
2541	if (isl_int_lt(bmap1->ineq[i1][0], bmap2->ineq[i2][0]))
2542	isl_int_set(bmap1->ineq[i1][0], bmap2->ineq[i2][0]);
2543	--i1;
2544	--i2;
2545	}
2546	for (; i1 >= 0; --i1)
2547	if (isl_basic_map_drop_inequality(bmap: bmap1, pos: i1) < 0)
2548	bmap1 = isl_basic_map_free(bmap: bmap1);
2549
2550	return bmap1;
2551	}
2552
2553	/* Drop all equalities from "bmap1" that do not also appear in "bmap2".
2554	* "bmap1" and "bmap2" are assumed to have the same (known)
2555	* integer divisions.
2556	*
2557	* Run through the equality constraints of "bmap1" and "bmap2".
2558	* Each constraint of "bmap1" without a matching constraint in "bmap2"
2559	* is removed.
2560	*/
2561	static __isl_give isl_basic_map *select_shared_equalities(
2562	__isl_take isl_basic_map *bmap1, __isl_keep isl_basic_map *bmap2)
2563	{
2564	int i1, i2;
2565	isl_size total;
2566
2567	bmap1 = isl_basic_map_cow(bmap: bmap1);
2568	total = isl_basic_map_dim(bmap: bmap1, type: isl_dim_all);
2569	if (total < 0 || !bmap2)
2570	return isl_basic_map_free(bmap: bmap1);
2571
2572	i1 = bmap1->n_eq - 1;
2573	i2 = bmap2->n_eq - 1;
2574	while (bmap1 && i1 >= 0 && i2 >= 0) {
2575	int last1, last2;
2576
2577	last1 = isl_seq_last_non_zero(p: bmap1->eq[i1] + 1, len: total);
2578	last2 = isl_seq_last_non_zero(p: bmap2->eq[i2] + 1, len: total);
2579	if (last1 > last2) {
2580	--i2;
2581	continue;
2582	}
2583	if (last1 < last2) {
2584	if (isl_basic_map_drop_equality(bmap: bmap1, pos: i1) < 0)
2585	bmap1 = isl_basic_map_free(bmap: bmap1);
2586	--i1;
2587	continue;
2588	}
2589	if (!isl_seq_eq(p1: bmap1->eq[i1], p2: bmap2->eq[i2], len: 1 + total)) {
2590	if (isl_basic_map_drop_equality(bmap: bmap1, pos: i1) < 0)
2591	bmap1 = isl_basic_map_free(bmap: bmap1);
2592	}
2593	--i1;
2594	--i2;
2595	}
2596	for (; i1 >= 0; --i1)
2597	if (isl_basic_map_drop_equality(bmap: bmap1, pos: i1) < 0)
2598	bmap1 = isl_basic_map_free(bmap: bmap1);
2599
2600	return bmap1;
2601	}
2602
2603	/* Compute a superset of "bmap1" and "bmap2" that is described
2604	* by only the constraints that appear in both "bmap1" and "bmap2".
2605	*
2606	* First drop constraints that involve unknown integer divisions
2607	* since it is not trivial to check whether two such integer divisions
2608	* in different basic maps are the same.
2609	* Then align the remaining (known) divs and sort the constraints.
2610	* Finally drop all inequalities and equalities from "bmap1" that
2611	* do not also appear in "bmap2".
2612	*/
2613	__isl_give isl_basic_map *isl_basic_map_plain_unshifted_simple_hull(
2614	__isl_take isl_basic_map *bmap1, __isl_take isl_basic_map *bmap2)
2615	{
2616	if (isl_basic_map_check_equal_space(bmap1, bmap2) < 0)
2617	goto error;
2618
2619	bmap1 = isl_basic_map_drop_constraints_involving_unknown_divs(bmap: bmap1);
2620	bmap2 = isl_basic_map_drop_constraints_involving_unknown_divs(bmap: bmap2);
2621	bmap1 = isl_basic_map_order_divs(bmap: bmap1);
2622	bmap2 = isl_basic_map_align_divs(dst: bmap2, src: bmap1);
2623	bmap1 = isl_basic_map_align_divs(dst: bmap1, src: bmap2);
2624	bmap1 = isl_basic_map_sort_constraints(bmap: bmap1);
2625	bmap2 = isl_basic_map_sort_constraints(bmap: bmap2);
2626
2627	bmap1 = select_shared_inequalities(bmap1, bmap2);
2628	bmap1 = select_shared_equalities(bmap1, bmap2);
2629
2630	isl_basic_map_free(bmap: bmap2);
2631	bmap1 = isl_basic_map_finalize(bmap: bmap1);
2632	return bmap1;
2633	error:
2634	isl_basic_map_free(bmap: bmap1);
2635	isl_basic_map_free(bmap: bmap2);
2636	return NULL;
2637	}
2638
2639	/* Compute a superset of the convex hull of "map" that is described
2640	* by only the constraints in the constituents of "map".
2641	* In particular, the result is composed of constraints that appear
2642	* in each of the basic maps of "map"
2643	*
2644	* Constraints that involve unknown integer divisions are dropped
2645	* since it is not trivial to check whether two such integer divisions
2646	* in different basic maps are the same.
2647	*
2648	* The hull is initialized from the first basic map and then
2649	* updated with respect to the other basic maps in turn.
2650	*/
2651	__isl_give isl_basic_map *isl_map_plain_unshifted_simple_hull(
2652	__isl_take isl_map *map)
2653	{
2654	int i;
2655	isl_basic_map *hull;
2656
2657	if (!map)
2658	return NULL;
2659	if (map->n <= 1)
2660	return map_simple_hull_trivial(map);
2661	map = isl_map_drop_constraints_involving_unknown_divs(map);
2662	hull = isl_basic_map_copy(bmap: map->p[0]);
2663	for (i = 1; i < map->n; ++i) {
2664	isl_basic_map *bmap_i;
2665
2666	bmap_i = isl_basic_map_copy(bmap: map->p[i]);
2667	hull = isl_basic_map_plain_unshifted_simple_hull(bmap1: hull, bmap2: bmap_i);
2668	}
2669
2670	isl_map_free(map);
2671	return hull;
2672	}
2673
2674	/* Compute a superset of the convex hull of "set" that is described
2675	* by only the constraints in the constituents of "set".
2676	* In particular, the result is composed of constraints that appear
2677	* in each of the basic sets of "set"
2678	*/
2679	__isl_give isl_basic_set *isl_set_plain_unshifted_simple_hull(
2680	__isl_take isl_set *set)
2681	{
2682	return isl_map_plain_unshifted_simple_hull(map: set);
2683	}
2684
2685	/* Check if "ineq" is a bound on "set" and, if so, add it to "hull".
2686	*
2687	* For each basic set in "set", we first check if the basic set
2688	* contains a translate of "ineq". If this translate is more relaxed,
2689	* then we assume that "ineq" is not a bound on this basic set.
2690	* Otherwise, we know that it is a bound.
2691	* If the basic set does not contain a translate of "ineq", then
2692	* we call is_bound to perform the test.
2693	*/
2694	static __isl_give isl_basic_set *add_bound_from_constraint(
2695	__isl_take isl_basic_set *hull, struct sh_data *data,
2696	__isl_keep isl_set *set, isl_int *ineq)
2697	{
2698	int i, k;
2699	isl_ctx *ctx;
2700	uint32_t c_hash;
2701	struct ineq_cmp_data v;
2702	isl_size total;
2703
2704	total = isl_basic_set_dim(bset: hull, type: isl_dim_all);
2705	if (total < 0 || !set)
2706	return isl_basic_set_free(bset: hull);
2707
2708	v.len = total;
2709	v.p = ineq;
2710	c_hash = isl_seq_get_hash(p: ineq + 1, len: v.len);
2711
2712	ctx = isl_basic_set_get_ctx(bset: hull);
2713	for (i = 0; i < set->n; ++i) {
2714	int bound;
2715	struct isl_hash_table_entry *entry;
2716
2717	entry = isl_hash_table_find(ctx, table: data->p[i].table,
2718	key_hash: c_hash, eq: &has_ineq, val: &v, reserve: 0);
2719	if (!entry)
2720	return isl_basic_set_free(bset: hull);
2721	if (entry != isl_hash_table_entry_none) {
2722	isl_int *ineq_i = entry->data;
2723	int neg, more_relaxed;
2724
2725	neg = isl_seq_is_neg(p1: ineq_i + 1, p2: ineq + 1, len: v.len);
2726	if (neg)
2727	isl_int_neg(ineq_i[0], ineq_i[0]);
2728	more_relaxed = isl_int_gt(ineq_i[0], ineq[0]);
2729	if (neg)
2730	isl_int_neg(ineq_i[0], ineq_i[0]);
2731	if (more_relaxed)
2732	break;
2733	else
2734	continue;
2735	}
2736	bound = is_bound(data, set, j: i, ineq, shift: 0);
2737	if (bound < 0)
2738	return isl_basic_set_free(bset: hull);
2739	if (!bound)
2740	break;
2741	}
2742	if (i < set->n)
2743	return hull;
2744
2745	k = isl_basic_set_alloc_inequality(bset: hull);
2746	if (k < 0)
2747	return isl_basic_set_free(bset: hull);
2748	isl_seq_cpy(dst: hull->ineq[k], src: ineq, len: 1 + v.len);
2749
2750	return hull;
2751	}
2752
2753	/* Compute a superset of the convex hull of "set" that is described
2754	* by only some of the "n_ineq" constraints in the list "ineq", where "set"
2755	* has no parameters or integer divisions.
2756	*
2757	* The inequalities in "ineq" are assumed to have been sorted such
2758	* that constraints with the same linear part appear together and
2759	* that among constraints with the same linear part, those with
2760	* smaller constant term appear first.
2761	*
2762	* We reuse the same data structure that is used by uset_simple_hull,
2763	* but we do not need the hull table since we will not consider the
2764	* same constraint more than once. We therefore allocate it with zero size.
2765	*
2766	* We run through the constraints and try to add them one by one,
2767	* skipping identical constraints. If we have added a constraint and
2768	* the next constraint is a more relaxed translate, then we skip this
2769	* next constraint as well.
2770	*/
2771	static __isl_give isl_basic_set *uset_unshifted_simple_hull_from_constraints(
2772	__isl_take isl_set *set, int n_ineq, isl_int **ineq)
2773	{
2774	int i;
2775	int last_added = 0;
2776	struct sh_data *data = NULL;
2777	isl_basic_set *hull = NULL;
2778	isl_size dim;
2779
2780	hull = isl_basic_set_alloc_space(space: isl_set_get_space(set), extra: 0, n_eq: 0, n_ineq);
2781	if (!hull)
2782	goto error;
2783
2784	data = sh_data_alloc(set, n_ineq: 0);
2785	if (!data)
2786	goto error;
2787
2788	dim = isl_set_dim(set, type: isl_dim_set);
2789	if (dim < 0)
2790	goto error;
2791	for (i = 0; i < n_ineq; ++i) {
2792	int hull_n_ineq = hull->n_ineq;
2793	int parallel;
2794
2795	parallel = i > 0 && isl_seq_eq(p1: ineq[i - 1] + 1, p2: ineq[i] + 1,
2796	len: dim);
2797	if (parallel &&
2798	(last_added || isl_int_eq(ineq[i - 1][0], ineq[i][0])))
2799	continue;
2800	hull = add_bound_from_constraint(hull, data, set, ineq: ineq[i]);
2801	if (!hull)
2802	goto error;
2803	last_added = hull->n_ineq > hull_n_ineq;
2804	}
2805
2806	sh_data_free(data);
2807	isl_set_free(set);
2808	return hull;
2809	error:
2810	sh_data_free(data);
2811	isl_set_free(set);
2812	isl_basic_set_free(bset: hull);
2813	return NULL;
2814	}
2815
2816	/* Collect pointers to all the inequalities in the elements of "list"
2817	* in "ineq". For equalities, store both a pointer to the equality and
2818	* a pointer to its opposite, which is first copied to "mat".
2819	* "ineq" and "mat" are assumed to have been preallocated to the right size
2820	* (the number of inequalities + 2 times the number of equalites and
2821	* the number of equalities, respectively).
2822	*/
2823	static __isl_give isl_mat *collect_inequalities(__isl_take isl_mat *mat,
2824	__isl_keep isl_basic_set_list *list, isl_int **ineq)
2825	{
2826	int i, j, n_eq, n_ineq;
2827	isl_size n;
2828
2829	n = isl_basic_set_list_n_basic_set(list);
2830	if (!mat || n < 0)
2831	return isl_mat_free(mat);
2832
2833	n_eq = 0;
2834	n_ineq = 0;
2835	for (i = 0; i < n; ++i) {
2836	isl_basic_set *bset;
2837	bset = isl_basic_set_list_get_basic_set(list, index: i);
2838	if (!bset)
2839	return isl_mat_free(mat);
2840	for (j = 0; j < bset->n_eq; ++j) {
2841	ineq[n_ineq++] = mat->row[n_eq];
2842	ineq[n_ineq++] = bset->eq[j];
2843	isl_seq_neg(dst: mat->row[n_eq++], src: bset->eq[j], len: mat->n_col);
2844	}
2845	for (j = 0; j < bset->n_ineq; ++j)
2846	ineq[n_ineq++] = bset->ineq[j];
2847	isl_basic_set_free(bset);
2848	}
2849
2850	return mat;
2851	}
2852
2853	/* Comparison routine for use as an isl_sort callback.
2854	*
2855	* Constraints with the same linear part are sorted together and
2856	* among constraints with the same linear part, those with smaller
2857	* constant term are sorted first.
2858	*/
2859	static int cmp_ineq(const void *a, const void *b, void *arg)
2860	{
2861	unsigned dim = *(unsigned *) arg;
2862	isl_int * const *ineq1 = a;
2863	isl_int * const *ineq2 = b;
2864	int cmp;
2865
2866	cmp = isl_seq_cmp(p1: (*ineq1) + 1, p2: (*ineq2) + 1, len: dim);
2867	if (cmp != 0)
2868	return cmp;
2869	return isl_int_cmp((*ineq1)[0], (*ineq2)[0]);
2870	}
2871
2872	/* Compute a superset of the convex hull of "set" that is described
2873	* by only constraints in the elements of "list", where "set" has
2874	* no parameters or integer divisions.
2875	*
2876	* We collect all the constraints in those elements and then
2877	* sort the constraints such that constraints with the same linear part
2878	* are sorted together and that those with smaller constant term are
2879	* sorted first.
2880	*/
2881	static __isl_give isl_basic_set *uset_unshifted_simple_hull_from_basic_set_list(
2882	__isl_take isl_set *set, __isl_take isl_basic_set_list *list)
2883	{
2884	int i, n_eq, n_ineq;
2885	isl_size n;
2886	isl_size dim;
2887	isl_ctx *ctx;
2888	isl_mat *mat = NULL;
2889	isl_int **ineq = NULL;
2890	isl_basic_set *hull;
2891
2892	n = isl_basic_set_list_n_basic_set(list);
2893	if (!set || n < 0)
2894	goto error;
2895	ctx = isl_set_get_ctx(set);
2896
2897	n_eq = 0;
2898	n_ineq = 0;
2899	for (i = 0; i < n; ++i) {
2900	isl_basic_set *bset;
2901	bset = isl_basic_set_list_get_basic_set(list, index: i);
2902	if (!bset)
2903	goto error;
2904	n_eq += bset->n_eq;
2905	n_ineq += 2 * bset->n_eq + bset->n_ineq;
2906	isl_basic_set_free(bset);
2907	}
2908
2909	ineq = isl_alloc_array(ctx, isl_int *, n_ineq);
2910	if (n_ineq > 0 && !ineq)
2911	goto error;
2912
2913	dim = isl_set_dim(set, type: isl_dim_set);
2914	if (dim < 0)
2915	goto error;
2916	mat = isl_mat_alloc(ctx, n_row: n_eq, n_col: 1 + dim);
2917	mat = collect_inequalities(mat, list, ineq);
2918	if (!mat)
2919	goto error;
2920
2921	if (isl_sort(pbase: ineq, total_elems: n_ineq, size: sizeof(ineq[0]), cmp: &cmp_ineq, arg: &dim) < 0)
2922	goto error;
2923
2924	hull = uset_unshifted_simple_hull_from_constraints(set, n_ineq, ineq);
2925
2926	isl_mat_free(mat);
2927	free(ptr: ineq);
2928	isl_basic_set_list_free(list);
2929	return hull;
2930	error:
2931	isl_mat_free(mat);
2932	free(ptr: ineq);
2933	isl_set_free(set);
2934	isl_basic_set_list_free(list);
2935	return NULL;
2936	}
2937
2938	/* Compute a superset of the convex hull of "map" that is described
2939	* by only constraints in the elements of "list".
2940	*
2941	* If the list is empty, then we can only describe the universe set.
2942	* If the input map is empty, then all constraints are valid, so
2943	* we return the intersection of the elements in "list".
2944	*
2945	* Otherwise, we align all divs and temporarily treat them
2946	* as regular variables, computing the unshifted simple hull in
2947	* uset_unshifted_simple_hull_from_basic_set_list.
2948	*/
2949	static __isl_give isl_basic_map *map_unshifted_simple_hull_from_basic_map_list(
2950	__isl_take isl_map *map, __isl_take isl_basic_map_list *list)
2951	{
2952	isl_size n;
2953	isl_basic_map *model;
2954	isl_basic_map *hull;
2955	isl_set *set;
2956	isl_basic_set_list *bset_list;
2957
2958	n = isl_basic_map_list_n_basic_map(list);
2959	if (!map || n < 0)
2960	goto error;
2961
2962	if (n == 0) {
2963	isl_space *space;
2964
2965	space = isl_map_get_space(map);
2966	isl_map_free(map);
2967	isl_basic_map_list_free(list);
2968	return isl_basic_map_universe(space);
2969	}
2970	if (isl_map_plain_is_empty(map)) {
2971	isl_map_free(map);
2972	return isl_basic_map_list_intersect(list);
2973	}
2974
2975	map = isl_map_align_divs_to_basic_map_list(map, list);
2976	if (!map)
2977	goto error;
2978	list = isl_basic_map_list_align_divs_to_basic_map(list, bmap: map->p[0]);
2979
2980	model = isl_basic_map_list_get_basic_map(list, index: 0);
2981
2982	set = isl_map_underlying_set(map);
2983	bset_list = isl_basic_map_list_underlying_set(list);
2984
2985	hull = uset_unshifted_simple_hull_from_basic_set_list(set, list: bset_list);
2986	hull = isl_basic_map_overlying_set(bset: hull, like: model);
2987
2988	return hull;
2989	error:
2990	isl_map_free(map);
2991	isl_basic_map_list_free(list);
2992	return NULL;
2993	}
2994
2995	/* Return a sequence of the basic maps that make up the maps in "list".
2996	*/
2997	static __isl_give isl_basic_map_list *collect_basic_maps(
2998	__isl_take isl_map_list *list)
2999	{
3000	int i;
3001	isl_size n;
3002	isl_ctx *ctx;
3003	isl_basic_map_list *bmap_list;
3004
3005	if (!list)
3006	return NULL;
3007	n = isl_map_list_n_map(list);
3008	ctx = isl_map_list_get_ctx(list);
3009	bmap_list = isl_basic_map_list_alloc(ctx, n: 0);
3010	if (n < 0)
3011	bmap_list = isl_basic_map_list_free(list: bmap_list);
3012
3013	for (i = 0; i < n; ++i) {
3014	isl_map *map;
3015	isl_basic_map_list *list_i;
3016
3017	map = isl_map_list_get_map(list, index: i);
3018	map = isl_map_compute_divs(map);
3019	list_i = isl_map_get_basic_map_list(map);
3020	isl_map_free(map);
3021	bmap_list = isl_basic_map_list_concat(list1: bmap_list, list2: list_i);
3022	}
3023
3024	isl_map_list_free(list);
3025	return bmap_list;
3026	}
3027
3028	/* Compute a superset of the convex hull of "map" that is described
3029	* by only constraints in the elements of "list".
3030	*
3031	* If "map" is the universe, then the convex hull (and therefore
3032	* any superset of the convexhull) is the universe as well.
3033	*
3034	* Otherwise, we collect all the basic maps in the map list and
3035	* continue with map_unshifted_simple_hull_from_basic_map_list.
3036	*/
3037	__isl_give isl_basic_map *isl_map_unshifted_simple_hull_from_map_list(
3038	__isl_take isl_map *map, __isl_take isl_map_list *list)
3039	{
3040	isl_basic_map_list *bmap_list;
3041	int is_universe;
3042
3043	is_universe = isl_map_plain_is_universe(map);
3044	if (is_universe < 0)
3045	map = isl_map_free(map);
3046	if (is_universe < 0 || is_universe) {
3047	isl_map_list_free(list);
3048	return isl_map_unshifted_simple_hull(map);
3049	}
3050
3051	bmap_list = collect_basic_maps(list);
3052	return map_unshifted_simple_hull_from_basic_map_list(map, list: bmap_list);
3053	}
3054
3055	/* Compute a superset of the convex hull of "set" that is described
3056	* by only constraints in the elements of "list".
3057	*/
3058	__isl_give isl_basic_set *isl_set_unshifted_simple_hull_from_set_list(
3059	__isl_take isl_set *set, __isl_take isl_set_list *list)
3060	{
3061	return isl_map_unshifted_simple_hull_from_map_list(map: set, list);
3062	}
3063
3064	/* Given a set "set", return parametric bounds on the dimension "dim".
3065	*/
3066	static __isl_give isl_basic_set *set_bounds(__isl_keep isl_set *set, int dim)
3067	{
3068	isl_size set_dim = isl_set_dim(set, type: isl_dim_set);
3069	if (set_dim < 0)
3070	return NULL;
3071	set = isl_set_copy(set);
3072	set = isl_set_eliminate_dims(set, first: dim + 1, n: set_dim - (dim + 1));
3073	set = isl_set_eliminate_dims(set, first: 0, n: dim);
3074	return isl_set_convex_hull(set);
3075	}
3076
3077	/* Computes a "simple hull" and then check if each dimension in the
3078	* resulting hull is bounded by a symbolic constant. If not, the
3079	* hull is intersected with the corresponding bounds on the whole set.
3080	*/
3081	__isl_give isl_basic_set *isl_set_bounded_simple_hull(__isl_take isl_set *set)
3082	{
3083	int i, j;
3084	struct isl_basic_set *hull;
3085	isl_size nparam, dim, total;
3086	unsigned left;
3087	int removed_divs = 0;
3088
3089	hull = isl_set_simple_hull(set: isl_set_copy(set));
3090	nparam = isl_basic_set_dim(bset: hull, type: isl_dim_param);
3091	dim = isl_basic_set_dim(bset: hull, type: isl_dim_set);
3092	total = isl_basic_set_dim(bset: hull, type: isl_dim_all);
3093	if (nparam < 0 || dim < 0 || total < 0)
3094	goto error;
3095
3096	for (i = 0; i < dim; ++i) {
3097	int lower = 0, upper = 0;
3098	struct isl_basic_set *bounds;
3099
3100	left = total - nparam - i - 1;
3101	for (j = 0; j < hull->n_eq; ++j) {
3102	if (isl_int_is_zero(hull->eq[j][1 + nparam + i]))
3103	continue;
3104	if (isl_seq_first_non_zero(p: hull->eq[j]+1+nparam+i+1,
3105	len: left) == -1)
3106	break;
3107	}
3108	if (j < hull->n_eq)
3109	continue;
3110
3111	for (j = 0; j < hull->n_ineq; ++j) {
3112	if (isl_int_is_zero(hull->ineq[j][1 + nparam + i]))
3113	continue;
3114	if (isl_seq_first_non_zero(p: hull->ineq[j]+1+nparam+i+1,
3115	len: left) != -1 ||
3116	isl_seq_first_non_zero(p: hull->ineq[j]+1+nparam,
3117	len: i) != -1)
3118	continue;
3119	if (isl_int_is_pos(hull->ineq[j][1 + nparam + i]))
3120	lower = 1;
3121	else
3122	upper = 1;
3123	if (lower && upper)
3124	break;
3125	}
3126
3127	if (lower && upper)
3128	continue;
3129
3130	if (!removed_divs) {
3131	set = isl_set_remove_divs(set);
3132	if (!set)
3133	goto error;
3134	removed_divs = 1;
3135	}
3136	bounds = set_bounds(set, dim: i);
3137	hull = isl_basic_set_intersect(bset1: hull, bset2: bounds);
3138	if (!hull)
3139	goto error;
3140	}
3141
3142	isl_set_free(set);
3143	return hull;
3144	error:
3145	isl_set_free(set);
3146	isl_basic_set_free(bset: hull);
3147	return NULL;
3148	}
3149