isl_constraint.c: use isl_basic_map_offset
[isl.git] / isl_convex_hull.c
blobaf10145ad2d0a3867cc13f7a2e0c19396bd15356
1 /*
2 * Copyright 2008-2009 Katholieke Universiteit Leuven
3 * Copyright 2014 INRIA Rocquencourt
5 * Use of this software is governed by the MIT license
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
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>
26 #include <bset_to_bmap.c>
27 #include <bset_from_bmap.c>
28 #include <set_to_map.c>
30 static __isl_give isl_basic_set *uset_convex_hull_wrap_bounded(
31 __isl_take isl_set *set);
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.
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.
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.
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.
51 __isl_give isl_basic_map *isl_basic_map_remove_redundancies(
52 __isl_take isl_basic_map *bmap)
54 struct isl_tab *tab;
56 if (!bmap)
57 return NULL;
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;
67 bmap = isl_basic_map_sort_constraints(bmap);
68 tab = isl_tab_from_basic_map(bmap, 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;
92 __isl_give isl_basic_set *isl_basic_set_remove_redundancies(
93 __isl_take isl_basic_set *bset)
95 return bset_from_bmap(
96 isl_basic_map_remove_redundancies(bset_to_bmap(bset)));
99 /* Remove redundant constraints in each of the basic maps.
101 __isl_give isl_map *isl_map_remove_redundancies(__isl_take isl_map *map)
103 return isl_map_inline_foreach_basic_map(map,
104 &isl_basic_map_remove_redundancies);
107 __isl_give isl_set *isl_set_remove_redundancies(__isl_take isl_set *set)
109 return isl_map_remove_redundancies(set);
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.
116 static int uset_is_bound(__isl_keep isl_set *set, isl_int *c, unsigned len)
118 int first;
119 int j;
120 isl_int opt;
121 isl_int opt_denom;
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;
129 if (ISL_F_ISSET(set->p[j], ISL_BASIC_SET_EMPTY))
130 continue;
132 res = isl_basic_set_solve_lp(set->p[j],
133 0, c, set->ctx->one, &opt, &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(set->p[j]);
140 if (!set->p[j])
141 goto error;
142 continue;
144 if (first || isl_int_is_neg(opt)) {
145 if (!isl_int_is_one(opt_denom))
146 isl_seq_scale(c, c, opt_denom, len);
147 isl_int_sub(c[0], c[0], opt);
149 first = 0;
151 isl_int_clear(opt);
152 isl_int_clear(opt_denom);
153 return j >= set->n;
154 error:
155 isl_int_clear(opt);
156 isl_int_clear(opt_denom);
157 return -1;
160 static __isl_give isl_set *isl_set_add_basic_set_equality(
161 __isl_take isl_set *set, isl_int *c)
163 int i;
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(set->p[i], c);
170 if (!set->p[i])
171 goto error;
173 return set;
174 error:
175 isl_set_free(set);
176 return NULL;
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
187 * then the resulting set is of dimension n*(1+d) and has as constraints
189 * [ a_i ]
190 * A_i [ x_i ] >= 0
192 * a_i >= 0
194 * \sum_i x_{i,1} = 1
196 static __isl_give isl_basic_set *wrap_constraints(__isl_keep isl_set *set)
198 struct isl_basic_set *lp;
199 unsigned n_eq;
200 unsigned n_ineq;
201 int i, j, k;
202 unsigned dim, lp_dim;
204 if (!set)
205 return NULL;
207 dim = 1 + isl_set_n_dim(set);
208 n_eq = 1;
209 n_ineq = set->n;
210 for (i = 0; i < set->n; ++i) {
211 n_eq += set->p[i]->n_eq;
212 n_ineq += set->p[i]->n_ineq;
214 lp = isl_basic_set_alloc(set->ctx, 0, dim * set->n, 0, n_eq, n_ineq);
215 lp = isl_basic_set_set_rational(lp);
216 if (!lp)
217 return NULL;
218 lp_dim = isl_basic_set_n_dim(lp);
219 k = isl_basic_set_alloc_equality(lp);
220 isl_int_set_si(lp->eq[k][0], -1);
221 for (i = 0; i < set->n; ++i) {
222 isl_int_set_si(lp->eq[k][1+dim*i], 0);
223 isl_int_set_si(lp->eq[k][1+dim*i+1], 1);
224 isl_seq_clr(lp->eq[k]+1+dim*i+2, dim-2);
226 for (i = 0; i < set->n; ++i) {
227 k = isl_basic_set_alloc_inequality(lp);
228 isl_seq_clr(lp->ineq[k], 1+lp_dim);
229 isl_int_set_si(lp->ineq[k][1+dim*i], 1);
231 for (j = 0; j < set->p[i]->n_eq; ++j) {
232 k = isl_basic_set_alloc_equality(lp);
233 isl_seq_clr(lp->eq[k], 1+dim*i);
234 isl_seq_cpy(lp->eq[k]+1+dim*i, set->p[i]->eq[j], dim);
235 isl_seq_clr(lp->eq[k]+1+dim*(i+1), dim*(set->n-i-1));
238 for (j = 0; j < set->p[i]->n_ineq; ++j) {
239 k = isl_basic_set_alloc_inequality(lp);
240 isl_seq_clr(lp->ineq[k], 1+dim*i);
241 isl_seq_cpy(lp->ineq[k]+1+dim*i, set->p[i]->ineq[j], dim);
242 isl_seq_clr(lp->ineq[k]+1+dim*(i+1), dim*(set->n-i-1));
245 return lp;
248 /* Given a facet "facet" of the convex hull of "set" and a facet "ridge"
249 * of that facet, compute the other facet of the convex hull that contains
250 * the ridge.
252 * We first transform the set such that the facet constraint becomes
254 * x_1 >= 0
256 * I.e., the facet lies in
258 * x_1 = 0
260 * and on that facet, the constraint that defines the ridge is
262 * x_2 >= 0
264 * (This transformation is not strictly needed, all that is needed is
265 * that the ridge contains the origin.)
267 * Since the ridge contains the origin, the cone of the convex hull
268 * will be of the form
270 * x_1 >= 0
271 * x_2 >= a x_1
273 * with this second constraint defining the new facet.
274 * The constant a is obtained by settting x_1 in the cone of the
275 * convex hull to 1 and minimizing x_2.
276 * Now, each element in the cone of the convex hull is the sum
277 * of elements in the cones of the basic sets.
278 * If a_i is the dilation factor of basic set i, then the problem
279 * we need to solve is
281 * min \sum_i x_{i,2}
282 * st
283 * \sum_i x_{i,1} = 1
284 * a_i >= 0
285 * [ a_i ]
286 * A [ x_i ] >= 0
288 * with
289 * [ 1 ]
290 * A_i [ x_i ] >= 0
292 * the constraints of each (transformed) basic set.
293 * If a = n/d, then the constraint defining the new facet (in the transformed
294 * space) is
296 * -n x_1 + d x_2 >= 0
298 * In the original space, we need to take the same combination of the
299 * corresponding constraints "facet" and "ridge".
301 * If a = -infty = "-1/0", then we just return the original facet constraint.
302 * This means that the facet is unbounded, but has a bounded intersection
303 * with the union of sets.
305 isl_int *isl_set_wrap_facet(__isl_keep isl_set *set,
306 isl_int *facet, isl_int *ridge)
308 int i;
309 isl_ctx *ctx;
310 struct isl_mat *T = NULL;
311 struct isl_basic_set *lp = NULL;
312 struct isl_vec *obj;
313 enum isl_lp_result res;
314 isl_int num, den;
315 unsigned dim;
317 if (!set)
318 return NULL;
319 ctx = set->ctx;
320 set = isl_set_copy(set);
321 set = isl_set_set_rational(set);
323 dim = 1 + isl_set_n_dim(set);
324 T = isl_mat_alloc(ctx, 3, dim);
325 if (!T)
326 goto error;
327 isl_int_set_si(T->row[0][0], 1);
328 isl_seq_clr(T->row[0]+1, dim - 1);
329 isl_seq_cpy(T->row[1], facet, dim);
330 isl_seq_cpy(T->row[2], ridge, dim);
331 T = isl_mat_right_inverse(T);
332 set = isl_set_preimage(set, T);
333 T = NULL;
334 if (!set)
335 goto error;
336 lp = wrap_constraints(set);
337 obj = isl_vec_alloc(ctx, 1 + dim*set->n);
338 if (!obj)
339 goto error;
340 isl_int_set_si(obj->block.data[0], 0);
341 for (i = 0; i < set->n; ++i) {
342 isl_seq_clr(obj->block.data + 1 + dim*i, 2);
343 isl_int_set_si(obj->block.data[1 + dim*i+2], 1);
344 isl_seq_clr(obj->block.data + 1 + dim*i+3, dim-3);
346 isl_int_init(num);
347 isl_int_init(den);
348 res = isl_basic_set_solve_lp(lp, 0,
349 obj->block.data, ctx->one, &num, &den, NULL);
350 if (res == isl_lp_ok) {
351 isl_int_neg(num, num);
352 isl_seq_combine(facet, num, facet, den, ridge, dim);
353 isl_seq_normalize(ctx, facet, dim);
355 isl_int_clear(num);
356 isl_int_clear(den);
357 isl_vec_free(obj);
358 isl_basic_set_free(lp);
359 isl_set_free(set);
360 if (res == isl_lp_error)
361 return NULL;
362 isl_assert(ctx, res == isl_lp_ok || res == isl_lp_unbounded,
363 return NULL);
364 return facet;
365 error:
366 isl_basic_set_free(lp);
367 isl_mat_free(T);
368 isl_set_free(set);
369 return NULL;
372 /* Compute the constraint of a facet of "set".
374 * We first compute the intersection with a bounding constraint
375 * that is orthogonal to one of the coordinate axes.
376 * If the affine hull of this intersection has only one equality,
377 * we have found a facet.
378 * Otherwise, we wrap the current bounding constraint around
379 * one of the equalities of the face (one that is not equal to
380 * the current bounding constraint).
381 * This process continues until we have found a facet.
382 * The dimension of the intersection increases by at least
383 * one on each iteration, so termination is guaranteed.
385 static __isl_give isl_mat *initial_facet_constraint(__isl_keep isl_set *set)
387 struct isl_set *slice = NULL;
388 struct isl_basic_set *face = NULL;
389 int i;
390 unsigned dim = isl_set_n_dim(set);
391 int is_bound;
392 isl_mat *bounds = NULL;
394 isl_assert(set->ctx, set->n > 0, goto error);
395 bounds = isl_mat_alloc(set->ctx, 1, 1 + dim);
396 if (!bounds)
397 return NULL;
399 isl_seq_clr(bounds->row[0], dim);
400 isl_int_set_si(bounds->row[0][1 + dim - 1], 1);
401 is_bound = uset_is_bound(set, bounds->row[0], 1 + dim);
402 if (is_bound < 0)
403 goto error;
404 isl_assert(set->ctx, is_bound, goto error);
405 isl_seq_normalize(set->ctx, bounds->row[0], 1 + dim);
406 bounds->n_row = 1;
408 for (;;) {
409 slice = isl_set_copy(set);
410 slice = isl_set_add_basic_set_equality(slice, bounds->row[0]);
411 face = isl_set_affine_hull(slice);
412 if (!face)
413 goto error;
414 if (face->n_eq == 1) {
415 isl_basic_set_free(face);
416 break;
418 for (i = 0; i < face->n_eq; ++i)
419 if (!isl_seq_eq(bounds->row[0], face->eq[i], 1 + dim) &&
420 !isl_seq_is_neg(bounds->row[0],
421 face->eq[i], 1 + dim))
422 break;
423 isl_assert(set->ctx, i < face->n_eq, goto error);
424 if (!isl_set_wrap_facet(set, bounds->row[0], face->eq[i]))
425 goto error;
426 isl_seq_normalize(set->ctx, bounds->row[0], bounds->n_col);
427 isl_basic_set_free(face);
430 return bounds;
431 error:
432 isl_basic_set_free(face);
433 isl_mat_free(bounds);
434 return NULL;
437 /* Given the bounding constraint "c" of a facet of the convex hull of "set",
438 * compute a hyperplane description of the facet, i.e., compute the facets
439 * of the facet.
441 * We compute an affine transformation that transforms the constraint
443 * [ 1 ]
444 * c [ x ] = 0
446 * to the constraint
448 * z_1 = 0
450 * by computing the right inverse U of a matrix that starts with the rows
452 * [ 1 0 ]
453 * [ c ]
455 * Then
456 * [ 1 ] [ 1 ]
457 * [ x ] = U [ z ]
458 * and
459 * [ 1 ] [ 1 ]
460 * [ z ] = Q [ x ]
462 * with Q = U^{-1}
463 * Since z_1 is zero, we can drop this variable as well as the corresponding
464 * column of U to obtain
466 * [ 1 ] [ 1 ]
467 * [ x ] = U' [ z' ]
468 * and
469 * [ 1 ] [ 1 ]
470 * [ z' ] = Q' [ x ]
472 * with Q' equal to Q, but without the corresponding row.
473 * After computing the facets of the facet in the z' space,
474 * we convert them back to the x space through Q.
476 static __isl_give isl_basic_set *compute_facet(__isl_keep isl_set *set,
477 isl_int *c)
479 struct isl_mat *m, *U, *Q;
480 struct isl_basic_set *facet = NULL;
481 struct isl_ctx *ctx;
482 unsigned dim;
484 ctx = set->ctx;
485 set = isl_set_copy(set);
486 dim = isl_set_n_dim(set);
487 m = isl_mat_alloc(set->ctx, 2, 1 + dim);
488 if (!m)
489 goto error;
490 isl_int_set_si(m->row[0][0], 1);
491 isl_seq_clr(m->row[0]+1, dim);
492 isl_seq_cpy(m->row[1], c, 1+dim);
493 U = isl_mat_right_inverse(m);
494 Q = isl_mat_right_inverse(isl_mat_copy(U));
495 U = isl_mat_drop_cols(U, 1, 1);
496 Q = isl_mat_drop_rows(Q, 1, 1);
497 set = isl_set_preimage(set, U);
498 facet = uset_convex_hull_wrap_bounded(set);
499 facet = isl_basic_set_preimage(facet, Q);
500 if (facet && facet->n_eq != 0)
501 isl_die(ctx, isl_error_internal, "unexpected equality",
502 return isl_basic_set_free(facet));
503 return facet;
504 error:
505 isl_basic_set_free(facet);
506 isl_set_free(set);
507 return NULL;
510 /* Given an initial facet constraint, compute the remaining facets.
511 * We do this by running through all facets found so far and computing
512 * the adjacent facets through wrapping, adding those facets that we
513 * hadn't already found before.
515 * For each facet we have found so far, we first compute its facets
516 * in the resulting convex hull. That is, we compute the ridges
517 * of the resulting convex hull contained in the facet.
518 * We also compute the corresponding facet in the current approximation
519 * of the convex hull. There is no need to wrap around the ridges
520 * in this facet since that would result in a facet that is already
521 * present in the current approximation.
523 * This function can still be significantly optimized by checking which of
524 * the facets of the basic sets are also facets of the convex hull and
525 * using all the facets so far to help in constructing the facets of the
526 * facets
527 * and/or
528 * using the technique in section "3.1 Ridge Generation" of
529 * "Extended Convex Hull" by Fukuda et al.
531 static __isl_give isl_basic_set *extend(__isl_take isl_basic_set *hull,
532 __isl_keep isl_set *set)
534 int i, j, f;
535 int k;
536 struct isl_basic_set *facet = NULL;
537 struct isl_basic_set *hull_facet = NULL;
538 unsigned dim;
540 if (!hull)
541 return NULL;
543 isl_assert(set->ctx, set->n > 0, goto error);
545 dim = isl_set_n_dim(set);
547 for (i = 0; i < hull->n_ineq; ++i) {
548 facet = compute_facet(set, hull->ineq[i]);
549 facet = isl_basic_set_add_eq(facet, hull->ineq[i]);
550 facet = isl_basic_set_gauss(facet, NULL);
551 facet = isl_basic_set_normalize_constraints(facet);
552 hull_facet = isl_basic_set_copy(hull);
553 hull_facet = isl_basic_set_add_eq(hull_facet, hull->ineq[i]);
554 hull_facet = isl_basic_set_gauss(hull_facet, NULL);
555 hull_facet = isl_basic_set_normalize_constraints(hull_facet);
556 if (!facet || !hull_facet)
557 goto error;
558 hull = isl_basic_set_cow(hull);
559 hull = isl_basic_set_extend_space(hull,
560 isl_space_copy(hull->dim), 0, 0, facet->n_ineq);
561 if (!hull)
562 goto error;
563 for (j = 0; j < facet->n_ineq; ++j) {
564 for (f = 0; f < hull_facet->n_ineq; ++f)
565 if (isl_seq_eq(facet->ineq[j],
566 hull_facet->ineq[f], 1 + dim))
567 break;
568 if (f < hull_facet->n_ineq)
569 continue;
570 k = isl_basic_set_alloc_inequality(hull);
571 if (k < 0)
572 goto error;
573 isl_seq_cpy(hull->ineq[k], hull->ineq[i], 1+dim);
574 if (!isl_set_wrap_facet(set, hull->ineq[k], facet->ineq[j]))
575 goto error;
577 isl_basic_set_free(hull_facet);
578 isl_basic_set_free(facet);
580 hull = isl_basic_set_simplify(hull);
581 hull = isl_basic_set_finalize(hull);
582 return hull;
583 error:
584 isl_basic_set_free(hull_facet);
585 isl_basic_set_free(facet);
586 isl_basic_set_free(hull);
587 return NULL;
590 /* Special case for computing the convex hull of a one dimensional set.
591 * We simply collect the lower and upper bounds of each basic set
592 * and the biggest of those.
594 static __isl_give isl_basic_set *convex_hull_1d(__isl_take isl_set *set)
596 struct isl_mat *c = NULL;
597 isl_int *lower = NULL;
598 isl_int *upper = NULL;
599 int i, j, k;
600 isl_int a, b;
601 struct isl_basic_set *hull;
603 for (i = 0; i < set->n; ++i) {
604 set->p[i] = isl_basic_set_simplify(set->p[i]);
605 if (!set->p[i])
606 goto error;
608 set = isl_set_remove_empty_parts(set);
609 if (!set)
610 goto error;
611 isl_assert(set->ctx, set->n > 0, goto error);
612 c = isl_mat_alloc(set->ctx, 2, 2);
613 if (!c)
614 goto error;
616 if (set->p[0]->n_eq > 0) {
617 isl_assert(set->ctx, set->p[0]->n_eq == 1, goto error);
618 lower = c->row[0];
619 upper = c->row[1];
620 if (isl_int_is_pos(set->p[0]->eq[0][1])) {
621 isl_seq_cpy(lower, set->p[0]->eq[0], 2);
622 isl_seq_neg(upper, set->p[0]->eq[0], 2);
623 } else {
624 isl_seq_neg(lower, set->p[0]->eq[0], 2);
625 isl_seq_cpy(upper, set->p[0]->eq[0], 2);
627 } else {
628 for (j = 0; j < set->p[0]->n_ineq; ++j) {
629 if (isl_int_is_pos(set->p[0]->ineq[j][1])) {
630 lower = c->row[0];
631 isl_seq_cpy(lower, set->p[0]->ineq[j], 2);
632 } else {
633 upper = c->row[1];
634 isl_seq_cpy(upper, set->p[0]->ineq[j], 2);
639 isl_int_init(a);
640 isl_int_init(b);
641 for (i = 0; i < set->n; ++i) {
642 struct isl_basic_set *bset = set->p[i];
643 int has_lower = 0;
644 int has_upper = 0;
646 for (j = 0; j < bset->n_eq; ++j) {
647 has_lower = 1;
648 has_upper = 1;
649 if (lower) {
650 isl_int_mul(a, lower[0], bset->eq[j][1]);
651 isl_int_mul(b, lower[1], bset->eq[j][0]);
652 if (isl_int_lt(a, b) && isl_int_is_pos(bset->eq[j][1]))
653 isl_seq_cpy(lower, bset->eq[j], 2);
654 if (isl_int_gt(a, b) && isl_int_is_neg(bset->eq[j][1]))
655 isl_seq_neg(lower, bset->eq[j], 2);
657 if (upper) {
658 isl_int_mul(a, upper[0], bset->eq[j][1]);
659 isl_int_mul(b, upper[1], bset->eq[j][0]);
660 if (isl_int_lt(a, b) && isl_int_is_pos(bset->eq[j][1]))
661 isl_seq_neg(upper, bset->eq[j], 2);
662 if (isl_int_gt(a, b) && isl_int_is_neg(bset->eq[j][1]))
663 isl_seq_cpy(upper, bset->eq[j], 2);
666 for (j = 0; j < bset->n_ineq; ++j) {
667 if (isl_int_is_pos(bset->ineq[j][1]))
668 has_lower = 1;
669 if (isl_int_is_neg(bset->ineq[j][1]))
670 has_upper = 1;
671 if (lower && isl_int_is_pos(bset->ineq[j][1])) {
672 isl_int_mul(a, lower[0], bset->ineq[j][1]);
673 isl_int_mul(b, lower[1], bset->ineq[j][0]);
674 if (isl_int_lt(a, b))
675 isl_seq_cpy(lower, bset->ineq[j], 2);
677 if (upper && isl_int_is_neg(bset->ineq[j][1])) {
678 isl_int_mul(a, upper[0], bset->ineq[j][1]);
679 isl_int_mul(b, upper[1], bset->ineq[j][0]);
680 if (isl_int_gt(a, b))
681 isl_seq_cpy(upper, bset->ineq[j], 2);
684 if (!has_lower)
685 lower = NULL;
686 if (!has_upper)
687 upper = NULL;
689 isl_int_clear(a);
690 isl_int_clear(b);
692 hull = isl_basic_set_alloc(set->ctx, 0, 1, 0, 0, 2);
693 hull = isl_basic_set_set_rational(hull);
694 if (!hull)
695 goto error;
696 if (lower) {
697 k = isl_basic_set_alloc_inequality(hull);
698 isl_seq_cpy(hull->ineq[k], lower, 2);
700 if (upper) {
701 k = isl_basic_set_alloc_inequality(hull);
702 isl_seq_cpy(hull->ineq[k], upper, 2);
704 hull = isl_basic_set_finalize(hull);
705 isl_set_free(set);
706 isl_mat_free(c);
707 return hull;
708 error:
709 isl_set_free(set);
710 isl_mat_free(c);
711 return NULL;
714 static __isl_give isl_basic_set *convex_hull_0d(__isl_take isl_set *set)
716 struct isl_basic_set *convex_hull;
718 if (!set)
719 return NULL;
721 if (isl_set_is_empty(set))
722 convex_hull = isl_basic_set_empty(isl_space_copy(set->dim));
723 else
724 convex_hull = isl_basic_set_universe(isl_space_copy(set->dim));
725 isl_set_free(set);
726 return convex_hull;
729 /* Compute the convex hull of a pair of basic sets without any parameters or
730 * integer divisions using Fourier-Motzkin elimination.
731 * The convex hull is the set of all points that can be written as
732 * the sum of points from both basic sets (in homogeneous coordinates).
733 * We set up the constraints in a space with dimensions for each of
734 * the three sets and then project out the dimensions corresponding
735 * to the two original basic sets, retaining only those corresponding
736 * to the convex hull.
738 static __isl_give isl_basic_set *convex_hull_pair_elim(
739 __isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)
741 int i, j, k;
742 struct isl_basic_set *bset[2];
743 struct isl_basic_set *hull = NULL;
744 unsigned dim;
746 if (!bset1 || !bset2)
747 goto error;
749 dim = isl_basic_set_n_dim(bset1);
750 hull = isl_basic_set_alloc(bset1->ctx, 0, 2 + 3 * dim, 0,
751 1 + dim + bset1->n_eq + bset2->n_eq,
752 2 + bset1->n_ineq + bset2->n_ineq);
753 bset[0] = bset1;
754 bset[1] = bset2;
755 for (i = 0; i < 2; ++i) {
756 for (j = 0; j < bset[i]->n_eq; ++j) {
757 k = isl_basic_set_alloc_equality(hull);
758 if (k < 0)
759 goto error;
760 isl_seq_clr(hull->eq[k], (i+1) * (1+dim));
761 isl_seq_clr(hull->eq[k]+(i+2)*(1+dim), (1-i)*(1+dim));
762 isl_seq_cpy(hull->eq[k]+(i+1)*(1+dim), bset[i]->eq[j],
763 1+dim);
765 for (j = 0; j < bset[i]->n_ineq; ++j) {
766 k = isl_basic_set_alloc_inequality(hull);
767 if (k < 0)
768 goto error;
769 isl_seq_clr(hull->ineq[k], (i+1) * (1+dim));
770 isl_seq_clr(hull->ineq[k]+(i+2)*(1+dim), (1-i)*(1+dim));
771 isl_seq_cpy(hull->ineq[k]+(i+1)*(1+dim),
772 bset[i]->ineq[j], 1+dim);
774 k = isl_basic_set_alloc_inequality(hull);
775 if (k < 0)
776 goto error;
777 isl_seq_clr(hull->ineq[k], 1+2+3*dim);
778 isl_int_set_si(hull->ineq[k][(i+1)*(1+dim)], 1);
780 for (j = 0; j < 1+dim; ++j) {
781 k = isl_basic_set_alloc_equality(hull);
782 if (k < 0)
783 goto error;
784 isl_seq_clr(hull->eq[k], 1+2+3*dim);
785 isl_int_set_si(hull->eq[k][j], -1);
786 isl_int_set_si(hull->eq[k][1+dim+j], 1);
787 isl_int_set_si(hull->eq[k][2*(1+dim)+j], 1);
789 hull = isl_basic_set_set_rational(hull);
790 hull = isl_basic_set_remove_dims(hull, isl_dim_set, dim, 2*(1+dim));
791 hull = isl_basic_set_remove_redundancies(hull);
792 isl_basic_set_free(bset1);
793 isl_basic_set_free(bset2);
794 return hull;
795 error:
796 isl_basic_set_free(bset1);
797 isl_basic_set_free(bset2);
798 isl_basic_set_free(hull);
799 return NULL;
802 /* Is the set bounded for each value of the parameters?
804 isl_bool isl_basic_set_is_bounded(__isl_keep isl_basic_set *bset)
806 struct isl_tab *tab;
807 isl_bool bounded;
809 if (!bset)
810 return isl_bool_error;
811 if (isl_basic_set_plain_is_empty(bset))
812 return isl_bool_true;
814 tab = isl_tab_from_recession_cone(bset, 1);
815 bounded = isl_tab_cone_is_bounded(tab);
816 isl_tab_free(tab);
817 return bounded;
820 /* Is the image bounded for each value of the parameters and
821 * the domain variables?
823 isl_bool isl_basic_map_image_is_bounded(__isl_keep isl_basic_map *bmap)
825 unsigned nparam = isl_basic_map_dim(bmap, isl_dim_param);
826 unsigned n_in = isl_basic_map_dim(bmap, isl_dim_in);
827 isl_bool bounded;
829 bmap = isl_basic_map_copy(bmap);
830 bmap = isl_basic_map_cow(bmap);
831 bmap = isl_basic_map_move_dims(bmap, isl_dim_param, nparam,
832 isl_dim_in, 0, n_in);
833 bounded = isl_basic_set_is_bounded(bset_from_bmap(bmap));
834 isl_basic_map_free(bmap);
836 return bounded;
839 /* Is the set bounded for each value of the parameters?
841 isl_bool isl_set_is_bounded(__isl_keep isl_set *set)
843 int i;
845 if (!set)
846 return isl_bool_error;
848 for (i = 0; i < set->n; ++i) {
849 isl_bool bounded = isl_basic_set_is_bounded(set->p[i]);
850 if (!bounded || bounded < 0)
851 return bounded;
853 return isl_bool_true;
856 /* Compute the lineality space of the convex hull of bset1 and bset2.
858 * We first compute the intersection of the recession cone of bset1
859 * with the negative of the recession cone of bset2 and then compute
860 * the linear hull of the resulting cone.
862 static __isl_give isl_basic_set *induced_lineality_space(
863 __isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)
865 int i, k;
866 struct isl_basic_set *lin = NULL;
867 unsigned dim;
869 if (!bset1 || !bset2)
870 goto error;
872 dim = isl_basic_set_total_dim(bset1);
873 lin = isl_basic_set_alloc_space(isl_basic_set_get_space(bset1), 0,
874 bset1->n_eq + bset2->n_eq,
875 bset1->n_ineq + bset2->n_ineq);
876 lin = isl_basic_set_set_rational(lin);
877 if (!lin)
878 goto error;
879 for (i = 0; i < bset1->n_eq; ++i) {
880 k = isl_basic_set_alloc_equality(lin);
881 if (k < 0)
882 goto error;
883 isl_int_set_si(lin->eq[k][0], 0);
884 isl_seq_cpy(lin->eq[k] + 1, bset1->eq[i] + 1, dim);
886 for (i = 0; i < bset1->n_ineq; ++i) {
887 k = isl_basic_set_alloc_inequality(lin);
888 if (k < 0)
889 goto error;
890 isl_int_set_si(lin->ineq[k][0], 0);
891 isl_seq_cpy(lin->ineq[k] + 1, bset1->ineq[i] + 1, dim);
893 for (i = 0; i < bset2->n_eq; ++i) {
894 k = isl_basic_set_alloc_equality(lin);
895 if (k < 0)
896 goto error;
897 isl_int_set_si(lin->eq[k][0], 0);
898 isl_seq_neg(lin->eq[k] + 1, bset2->eq[i] + 1, dim);
900 for (i = 0; i < bset2->n_ineq; ++i) {
901 k = isl_basic_set_alloc_inequality(lin);
902 if (k < 0)
903 goto error;
904 isl_int_set_si(lin->ineq[k][0], 0);
905 isl_seq_neg(lin->ineq[k] + 1, bset2->ineq[i] + 1, dim);
908 isl_basic_set_free(bset1);
909 isl_basic_set_free(bset2);
910 return isl_basic_set_affine_hull(lin);
911 error:
912 isl_basic_set_free(lin);
913 isl_basic_set_free(bset1);
914 isl_basic_set_free(bset2);
915 return NULL;
918 static __isl_give isl_basic_set *uset_convex_hull(__isl_take isl_set *set);
920 /* Given a set and a linear space "lin" of dimension n > 0,
921 * project the linear space from the set, compute the convex hull
922 * and then map the set back to the original space.
924 * Let
926 * M x = 0
928 * describe the linear space. We first compute the Hermite normal
929 * form H = M U of M = H Q, to obtain
931 * H Q x = 0
933 * The last n rows of H will be zero, so the last n variables of x' = Q x
934 * are the one we want to project out. We do this by transforming each
935 * basic set A x >= b to A U x' >= b and then removing the last n dimensions.
936 * After computing the convex hull in x'_1, i.e., A' x'_1 >= b',
937 * we transform the hull back to the original space as A' Q_1 x >= b',
938 * with Q_1 all but the last n rows of Q.
940 static __isl_give isl_basic_set *modulo_lineality(__isl_take isl_set *set,
941 __isl_take isl_basic_set *lin)
943 unsigned total = isl_basic_set_total_dim(lin);
944 unsigned lin_dim;
945 struct isl_basic_set *hull;
946 struct isl_mat *M, *U, *Q;
948 if (!set || !lin)
949 goto error;
950 lin_dim = total - lin->n_eq;
951 M = isl_mat_sub_alloc6(set->ctx, lin->eq, 0, lin->n_eq, 1, total);
952 M = isl_mat_left_hermite(M, 0, &U, &Q);
953 if (!M)
954 goto error;
955 isl_mat_free(M);
956 isl_basic_set_free(lin);
958 Q = isl_mat_drop_rows(Q, Q->n_row - lin_dim, lin_dim);
960 U = isl_mat_lin_to_aff(U);
961 Q = isl_mat_lin_to_aff(Q);
963 set = isl_set_preimage(set, U);
964 set = isl_set_remove_dims(set, isl_dim_set, total - lin_dim, lin_dim);
965 hull = uset_convex_hull(set);
966 hull = isl_basic_set_preimage(hull, Q);
968 return hull;
969 error:
970 isl_basic_set_free(lin);
971 isl_set_free(set);
972 return NULL;
975 /* Given two polyhedra with as constraints h_{ij} x >= 0 in homegeneous space,
976 * set up an LP for solving
978 * \sum_j \alpha_{1j} h_{1j} = \sum_j \alpha_{2j} h_{2j}
980 * \alpha{i0} corresponds to the (implicit) positivity constraint 1 >= 0
981 * The next \alpha{ij} correspond to the equalities and come in pairs.
982 * The final \alpha{ij} correspond to the inequalities.
984 static __isl_give isl_basic_set *valid_direction_lp(
985 __isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)
987 isl_space *dim;
988 struct isl_basic_set *lp;
989 unsigned d;
990 int n;
991 int i, j, k;
993 if (!bset1 || !bset2)
994 goto error;
995 d = 1 + isl_basic_set_total_dim(bset1);
996 n = 2 +
997 2 * bset1->n_eq + bset1->n_ineq + 2 * bset2->n_eq + bset2->n_ineq;
998 dim = isl_space_set_alloc(bset1->ctx, 0, n);
999 lp = isl_basic_set_alloc_space(dim, 0, d, n);
1000 if (!lp)
1001 goto error;
1002 for (i = 0; i < n; ++i) {
1003 k = isl_basic_set_alloc_inequality(lp);
1004 if (k < 0)
1005 goto error;
1006 isl_seq_clr(lp->ineq[k] + 1, n);
1007 isl_int_set_si(lp->ineq[k][0], -1);
1008 isl_int_set_si(lp->ineq[k][1 + i], 1);
1010 for (i = 0; i < d; ++i) {
1011 k = isl_basic_set_alloc_equality(lp);
1012 if (k < 0)
1013 goto error;
1014 n = 0;
1015 isl_int_set_si(lp->eq[k][n], 0); n++;
1016 /* positivity constraint 1 >= 0 */
1017 isl_int_set_si(lp->eq[k][n], i == 0); n++;
1018 for (j = 0; j < bset1->n_eq; ++j) {
1019 isl_int_set(lp->eq[k][n], bset1->eq[j][i]); n++;
1020 isl_int_neg(lp->eq[k][n], bset1->eq[j][i]); n++;
1022 for (j = 0; j < bset1->n_ineq; ++j) {
1023 isl_int_set(lp->eq[k][n], bset1->ineq[j][i]); n++;
1025 /* positivity constraint 1 >= 0 */
1026 isl_int_set_si(lp->eq[k][n], -(i == 0)); n++;
1027 for (j = 0; j < bset2->n_eq; ++j) {
1028 isl_int_neg(lp->eq[k][n], bset2->eq[j][i]); n++;
1029 isl_int_set(lp->eq[k][n], bset2->eq[j][i]); n++;
1031 for (j = 0; j < bset2->n_ineq; ++j) {
1032 isl_int_neg(lp->eq[k][n], bset2->ineq[j][i]); n++;
1035 lp = isl_basic_set_gauss(lp, NULL);
1036 isl_basic_set_free(bset1);
1037 isl_basic_set_free(bset2);
1038 return lp;
1039 error:
1040 isl_basic_set_free(bset1);
1041 isl_basic_set_free(bset2);
1042 return NULL;
1045 /* Compute a vector s in the homogeneous space such that <s, r> > 0
1046 * for all rays in the homogeneous space of the two cones that correspond
1047 * to the input polyhedra bset1 and bset2.
1049 * We compute s as a vector that satisfies
1051 * s = \sum_j \alpha_{ij} h_{ij} for i = 1,2 (*)
1053 * with h_{ij} the normals of the facets of polyhedron i
1054 * (including the "positivity constraint" 1 >= 0) and \alpha_{ij}
1055 * strictly positive numbers. For simplicity we impose \alpha_{ij} >= 1.
1056 * We first set up an LP with as variables the \alpha{ij}.
1057 * In this formulation, for each polyhedron i,
1058 * the first constraint is the positivity constraint, followed by pairs
1059 * of variables for the equalities, followed by variables for the inequalities.
1060 * We then simply pick a feasible solution and compute s using (*).
1062 * Note that we simply pick any valid direction and make no attempt
1063 * to pick a "good" or even the "best" valid direction.
1065 static __isl_give isl_vec *valid_direction(
1066 __isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)
1068 struct isl_basic_set *lp;
1069 struct isl_tab *tab;
1070 struct isl_vec *sample = NULL;
1071 struct isl_vec *dir;
1072 unsigned d;
1073 int i;
1074 int n;
1076 if (!bset1 || !bset2)
1077 goto error;
1078 lp = valid_direction_lp(isl_basic_set_copy(bset1),
1079 isl_basic_set_copy(bset2));
1080 tab = isl_tab_from_basic_set(lp, 0);
1081 sample = isl_tab_get_sample_value(tab);
1082 isl_tab_free(tab);
1083 isl_basic_set_free(lp);
1084 if (!sample)
1085 goto error;
1086 d = isl_basic_set_total_dim(bset1);
1087 dir = isl_vec_alloc(bset1->ctx, 1 + d);
1088 if (!dir)
1089 goto error;
1090 isl_seq_clr(dir->block.data + 1, dir->size - 1);
1091 n = 1;
1092 /* positivity constraint 1 >= 0 */
1093 isl_int_set(dir->block.data[0], sample->block.data[n]); n++;
1094 for (i = 0; i < bset1->n_eq; ++i) {
1095 isl_int_sub(sample->block.data[n],
1096 sample->block.data[n], sample->block.data[n+1]);
1097 isl_seq_combine(dir->block.data,
1098 bset1->ctx->one, dir->block.data,
1099 sample->block.data[n], bset1->eq[i], 1 + d);
1101 n += 2;
1103 for (i = 0; i < bset1->n_ineq; ++i)
1104 isl_seq_combine(dir->block.data,
1105 bset1->ctx->one, dir->block.data,
1106 sample->block.data[n++], bset1->ineq[i], 1 + d);
1107 isl_vec_free(sample);
1108 isl_seq_normalize(bset1->ctx, dir->el, dir->size);
1109 isl_basic_set_free(bset1);
1110 isl_basic_set_free(bset2);
1111 return dir;
1112 error:
1113 isl_vec_free(sample);
1114 isl_basic_set_free(bset1);
1115 isl_basic_set_free(bset2);
1116 return NULL;
1119 /* Given a polyhedron b_i + A_i x >= 0 and a map T = S^{-1},
1120 * compute b_i' + A_i' x' >= 0, with
1122 * [ b_i A_i ] [ y' ] [ y' ]
1123 * [ 1 0 ] S^{-1} [ x' ] >= 0 or [ b_i' A_i' ] [ x' ] >= 0
1125 * In particular, add the "positivity constraint" and then perform
1126 * the mapping.
1128 static __isl_give isl_basic_set *homogeneous_map(__isl_take isl_basic_set *bset,
1129 __isl_take isl_mat *T)
1131 int k;
1133 if (!bset)
1134 goto error;
1135 bset = isl_basic_set_extend_constraints(bset, 0, 1);
1136 k = isl_basic_set_alloc_inequality(bset);
1137 if (k < 0)
1138 goto error;
1139 isl_seq_clr(bset->ineq[k] + 1, isl_basic_set_total_dim(bset));
1140 isl_int_set_si(bset->ineq[k][0], 1);
1141 bset = isl_basic_set_preimage(bset, T);
1142 return bset;
1143 error:
1144 isl_mat_free(T);
1145 isl_basic_set_free(bset);
1146 return NULL;
1149 /* Compute the convex hull of a pair of basic sets without any parameters or
1150 * integer divisions, where the convex hull is known to be pointed,
1151 * but the basic sets may be unbounded.
1153 * We turn this problem into the computation of a convex hull of a pair
1154 * _bounded_ polyhedra by "changing the direction of the homogeneous
1155 * dimension". This idea is due to Matthias Koeppe.
1157 * Consider the cones in homogeneous space that correspond to the
1158 * input polyhedra. The rays of these cones are also rays of the
1159 * polyhedra if the coordinate that corresponds to the homogeneous
1160 * dimension is zero. That is, if the inner product of the rays
1161 * with the homogeneous direction is zero.
1162 * The cones in the homogeneous space can also be considered to
1163 * correspond to other pairs of polyhedra by chosing a different
1164 * homogeneous direction. To ensure that both of these polyhedra
1165 * are bounded, we need to make sure that all rays of the cones
1166 * correspond to vertices and not to rays.
1167 * Let s be a direction such that <s, r> > 0 for all rays r of both cones.
1168 * Then using s as a homogeneous direction, we obtain a pair of polytopes.
1169 * The vector s is computed in valid_direction.
1171 * Note that we need to consider _all_ rays of the cones and not just
1172 * the rays that correspond to rays in the polyhedra. If we were to
1173 * only consider those rays and turn them into vertices, then we
1174 * may inadvertently turn some vertices into rays.
1176 * The standard homogeneous direction is the unit vector in the 0th coordinate.
1177 * We therefore transform the two polyhedra such that the selected
1178 * direction is mapped onto this standard direction and then proceed
1179 * with the normal computation.
1180 * Let S be a non-singular square matrix with s as its first row,
1181 * then we want to map the polyhedra to the space
1183 * [ y' ] [ y ] [ y ] [ y' ]
1184 * [ x' ] = S [ x ] i.e., [ x ] = S^{-1} [ x' ]
1186 * We take S to be the unimodular completion of s to limit the growth
1187 * of the coefficients in the following computations.
1189 * Let b_i + A_i x >= 0 be the constraints of polyhedron i.
1190 * We first move to the homogeneous dimension
1192 * b_i y + A_i x >= 0 [ b_i A_i ] [ y ] [ 0 ]
1193 * y >= 0 or [ 1 0 ] [ x ] >= [ 0 ]
1195 * Then we change directoin
1197 * [ b_i A_i ] [ y' ] [ y' ]
1198 * [ 1 0 ] S^{-1} [ x' ] >= 0 or [ b_i' A_i' ] [ x' ] >= 0
1200 * Then we compute the convex hull of the polytopes b_i' + A_i' x' >= 0
1201 * resulting in b' + A' x' >= 0, which we then convert back
1203 * [ y ] [ y ]
1204 * [ b' A' ] S [ x ] >= 0 or [ b A ] [ x ] >= 0
1206 * The polyhedron b + A x >= 0 is then the convex hull of the input polyhedra.
1208 static __isl_give isl_basic_set *convex_hull_pair_pointed(
1209 __isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)
1211 struct isl_ctx *ctx = NULL;
1212 struct isl_vec *dir = NULL;
1213 struct isl_mat *T = NULL;
1214 struct isl_mat *T2 = NULL;
1215 struct isl_basic_set *hull;
1216 struct isl_set *set;
1218 if (!bset1 || !bset2)
1219 goto error;
1220 ctx = isl_basic_set_get_ctx(bset1);
1221 dir = valid_direction(isl_basic_set_copy(bset1),
1222 isl_basic_set_copy(bset2));
1223 if (!dir)
1224 goto error;
1225 T = isl_mat_alloc(ctx, dir->size, dir->size);
1226 if (!T)
1227 goto error;
1228 isl_seq_cpy(T->row[0], dir->block.data, dir->size);
1229 T = isl_mat_unimodular_complete(T, 1);
1230 T2 = isl_mat_right_inverse(isl_mat_copy(T));
1232 bset1 = homogeneous_map(bset1, isl_mat_copy(T2));
1233 bset2 = homogeneous_map(bset2, T2);
1234 set = isl_set_alloc_space(isl_basic_set_get_space(bset1), 2, 0);
1235 set = isl_set_add_basic_set(set, bset1);
1236 set = isl_set_add_basic_set(set, bset2);
1237 hull = uset_convex_hull(set);
1238 hull = isl_basic_set_preimage(hull, T);
1240 isl_vec_free(dir);
1242 return hull;
1243 error:
1244 isl_vec_free(dir);
1245 isl_basic_set_free(bset1);
1246 isl_basic_set_free(bset2);
1247 return NULL;
1250 static __isl_give isl_basic_set *uset_convex_hull_wrap(__isl_take isl_set *set);
1251 static __isl_give isl_basic_set *modulo_affine_hull(
1252 __isl_take isl_set *set, __isl_take isl_basic_set *affine_hull);
1254 /* Compute the convex hull of a pair of basic sets without any parameters or
1255 * integer divisions.
1257 * This function is called from uset_convex_hull_unbounded, which
1258 * means that the complete convex hull is unbounded. Some pairs
1259 * of basic sets may still be bounded, though.
1260 * They may even lie inside a lower dimensional space, in which
1261 * case they need to be handled inside their affine hull since
1262 * the main algorithm assumes that the result is full-dimensional.
1264 * If the convex hull of the two basic sets would have a non-trivial
1265 * lineality space, we first project out this lineality space.
1267 static __isl_give isl_basic_set *convex_hull_pair(
1268 __isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)
1270 isl_basic_set *lin, *aff;
1271 isl_bool bounded1, bounded2;
1273 if (bset1->ctx->opt->convex == ISL_CONVEX_HULL_FM)
1274 return convex_hull_pair_elim(bset1, bset2);
1276 aff = isl_set_affine_hull(isl_basic_set_union(isl_basic_set_copy(bset1),
1277 isl_basic_set_copy(bset2)));
1278 if (!aff)
1279 goto error;
1280 if (aff->n_eq != 0)
1281 return modulo_affine_hull(isl_basic_set_union(bset1, bset2), aff);
1282 isl_basic_set_free(aff);
1284 bounded1 = isl_basic_set_is_bounded(bset1);
1285 bounded2 = isl_basic_set_is_bounded(bset2);
1287 if (bounded1 < 0 || bounded2 < 0)
1288 goto error;
1290 if (bounded1 && bounded2)
1291 return uset_convex_hull_wrap(isl_basic_set_union(bset1, bset2));
1293 if (bounded1 || bounded2)
1294 return convex_hull_pair_pointed(bset1, bset2);
1296 lin = induced_lineality_space(isl_basic_set_copy(bset1),
1297 isl_basic_set_copy(bset2));
1298 if (!lin)
1299 goto error;
1300 if (isl_basic_set_plain_is_universe(lin)) {
1301 isl_basic_set_free(bset1);
1302 isl_basic_set_free(bset2);
1303 return lin;
1305 if (lin->n_eq < isl_basic_set_total_dim(lin)) {
1306 struct isl_set *set;
1307 set = isl_set_alloc_space(isl_basic_set_get_space(bset1), 2, 0);
1308 set = isl_set_add_basic_set(set, bset1);
1309 set = isl_set_add_basic_set(set, bset2);
1310 return modulo_lineality(set, lin);
1312 isl_basic_set_free(lin);
1314 return convex_hull_pair_pointed(bset1, bset2);
1315 error:
1316 isl_basic_set_free(bset1);
1317 isl_basic_set_free(bset2);
1318 return NULL;
1321 /* Compute the lineality space of a basic set.
1322 * We basically just drop the constants and turn every inequality
1323 * into an equality.
1324 * Any explicit representations of local variables are removed
1325 * because they may no longer be valid representations
1326 * in the lineality space.
1328 __isl_give isl_basic_set *isl_basic_set_lineality_space(
1329 __isl_take isl_basic_set *bset)
1331 int i, k;
1332 struct isl_basic_set *lin = NULL;
1333 unsigned n_div, dim;
1335 if (!bset)
1336 goto error;
1337 n_div = isl_basic_set_dim(bset, isl_dim_div);
1338 dim = isl_basic_set_total_dim(bset);
1340 lin = isl_basic_set_alloc_space(isl_basic_set_get_space(bset),
1341 n_div, dim, 0);
1342 for (i = 0; i < n_div; ++i)
1343 if (isl_basic_set_alloc_div(lin) < 0)
1344 goto error;
1345 if (!lin)
1346 goto error;
1347 for (i = 0; i < bset->n_eq; ++i) {
1348 k = isl_basic_set_alloc_equality(lin);
1349 if (k < 0)
1350 goto error;
1351 isl_int_set_si(lin->eq[k][0], 0);
1352 isl_seq_cpy(lin->eq[k] + 1, bset->eq[i] + 1, dim);
1354 lin = isl_basic_set_gauss(lin, NULL);
1355 if (!lin)
1356 goto error;
1357 for (i = 0; i < bset->n_ineq && lin->n_eq < dim; ++i) {
1358 k = isl_basic_set_alloc_equality(lin);
1359 if (k < 0)
1360 goto error;
1361 isl_int_set_si(lin->eq[k][0], 0);
1362 isl_seq_cpy(lin->eq[k] + 1, bset->ineq[i] + 1, dim);
1363 lin = isl_basic_set_gauss(lin, NULL);
1364 if (!lin)
1365 goto error;
1367 isl_basic_set_free(bset);
1368 return lin;
1369 error:
1370 isl_basic_set_free(lin);
1371 isl_basic_set_free(bset);
1372 return NULL;
1375 /* Compute the (linear) hull of the lineality spaces of the basic sets in the
1376 * set "set".
1378 __isl_give isl_basic_set *isl_set_combined_lineality_space(
1379 __isl_take isl_set *set)
1381 int i;
1382 struct isl_set *lin = NULL;
1384 if (!set)
1385 return NULL;
1386 if (set->n == 0) {
1387 isl_space *space = isl_set_get_space(set);
1388 isl_set_free(set);
1389 return isl_basic_set_empty(space);
1392 lin = isl_set_alloc_space(isl_set_get_space(set), set->n, 0);
1393 for (i = 0; i < set->n; ++i)
1394 lin = isl_set_add_basic_set(lin,
1395 isl_basic_set_lineality_space(isl_basic_set_copy(set->p[i])));
1396 isl_set_free(set);
1397 return isl_set_affine_hull(lin);
1400 /* Compute the convex hull of a set without any parameters or
1401 * integer divisions.
1402 * In each step, we combined two basic sets until only one
1403 * basic set is left.
1404 * The input basic sets are assumed not to have a non-trivial
1405 * lineality space. If any of the intermediate results has
1406 * a non-trivial lineality space, it is projected out.
1408 static __isl_give isl_basic_set *uset_convex_hull_unbounded(
1409 __isl_take isl_set *set)
1411 isl_basic_set_list *list;
1413 list = isl_set_get_basic_set_list(set);
1414 isl_set_free(set);
1416 while (list) {
1417 int n;
1418 struct isl_basic_set *t;
1419 isl_basic_set *bset1, *bset2;
1421 n = isl_basic_set_list_n_basic_set(list);
1422 if (n < 2)
1423 isl_die(isl_basic_set_list_get_ctx(list),
1424 isl_error_internal,
1425 "expecting at least two elements", goto error);
1426 bset1 = isl_basic_set_list_get_basic_set(list, n - 1);
1427 bset2 = isl_basic_set_list_get_basic_set(list, n - 2);
1428 bset1 = convex_hull_pair(bset1, bset2);
1429 if (n == 2) {
1430 isl_basic_set_list_free(list);
1431 return bset1;
1433 bset1 = isl_basic_set_underlying_set(bset1);
1434 list = isl_basic_set_list_drop(list, n - 2, 2);
1435 list = isl_basic_set_list_add(list, bset1);
1437 t = isl_basic_set_list_get_basic_set(list, n - 2);
1438 t = isl_basic_set_lineality_space(t);
1439 if (!t)
1440 goto error;
1441 if (isl_basic_set_plain_is_universe(t)) {
1442 isl_basic_set_list_free(list);
1443 return t;
1445 if (t->n_eq < isl_basic_set_total_dim(t)) {
1446 set = isl_basic_set_list_union(list);
1447 return modulo_lineality(set, t);
1449 isl_basic_set_free(t);
1452 return NULL;
1453 error:
1454 isl_basic_set_list_free(list);
1455 return NULL;
1458 /* Compute an initial hull for wrapping containing a single initial
1459 * facet.
1460 * This function assumes that the given set is bounded.
1462 static __isl_give isl_basic_set *initial_hull(__isl_take isl_basic_set *hull,
1463 __isl_keep isl_set *set)
1465 struct isl_mat *bounds = NULL;
1466 unsigned dim;
1467 int k;
1469 if (!hull)
1470 goto error;
1471 bounds = initial_facet_constraint(set);
1472 if (!bounds)
1473 goto error;
1474 k = isl_basic_set_alloc_inequality(hull);
1475 if (k < 0)
1476 goto error;
1477 dim = isl_set_n_dim(set);
1478 isl_assert(set->ctx, 1 + dim == bounds->n_col, goto error);
1479 isl_seq_cpy(hull->ineq[k], bounds->row[0], bounds->n_col);
1480 isl_mat_free(bounds);
1482 return hull;
1483 error:
1484 isl_basic_set_free(hull);
1485 isl_mat_free(bounds);
1486 return NULL;
1489 struct max_constraint {
1490 struct isl_mat *c;
1491 int count;
1492 int ineq;
1495 static int max_constraint_equal(const void *entry, const void *val)
1497 struct max_constraint *a = (struct max_constraint *)entry;
1498 isl_int *b = (isl_int *)val;
1500 return isl_seq_eq(a->c->row[0] + 1, b, a->c->n_col - 1);
1503 static void update_constraint(struct isl_ctx *ctx, struct isl_hash_table *table,
1504 isl_int *con, unsigned len, int n, int ineq)
1506 struct isl_hash_table_entry *entry;
1507 struct max_constraint *c;
1508 uint32_t c_hash;
1510 c_hash = isl_seq_get_hash(con + 1, len);
1511 entry = isl_hash_table_find(ctx, table, c_hash, max_constraint_equal,
1512 con + 1, 0);
1513 if (!entry)
1514 return;
1515 c = entry->data;
1516 if (c->count < n) {
1517 isl_hash_table_remove(ctx, table, entry);
1518 return;
1520 c->count++;
1521 if (isl_int_gt(c->c->row[0][0], con[0]))
1522 return;
1523 if (isl_int_eq(c->c->row[0][0], con[0])) {
1524 if (ineq)
1525 c->ineq = ineq;
1526 return;
1528 c->c = isl_mat_cow(c->c);
1529 isl_int_set(c->c->row[0][0], con[0]);
1530 c->ineq = ineq;
1533 /* Check whether the constraint hash table "table" contains the constraint
1534 * "con".
1536 static int has_constraint(struct isl_ctx *ctx, struct isl_hash_table *table,
1537 isl_int *con, unsigned len, int n)
1539 struct isl_hash_table_entry *entry;
1540 struct max_constraint *c;
1541 uint32_t c_hash;
1543 c_hash = isl_seq_get_hash(con + 1, len);
1544 entry = isl_hash_table_find(ctx, table, c_hash, max_constraint_equal,
1545 con + 1, 0);
1546 if (!entry)
1547 return 0;
1548 c = entry->data;
1549 if (c->count < n)
1550 return 0;
1551 return isl_int_eq(c->c->row[0][0], con[0]);
1554 /* Check for inequality constraints of a basic set without equalities
1555 * such that the same or more stringent copies of the constraint appear
1556 * in all of the basic sets. Such constraints are necessarily facet
1557 * constraints of the convex hull.
1559 * If the resulting basic set is by chance identical to one of
1560 * the basic sets in "set", then we know that this basic set contains
1561 * all other basic sets and is therefore the convex hull of set.
1562 * In this case we set *is_hull to 1.
1564 static __isl_give isl_basic_set *common_constraints(
1565 __isl_take isl_basic_set *hull, __isl_keep isl_set *set, int *is_hull)
1567 int i, j, s, n;
1568 int min_constraints;
1569 int best;
1570 struct max_constraint *constraints = NULL;
1571 struct isl_hash_table *table = NULL;
1572 unsigned total;
1574 *is_hull = 0;
1576 for (i = 0; i < set->n; ++i)
1577 if (set->p[i]->n_eq == 0)
1578 break;
1579 if (i >= set->n)
1580 return hull;
1581 min_constraints = set->p[i]->n_ineq;
1582 best = i;
1583 for (i = best + 1; i < set->n; ++i) {
1584 if (set->p[i]->n_eq != 0)
1585 continue;
1586 if (set->p[i]->n_ineq >= min_constraints)
1587 continue;
1588 min_constraints = set->p[i]->n_ineq;
1589 best = i;
1591 constraints = isl_calloc_array(hull->ctx, struct max_constraint,
1592 min_constraints);
1593 if (!constraints)
1594 return hull;
1595 table = isl_alloc_type(hull->ctx, struct isl_hash_table);
1596 if (isl_hash_table_init(hull->ctx, table, min_constraints))
1597 goto error;
1599 total = isl_space_dim(set->dim, isl_dim_all);
1600 for (i = 0; i < set->p[best]->n_ineq; ++i) {
1601 constraints[i].c = isl_mat_sub_alloc6(hull->ctx,
1602 set->p[best]->ineq + i, 0, 1, 0, 1 + total);
1603 if (!constraints[i].c)
1604 goto error;
1605 constraints[i].ineq = 1;
1607 for (i = 0; i < min_constraints; ++i) {
1608 struct isl_hash_table_entry *entry;
1609 uint32_t c_hash;
1610 c_hash = isl_seq_get_hash(constraints[i].c->row[0] + 1, total);
1611 entry = isl_hash_table_find(hull->ctx, table, c_hash,
1612 max_constraint_equal, constraints[i].c->row[0] + 1, 1);
1613 if (!entry)
1614 goto error;
1615 isl_assert(hull->ctx, !entry->data, goto error);
1616 entry->data = &constraints[i];
1619 n = 0;
1620 for (s = 0; s < set->n; ++s) {
1621 if (s == best)
1622 continue;
1624 for (i = 0; i < set->p[s]->n_eq; ++i) {
1625 isl_int *eq = set->p[s]->eq[i];
1626 for (j = 0; j < 2; ++j) {
1627 isl_seq_neg(eq, eq, 1 + total);
1628 update_constraint(hull->ctx, table,
1629 eq, total, n, 0);
1632 for (i = 0; i < set->p[s]->n_ineq; ++i) {
1633 isl_int *ineq = set->p[s]->ineq[i];
1634 update_constraint(hull->ctx, table, ineq, total, n,
1635 set->p[s]->n_eq == 0);
1637 ++n;
1640 for (i = 0; i < min_constraints; ++i) {
1641 if (constraints[i].count < n)
1642 continue;
1643 if (!constraints[i].ineq)
1644 continue;
1645 j = isl_basic_set_alloc_inequality(hull);
1646 if (j < 0)
1647 goto error;
1648 isl_seq_cpy(hull->ineq[j], constraints[i].c->row[0], 1 + total);
1651 for (s = 0; s < set->n; ++s) {
1652 if (set->p[s]->n_eq)
1653 continue;
1654 if (set->p[s]->n_ineq != hull->n_ineq)
1655 continue;
1656 for (i = 0; i < set->p[s]->n_ineq; ++i) {
1657 isl_int *ineq = set->p[s]->ineq[i];
1658 if (!has_constraint(hull->ctx, table, ineq, total, n))
1659 break;
1661 if (i == set->p[s]->n_ineq)
1662 *is_hull = 1;
1665 isl_hash_table_clear(table);
1666 for (i = 0; i < min_constraints; ++i)
1667 isl_mat_free(constraints[i].c);
1668 free(constraints);
1669 free(table);
1670 return hull;
1671 error:
1672 isl_hash_table_clear(table);
1673 free(table);
1674 if (constraints)
1675 for (i = 0; i < min_constraints; ++i)
1676 isl_mat_free(constraints[i].c);
1677 free(constraints);
1678 return hull;
1681 /* Create a template for the convex hull of "set" and fill it up
1682 * obvious facet constraints, if any. If the result happens to
1683 * be the convex hull of "set" then *is_hull is set to 1.
1685 static __isl_give isl_basic_set *proto_hull(__isl_keep isl_set *set,
1686 int *is_hull)
1688 struct isl_basic_set *hull;
1689 unsigned n_ineq;
1690 int i;
1692 n_ineq = 1;
1693 for (i = 0; i < set->n; ++i) {
1694 n_ineq += set->p[i]->n_eq;
1695 n_ineq += set->p[i]->n_ineq;
1697 hull = isl_basic_set_alloc_space(isl_space_copy(set->dim), 0, 0, n_ineq);
1698 hull = isl_basic_set_set_rational(hull);
1699 if (!hull)
1700 return NULL;
1701 return common_constraints(hull, set, is_hull);
1704 static __isl_give isl_basic_set *uset_convex_hull_wrap(__isl_take isl_set *set)
1706 struct isl_basic_set *hull;
1707 int is_hull;
1709 hull = proto_hull(set, &is_hull);
1710 if (hull && !is_hull) {
1711 if (hull->n_ineq == 0)
1712 hull = initial_hull(hull, set);
1713 hull = extend(hull, set);
1715 isl_set_free(set);
1717 return hull;
1720 /* Compute the convex hull of a set without any parameters or
1721 * integer divisions. Depending on whether the set is bounded,
1722 * we pass control to the wrapping based convex hull or
1723 * the Fourier-Motzkin elimination based convex hull.
1724 * We also handle a few special cases before checking the boundedness.
1726 static __isl_give isl_basic_set *uset_convex_hull(__isl_take isl_set *set)
1728 isl_bool bounded;
1729 struct isl_basic_set *convex_hull = NULL;
1730 struct isl_basic_set *lin;
1732 if (isl_set_n_dim(set) == 0)
1733 return convex_hull_0d(set);
1735 set = isl_set_coalesce(set);
1736 set = isl_set_set_rational(set);
1738 if (!set)
1739 return NULL;
1740 if (set->n == 1) {
1741 convex_hull = isl_basic_set_copy(set->p[0]);
1742 isl_set_free(set);
1743 return convex_hull;
1745 if (isl_set_n_dim(set) == 1)
1746 return convex_hull_1d(set);
1748 bounded = isl_set_is_bounded(set);
1749 if (bounded < 0)
1750 goto error;
1751 if (bounded && set->ctx->opt->convex == ISL_CONVEX_HULL_WRAP)
1752 return uset_convex_hull_wrap(set);
1754 lin = isl_set_combined_lineality_space(isl_set_copy(set));
1755 if (!lin)
1756 goto error;
1757 if (isl_basic_set_plain_is_universe(lin)) {
1758 isl_set_free(set);
1759 return lin;
1761 if (lin->n_eq < isl_basic_set_total_dim(lin))
1762 return modulo_lineality(set, lin);
1763 isl_basic_set_free(lin);
1765 return uset_convex_hull_unbounded(set);
1766 error:
1767 isl_set_free(set);
1768 isl_basic_set_free(convex_hull);
1769 return NULL;
1772 /* This is the core procedure, where "set" is a "pure" set, i.e.,
1773 * without parameters or divs and where the convex hull of set is
1774 * known to be full-dimensional.
1776 static __isl_give isl_basic_set *uset_convex_hull_wrap_bounded(
1777 __isl_take isl_set *set)
1779 struct isl_basic_set *convex_hull = NULL;
1781 if (!set)
1782 goto error;
1784 if (isl_set_n_dim(set) == 0) {
1785 convex_hull = isl_basic_set_universe(isl_space_copy(set->dim));
1786 isl_set_free(set);
1787 convex_hull = isl_basic_set_set_rational(convex_hull);
1788 return convex_hull;
1791 set = isl_set_set_rational(set);
1792 set = isl_set_coalesce(set);
1793 if (!set)
1794 goto error;
1795 if (set->n == 1) {
1796 convex_hull = isl_basic_set_copy(set->p[0]);
1797 isl_set_free(set);
1798 convex_hull = isl_basic_map_remove_redundancies(convex_hull);
1799 return convex_hull;
1801 if (isl_set_n_dim(set) == 1)
1802 return convex_hull_1d(set);
1804 return uset_convex_hull_wrap(set);
1805 error:
1806 isl_set_free(set);
1807 return NULL;
1810 /* Compute the convex hull of set "set" with affine hull "affine_hull",
1811 * We first remove the equalities (transforming the set), compute the
1812 * convex hull of the transformed set and then add the equalities back
1813 * (after performing the inverse transformation.
1815 static __isl_give isl_basic_set *modulo_affine_hull(
1816 __isl_take isl_set *set, __isl_take isl_basic_set *affine_hull)
1818 struct isl_mat *T;
1819 struct isl_mat *T2;
1820 struct isl_basic_set *dummy;
1821 struct isl_basic_set *convex_hull;
1823 dummy = isl_basic_set_remove_equalities(
1824 isl_basic_set_copy(affine_hull), &T, &T2);
1825 if (!dummy)
1826 goto error;
1827 isl_basic_set_free(dummy);
1828 set = isl_set_preimage(set, T);
1829 convex_hull = uset_convex_hull(set);
1830 convex_hull = isl_basic_set_preimage(convex_hull, T2);
1831 convex_hull = isl_basic_set_intersect(convex_hull, affine_hull);
1832 return convex_hull;
1833 error:
1834 isl_mat_free(T);
1835 isl_mat_free(T2);
1836 isl_basic_set_free(affine_hull);
1837 isl_set_free(set);
1838 return NULL;
1841 /* Return an empty basic map living in the same space as "map".
1843 static __isl_give isl_basic_map *replace_map_by_empty_basic_map(
1844 __isl_take isl_map *map)
1846 isl_space *space;
1848 space = isl_map_get_space(map);
1849 isl_map_free(map);
1850 return isl_basic_map_empty(space);
1853 /* Compute the convex hull of a map.
1855 * The implementation was inspired by "Extended Convex Hull" by Fukuda et al.,
1856 * specifically, the wrapping of facets to obtain new facets.
1858 __isl_give isl_basic_map *isl_map_convex_hull(__isl_take isl_map *map)
1860 struct isl_basic_set *bset;
1861 struct isl_basic_map *model = NULL;
1862 struct isl_basic_set *affine_hull = NULL;
1863 struct isl_basic_map *convex_hull = NULL;
1864 struct isl_set *set = NULL;
1866 map = isl_map_detect_equalities(map);
1867 map = isl_map_align_divs_internal(map);
1868 if (!map)
1869 goto error;
1871 if (map->n == 0)
1872 return replace_map_by_empty_basic_map(map);
1874 model = isl_basic_map_copy(map->p[0]);
1875 set = isl_map_underlying_set(map);
1876 if (!set)
1877 goto error;
1879 affine_hull = isl_set_affine_hull(isl_set_copy(set));
1880 if (!affine_hull)
1881 goto error;
1882 if (affine_hull->n_eq != 0)
1883 bset = modulo_affine_hull(set, affine_hull);
1884 else {
1885 isl_basic_set_free(affine_hull);
1886 bset = uset_convex_hull(set);
1889 convex_hull = isl_basic_map_overlying_set(bset, model);
1890 if (!convex_hull)
1891 return NULL;
1893 ISL_F_SET(convex_hull, ISL_BASIC_MAP_NO_IMPLICIT);
1894 ISL_F_SET(convex_hull, ISL_BASIC_MAP_ALL_EQUALITIES);
1895 ISL_F_CLR(convex_hull, ISL_BASIC_MAP_RATIONAL);
1896 return convex_hull;
1897 error:
1898 isl_set_free(set);
1899 isl_basic_map_free(model);
1900 return NULL;
1903 struct isl_basic_set *isl_set_convex_hull(struct isl_set *set)
1905 return bset_from_bmap(isl_map_convex_hull(set_to_map(set)));
1908 __isl_give isl_basic_map *isl_map_polyhedral_hull(__isl_take isl_map *map)
1910 isl_basic_map *hull;
1912 hull = isl_map_convex_hull(map);
1913 return isl_basic_map_remove_divs(hull);
1916 __isl_give isl_basic_set *isl_set_polyhedral_hull(__isl_take isl_set *set)
1918 return bset_from_bmap(isl_map_polyhedral_hull(set_to_map(set)));
1921 struct sh_data_entry {
1922 struct isl_hash_table *table;
1923 struct isl_tab *tab;
1926 /* Holds the data needed during the simple hull computation.
1927 * In particular,
1928 * n the number of basic sets in the original set
1929 * hull_table a hash table of already computed constraints
1930 * in the simple hull
1931 * p for each basic set,
1932 * table a hash table of the constraints
1933 * tab the tableau corresponding to the basic set
1935 struct sh_data {
1936 struct isl_ctx *ctx;
1937 unsigned n;
1938 struct isl_hash_table *hull_table;
1939 struct sh_data_entry p[1];
1942 static void sh_data_free(struct sh_data *data)
1944 int i;
1946 if (!data)
1947 return;
1948 isl_hash_table_free(data->ctx, data->hull_table);
1949 for (i = 0; i < data->n; ++i) {
1950 isl_hash_table_free(data->ctx, data->p[i].table);
1951 isl_tab_free(data->p[i].tab);
1953 free(data);
1956 struct ineq_cmp_data {
1957 unsigned len;
1958 isl_int *p;
1961 static int has_ineq(const void *entry, const void *val)
1963 isl_int *row = (isl_int *)entry;
1964 struct ineq_cmp_data *v = (struct ineq_cmp_data *)val;
1966 return isl_seq_eq(row + 1, v->p + 1, v->len) ||
1967 isl_seq_is_neg(row + 1, v->p + 1, v->len);
1970 static int hash_ineq(struct isl_ctx *ctx, struct isl_hash_table *table,
1971 isl_int *ineq, unsigned len)
1973 uint32_t c_hash;
1974 struct ineq_cmp_data v;
1975 struct isl_hash_table_entry *entry;
1977 v.len = len;
1978 v.p = ineq;
1979 c_hash = isl_seq_get_hash(ineq + 1, len);
1980 entry = isl_hash_table_find(ctx, table, c_hash, has_ineq, &v, 1);
1981 if (!entry)
1982 return - 1;
1983 entry->data = ineq;
1984 return 0;
1987 /* Fill hash table "table" with the constraints of "bset".
1988 * Equalities are added as two inequalities.
1989 * The value in the hash table is a pointer to the (in)equality of "bset".
1991 static int hash_basic_set(struct isl_hash_table *table,
1992 __isl_keep isl_basic_set *bset)
1994 int i, j;
1995 unsigned dim = isl_basic_set_total_dim(bset);
1997 for (i = 0; i < bset->n_eq; ++i) {
1998 for (j = 0; j < 2; ++j) {
1999 isl_seq_neg(bset->eq[i], bset->eq[i], 1 + dim);
2000 if (hash_ineq(bset->ctx, table, bset->eq[i], dim) < 0)
2001 return -1;
2004 for (i = 0; i < bset->n_ineq; ++i) {
2005 if (hash_ineq(bset->ctx, table, bset->ineq[i], dim) < 0)
2006 return -1;
2008 return 0;
2011 static struct sh_data *sh_data_alloc(__isl_keep isl_set *set, unsigned n_ineq)
2013 struct sh_data *data;
2014 int i;
2016 data = isl_calloc(set->ctx, struct sh_data,
2017 sizeof(struct sh_data) +
2018 (set->n - 1) * sizeof(struct sh_data_entry));
2019 if (!data)
2020 return NULL;
2021 data->ctx = set->ctx;
2022 data->n = set->n;
2023 data->hull_table = isl_hash_table_alloc(set->ctx, n_ineq);
2024 if (!data->hull_table)
2025 goto error;
2026 for (i = 0; i < set->n; ++i) {
2027 data->p[i].table = isl_hash_table_alloc(set->ctx,
2028 2 * set->p[i]->n_eq + set->p[i]->n_ineq);
2029 if (!data->p[i].table)
2030 goto error;
2031 if (hash_basic_set(data->p[i].table, set->p[i]) < 0)
2032 goto error;
2034 return data;
2035 error:
2036 sh_data_free(data);
2037 return NULL;
2040 /* Check if inequality "ineq" is a bound for basic set "j" or if
2041 * it can be relaxed (by increasing the constant term) to become
2042 * a bound for that basic set. In the latter case, the constant
2043 * term is updated.
2044 * Relaxation of the constant term is only allowed if "shift" is set.
2046 * Return 1 if "ineq" is a bound
2047 * 0 if "ineq" may attain arbitrarily small values on basic set "j"
2048 * -1 if some error occurred
2050 static int is_bound(struct sh_data *data, __isl_keep isl_set *set, int j,
2051 isl_int *ineq, int shift)
2053 enum isl_lp_result res;
2054 isl_int opt;
2056 if (!data->p[j].tab) {
2057 data->p[j].tab = isl_tab_from_basic_set(set->p[j], 0);
2058 if (!data->p[j].tab)
2059 return -1;
2062 isl_int_init(opt);
2064 res = isl_tab_min(data->p[j].tab, ineq, data->ctx->one,
2065 &opt, NULL, 0);
2066 if (res == isl_lp_ok && isl_int_is_neg(opt)) {
2067 if (shift)
2068 isl_int_sub(ineq[0], ineq[0], opt);
2069 else
2070 res = isl_lp_unbounded;
2073 isl_int_clear(opt);
2075 return (res == isl_lp_ok || res == isl_lp_empty) ? 1 :
2076 res == isl_lp_unbounded ? 0 : -1;
2079 /* Set the constant term of "ineq" to the maximum of those of the constraints
2080 * in the basic sets of "set" following "i" that are parallel to "ineq".
2081 * That is, if any of the basic sets of "set" following "i" have a more
2082 * relaxed copy of "ineq", then replace "ineq" by the most relaxed copy.
2083 * "c_hash" is the hash value of the linear part of "ineq".
2084 * "v" has been set up for use by has_ineq.
2086 * Note that the two inequality constraints corresponding to an equality are
2087 * represented by the same inequality constraint in data->p[j].table
2088 * (but with different hash values). This means the constraint (or at
2089 * least its constant term) may need to be temporarily negated to get
2090 * the actually hashed constraint.
2092 static void set_max_constant_term(struct sh_data *data, __isl_keep isl_set *set,
2093 int i, isl_int *ineq, uint32_t c_hash, struct ineq_cmp_data *v)
2095 int j;
2096 isl_ctx *ctx;
2097 struct isl_hash_table_entry *entry;
2099 ctx = isl_set_get_ctx(set);
2100 for (j = i + 1; j < set->n; ++j) {
2101 int neg;
2102 isl_int *ineq_j;
2104 entry = isl_hash_table_find(ctx, data->p[j].table,
2105 c_hash, &has_ineq, v, 0);
2106 if (!entry)
2107 continue;
2109 ineq_j = entry->data;
2110 neg = isl_seq_is_neg(ineq_j + 1, ineq + 1, v->len);
2111 if (neg)
2112 isl_int_neg(ineq_j[0], ineq_j[0]);
2113 if (isl_int_gt(ineq_j[0], ineq[0]))
2114 isl_int_set(ineq[0], ineq_j[0]);
2115 if (neg)
2116 isl_int_neg(ineq_j[0], ineq_j[0]);
2120 /* Check if inequality "ineq" from basic set "i" is or can be relaxed to
2121 * become a bound on the whole set. If so, add the (relaxed) inequality
2122 * to "hull". Relaxation is only allowed if "shift" is set.
2124 * We first check if "hull" already contains a translate of the inequality.
2125 * If so, we are done.
2126 * Then, we check if any of the previous basic sets contains a translate
2127 * of the inequality. If so, then we have already considered this
2128 * inequality and we are done.
2129 * Otherwise, for each basic set other than "i", we check if the inequality
2130 * is a bound on the basic set, but first replace the constant term
2131 * by the maximal value of any translate of the inequality in any
2132 * of the following basic sets.
2133 * For previous basic sets, we know that they do not contain a translate
2134 * of the inequality, so we directly call is_bound.
2135 * For following basic sets, we first check if a translate of the
2136 * inequality appears in its description. If so, the constant term
2137 * of the inequality has already been updated with respect to this
2138 * translate and the inequality is therefore known to be a bound
2139 * of this basic set.
2141 static __isl_give isl_basic_set *add_bound(__isl_take isl_basic_set *hull,
2142 struct sh_data *data, __isl_keep isl_set *set, int i, isl_int *ineq,
2143 int shift)
2145 uint32_t c_hash;
2146 struct ineq_cmp_data v;
2147 struct isl_hash_table_entry *entry;
2148 int j, k;
2150 if (!hull)
2151 return NULL;
2153 v.len = isl_basic_set_total_dim(hull);
2154 v.p = ineq;
2155 c_hash = isl_seq_get_hash(ineq + 1, v.len);
2157 entry = isl_hash_table_find(hull->ctx, data->hull_table, c_hash,
2158 has_ineq, &v, 0);
2159 if (entry)
2160 return hull;
2162 for (j = 0; j < i; ++j) {
2163 entry = isl_hash_table_find(hull->ctx, data->p[j].table,
2164 c_hash, has_ineq, &v, 0);
2165 if (entry)
2166 break;
2168 if (j < i)
2169 return hull;
2171 k = isl_basic_set_alloc_inequality(hull);
2172 if (k < 0)
2173 goto error;
2174 isl_seq_cpy(hull->ineq[k], ineq, 1 + v.len);
2176 set_max_constant_term(data, set, i, hull->ineq[k], c_hash, &v);
2177 for (j = 0; j < i; ++j) {
2178 int bound;
2179 bound = is_bound(data, set, j, hull->ineq[k], shift);
2180 if (bound < 0)
2181 goto error;
2182 if (!bound)
2183 break;
2185 if (j < i) {
2186 isl_basic_set_free_inequality(hull, 1);
2187 return hull;
2190 for (j = i + 1; j < set->n; ++j) {
2191 int bound;
2192 entry = isl_hash_table_find(hull->ctx, data->p[j].table,
2193 c_hash, has_ineq, &v, 0);
2194 if (entry)
2195 continue;
2196 bound = is_bound(data, set, j, hull->ineq[k], shift);
2197 if (bound < 0)
2198 goto error;
2199 if (!bound)
2200 break;
2202 if (j < set->n) {
2203 isl_basic_set_free_inequality(hull, 1);
2204 return hull;
2207 entry = isl_hash_table_find(hull->ctx, data->hull_table, c_hash,
2208 has_ineq, &v, 1);
2209 if (!entry)
2210 goto error;
2211 entry->data = hull->ineq[k];
2213 return hull;
2214 error:
2215 isl_basic_set_free(hull);
2216 return NULL;
2219 /* Check if any inequality from basic set "i" is or can be relaxed to
2220 * become a bound on the whole set. If so, add the (relaxed) inequality
2221 * to "hull". Relaxation is only allowed if "shift" is set.
2223 static __isl_give isl_basic_set *add_bounds(__isl_take isl_basic_set *bset,
2224 struct sh_data *data, __isl_keep isl_set *set, int i, int shift)
2226 int j, k;
2227 unsigned dim = isl_basic_set_total_dim(bset);
2229 for (j = 0; j < set->p[i]->n_eq; ++j) {
2230 for (k = 0; k < 2; ++k) {
2231 isl_seq_neg(set->p[i]->eq[j], set->p[i]->eq[j], 1+dim);
2232 bset = add_bound(bset, data, set, i, set->p[i]->eq[j],
2233 shift);
2236 for (j = 0; j < set->p[i]->n_ineq; ++j)
2237 bset = add_bound(bset, data, set, i, set->p[i]->ineq[j], shift);
2238 return bset;
2241 /* Compute a superset of the convex hull of set that is described
2242 * by only (translates of) the constraints in the constituents of set.
2243 * Translation is only allowed if "shift" is set.
2245 static __isl_give isl_basic_set *uset_simple_hull(__isl_take isl_set *set,
2246 int shift)
2248 struct sh_data *data = NULL;
2249 struct isl_basic_set *hull = NULL;
2250 unsigned n_ineq;
2251 int i;
2253 if (!set)
2254 return NULL;
2256 n_ineq = 0;
2257 for (i = 0; i < set->n; ++i) {
2258 if (!set->p[i])
2259 goto error;
2260 n_ineq += 2 * set->p[i]->n_eq + set->p[i]->n_ineq;
2263 hull = isl_basic_set_alloc_space(isl_space_copy(set->dim), 0, 0, n_ineq);
2264 if (!hull)
2265 goto error;
2267 data = sh_data_alloc(set, n_ineq);
2268 if (!data)
2269 goto error;
2271 for (i = 0; i < set->n; ++i)
2272 hull = add_bounds(hull, data, set, i, shift);
2274 sh_data_free(data);
2275 isl_set_free(set);
2277 return hull;
2278 error:
2279 sh_data_free(data);
2280 isl_basic_set_free(hull);
2281 isl_set_free(set);
2282 return NULL;
2285 /* Compute a superset of the convex hull of map that is described
2286 * by only (translates of) the constraints in the constituents of map.
2287 * Handle trivial cases where map is NULL or contains at most one disjunct.
2289 static __isl_give isl_basic_map *map_simple_hull_trivial(
2290 __isl_take isl_map *map)
2292 isl_basic_map *hull;
2294 if (!map)
2295 return NULL;
2296 if (map->n == 0)
2297 return replace_map_by_empty_basic_map(map);
2299 hull = isl_basic_map_copy(map->p[0]);
2300 isl_map_free(map);
2301 return hull;
2304 /* Return a copy of the simple hull cached inside "map".
2305 * "shift" determines whether to return the cached unshifted or shifted
2306 * simple hull.
2308 static __isl_give isl_basic_map *cached_simple_hull(__isl_take isl_map *map,
2309 int shift)
2311 isl_basic_map *hull;
2313 hull = isl_basic_map_copy(map->cached_simple_hull[shift]);
2314 isl_map_free(map);
2316 return hull;
2319 /* Compute a superset of the convex hull of map that is described
2320 * by only (translates of) the constraints in the constituents of map.
2321 * Translation is only allowed if "shift" is set.
2323 * The constraints are sorted while removing redundant constraints
2324 * in order to indicate a preference of which constraints should
2325 * be preserved. In particular, pairs of constraints that are
2326 * sorted together are preferred to either both be preserved
2327 * or both be removed. The sorting is performed inside
2328 * isl_basic_map_remove_redundancies.
2330 * The result of the computation is stored in map->cached_simple_hull[shift]
2331 * such that it can be reused in subsequent calls. The cache is cleared
2332 * whenever the map is modified (in isl_map_cow).
2333 * Note that the results need to be stored in the input map for there
2334 * to be any chance that they may get reused. In particular, they
2335 * are stored in a copy of the input map that is saved before
2336 * the integer division alignment.
2338 static __isl_give isl_basic_map *map_simple_hull(__isl_take isl_map *map,
2339 int shift)
2341 struct isl_set *set = NULL;
2342 struct isl_basic_map *model = NULL;
2343 struct isl_basic_map *hull;
2344 struct isl_basic_map *affine_hull;
2345 struct isl_basic_set *bset = NULL;
2346 isl_map *input;
2348 if (!map || map->n <= 1)
2349 return map_simple_hull_trivial(map);
2351 if (map->cached_simple_hull[shift])
2352 return cached_simple_hull(map, shift);
2354 map = isl_map_detect_equalities(map);
2355 if (!map || map->n <= 1)
2356 return map_simple_hull_trivial(map);
2357 affine_hull = isl_map_affine_hull(isl_map_copy(map));
2358 input = isl_map_copy(map);
2359 map = isl_map_align_divs_internal(map);
2360 model = map ? isl_basic_map_copy(map->p[0]) : NULL;
2362 set = isl_map_underlying_set(map);
2364 bset = uset_simple_hull(set, shift);
2366 hull = isl_basic_map_overlying_set(bset, model);
2368 hull = isl_basic_map_intersect(hull, affine_hull);
2369 hull = isl_basic_map_remove_redundancies(hull);
2371 if (hull) {
2372 ISL_F_SET(hull, ISL_BASIC_MAP_NO_IMPLICIT);
2373 ISL_F_SET(hull, ISL_BASIC_MAP_ALL_EQUALITIES);
2376 hull = isl_basic_map_finalize(hull);
2377 if (input)
2378 input->cached_simple_hull[shift] = isl_basic_map_copy(hull);
2379 isl_map_free(input);
2381 return hull;
2384 /* Compute a superset of the convex hull of map that is described
2385 * by only translates of the constraints in the constituents of map.
2387 __isl_give isl_basic_map *isl_map_simple_hull(__isl_take isl_map *map)
2389 return map_simple_hull(map, 1);
2392 struct isl_basic_set *isl_set_simple_hull(struct isl_set *set)
2394 return bset_from_bmap(isl_map_simple_hull(set_to_map(set)));
2397 /* Compute a superset of the convex hull of map that is described
2398 * by only the constraints in the constituents of map.
2400 __isl_give isl_basic_map *isl_map_unshifted_simple_hull(
2401 __isl_take isl_map *map)
2403 return map_simple_hull(map, 0);
2406 __isl_give isl_basic_set *isl_set_unshifted_simple_hull(
2407 __isl_take isl_set *set)
2409 return isl_map_unshifted_simple_hull(set);
2412 /* Drop all inequalities from "bmap1" that do not also appear in "bmap2".
2413 * A constraint that appears with different constant terms
2414 * in "bmap1" and "bmap2" is also kept, with the least restrictive
2415 * (i.e., greatest) constant term.
2416 * "bmap1" and "bmap2" are assumed to have the same (known)
2417 * integer divisions.
2418 * The constraints of both "bmap1" and "bmap2" are assumed
2419 * to have been sorted using isl_basic_map_sort_constraints.
2421 * Run through the inequality constraints of "bmap1" and "bmap2"
2422 * in sorted order.
2423 * Each constraint of "bmap1" without a matching constraint in "bmap2"
2424 * is removed.
2425 * If a match is found, the constraint is kept. If needed, the constant
2426 * term of the constraint is adjusted.
2428 static __isl_give isl_basic_map *select_shared_inequalities(
2429 __isl_take isl_basic_map *bmap1, __isl_keep isl_basic_map *bmap2)
2431 int i1, i2;
2433 bmap1 = isl_basic_map_cow(bmap1);
2434 if (!bmap1 || !bmap2)
2435 return isl_basic_map_free(bmap1);
2437 i1 = bmap1->n_ineq - 1;
2438 i2 = bmap2->n_ineq - 1;
2439 while (bmap1 && i1 >= 0 && i2 >= 0) {
2440 int cmp;
2442 cmp = isl_basic_map_constraint_cmp(bmap1, bmap1->ineq[i1],
2443 bmap2->ineq[i2]);
2444 if (cmp < 0) {
2445 --i2;
2446 continue;
2448 if (cmp > 0) {
2449 if (isl_basic_map_drop_inequality(bmap1, i1) < 0)
2450 bmap1 = isl_basic_map_free(bmap1);
2451 --i1;
2452 continue;
2454 if (isl_int_lt(bmap1->ineq[i1][0], bmap2->ineq[i2][0]))
2455 isl_int_set(bmap1->ineq[i1][0], bmap2->ineq[i2][0]);
2456 --i1;
2457 --i2;
2459 for (; i1 >= 0; --i1)
2460 if (isl_basic_map_drop_inequality(bmap1, i1) < 0)
2461 bmap1 = isl_basic_map_free(bmap1);
2463 return bmap1;
2466 /* Drop all equalities from "bmap1" that do not also appear in "bmap2".
2467 * "bmap1" and "bmap2" are assumed to have the same (known)
2468 * integer divisions.
2470 * Run through the equality constraints of "bmap1" and "bmap2".
2471 * Each constraint of "bmap1" without a matching constraint in "bmap2"
2472 * is removed.
2474 static __isl_give isl_basic_map *select_shared_equalities(
2475 __isl_take isl_basic_map *bmap1, __isl_keep isl_basic_map *bmap2)
2477 int i1, i2;
2478 unsigned total;
2480 bmap1 = isl_basic_map_cow(bmap1);
2481 if (!bmap1 || !bmap2)
2482 return isl_basic_map_free(bmap1);
2484 total = isl_basic_map_total_dim(bmap1);
2486 i1 = bmap1->n_eq - 1;
2487 i2 = bmap2->n_eq - 1;
2488 while (bmap1 && i1 >= 0 && i2 >= 0) {
2489 int last1, last2;
2491 last1 = isl_seq_last_non_zero(bmap1->eq[i1] + 1, total);
2492 last2 = isl_seq_last_non_zero(bmap2->eq[i2] + 1, total);
2493 if (last1 > last2) {
2494 --i2;
2495 continue;
2497 if (last1 < last2) {
2498 if (isl_basic_map_drop_equality(bmap1, i1) < 0)
2499 bmap1 = isl_basic_map_free(bmap1);
2500 --i1;
2501 continue;
2503 if (!isl_seq_eq(bmap1->eq[i1], bmap2->eq[i2], 1 + total)) {
2504 if (isl_basic_map_drop_equality(bmap1, i1) < 0)
2505 bmap1 = isl_basic_map_free(bmap1);
2507 --i1;
2508 --i2;
2510 for (; i1 >= 0; --i1)
2511 if (isl_basic_map_drop_equality(bmap1, i1) < 0)
2512 bmap1 = isl_basic_map_free(bmap1);
2514 return bmap1;
2517 /* Compute a superset of "bmap1" and "bmap2" that is described
2518 * by only the constraints that appear in both "bmap1" and "bmap2".
2520 * First drop constraints that involve unknown integer divisions
2521 * since it is not trivial to check whether two such integer divisions
2522 * in different basic maps are the same.
2523 * Then align the remaining (known) divs and sort the constraints.
2524 * Finally drop all inequalities and equalities from "bmap1" that
2525 * do not also appear in "bmap2".
2527 __isl_give isl_basic_map *isl_basic_map_plain_unshifted_simple_hull(
2528 __isl_take isl_basic_map *bmap1, __isl_take isl_basic_map *bmap2)
2530 bmap1 = isl_basic_map_drop_constraint_involving_unknown_divs(bmap1);
2531 bmap2 = isl_basic_map_drop_constraint_involving_unknown_divs(bmap2);
2532 bmap2 = isl_basic_map_align_divs(bmap2, bmap1);
2533 bmap1 = isl_basic_map_align_divs(bmap1, bmap2);
2534 bmap1 = isl_basic_map_gauss(bmap1, NULL);
2535 bmap2 = isl_basic_map_gauss(bmap2, NULL);
2536 bmap1 = isl_basic_map_sort_constraints(bmap1);
2537 bmap2 = isl_basic_map_sort_constraints(bmap2);
2539 bmap1 = select_shared_inequalities(bmap1, bmap2);
2540 bmap1 = select_shared_equalities(bmap1, bmap2);
2542 isl_basic_map_free(bmap2);
2543 bmap1 = isl_basic_map_finalize(bmap1);
2544 return bmap1;
2547 /* Compute a superset of the convex hull of "map" that is described
2548 * by only the constraints in the constituents of "map".
2549 * In particular, the result is composed of constraints that appear
2550 * in each of the basic maps of "map"
2552 * Constraints that involve unknown integer divisions are dropped
2553 * since it is not trivial to check whether two such integer divisions
2554 * in different basic maps are the same.
2556 * The hull is initialized from the first basic map and then
2557 * updated with respect to the other basic maps in turn.
2559 __isl_give isl_basic_map *isl_map_plain_unshifted_simple_hull(
2560 __isl_take isl_map *map)
2562 int i;
2563 isl_basic_map *hull;
2565 if (!map)
2566 return NULL;
2567 if (map->n <= 1)
2568 return map_simple_hull_trivial(map);
2569 map = isl_map_drop_constraint_involving_unknown_divs(map);
2570 hull = isl_basic_map_copy(map->p[0]);
2571 for (i = 1; i < map->n; ++i) {
2572 isl_basic_map *bmap_i;
2574 bmap_i = isl_basic_map_copy(map->p[i]);
2575 hull = isl_basic_map_plain_unshifted_simple_hull(hull, bmap_i);
2578 isl_map_free(map);
2579 return hull;
2582 /* Compute a superset of the convex hull of "set" that is described
2583 * by only the constraints in the constituents of "set".
2584 * In particular, the result is composed of constraints that appear
2585 * in each of the basic sets of "set"
2587 __isl_give isl_basic_set *isl_set_plain_unshifted_simple_hull(
2588 __isl_take isl_set *set)
2590 return isl_map_plain_unshifted_simple_hull(set);
2593 /* Check if "ineq" is a bound on "set" and, if so, add it to "hull".
2595 * For each basic set in "set", we first check if the basic set
2596 * contains a translate of "ineq". If this translate is more relaxed,
2597 * then we assume that "ineq" is not a bound on this basic set.
2598 * Otherwise, we know that it is a bound.
2599 * If the basic set does not contain a translate of "ineq", then
2600 * we call is_bound to perform the test.
2602 static __isl_give isl_basic_set *add_bound_from_constraint(
2603 __isl_take isl_basic_set *hull, struct sh_data *data,
2604 __isl_keep isl_set *set, isl_int *ineq)
2606 int i, k;
2607 isl_ctx *ctx;
2608 uint32_t c_hash;
2609 struct ineq_cmp_data v;
2611 if (!hull || !set)
2612 return isl_basic_set_free(hull);
2614 v.len = isl_basic_set_total_dim(hull);
2615 v.p = ineq;
2616 c_hash = isl_seq_get_hash(ineq + 1, v.len);
2618 ctx = isl_basic_set_get_ctx(hull);
2619 for (i = 0; i < set->n; ++i) {
2620 int bound;
2621 struct isl_hash_table_entry *entry;
2623 entry = isl_hash_table_find(ctx, data->p[i].table,
2624 c_hash, &has_ineq, &v, 0);
2625 if (entry) {
2626 isl_int *ineq_i = entry->data;
2627 int neg, more_relaxed;
2629 neg = isl_seq_is_neg(ineq_i + 1, ineq + 1, v.len);
2630 if (neg)
2631 isl_int_neg(ineq_i[0], ineq_i[0]);
2632 more_relaxed = isl_int_gt(ineq_i[0], ineq[0]);
2633 if (neg)
2634 isl_int_neg(ineq_i[0], ineq_i[0]);
2635 if (more_relaxed)
2636 break;
2637 else
2638 continue;
2640 bound = is_bound(data, set, i, ineq, 0);
2641 if (bound < 0)
2642 return isl_basic_set_free(hull);
2643 if (!bound)
2644 break;
2646 if (i < set->n)
2647 return hull;
2649 k = isl_basic_set_alloc_inequality(hull);
2650 if (k < 0)
2651 return isl_basic_set_free(hull);
2652 isl_seq_cpy(hull->ineq[k], ineq, 1 + v.len);
2654 return hull;
2657 /* Compute a superset of the convex hull of "set" that is described
2658 * by only some of the "n_ineq" constraints in the list "ineq", where "set"
2659 * has no parameters or integer divisions.
2661 * The inequalities in "ineq" are assumed to have been sorted such
2662 * that constraints with the same linear part appear together and
2663 * that among constraints with the same linear part, those with
2664 * smaller constant term appear first.
2666 * We reuse the same data structure that is used by uset_simple_hull,
2667 * but we do not need the hull table since we will not consider the
2668 * same constraint more than once. We therefore allocate it with zero size.
2670 * We run through the constraints and try to add them one by one,
2671 * skipping identical constraints. If we have added a constraint and
2672 * the next constraint is a more relaxed translate, then we skip this
2673 * next constraint as well.
2675 static __isl_give isl_basic_set *uset_unshifted_simple_hull_from_constraints(
2676 __isl_take isl_set *set, int n_ineq, isl_int **ineq)
2678 int i;
2679 int last_added = 0;
2680 struct sh_data *data = NULL;
2681 isl_basic_set *hull = NULL;
2682 unsigned dim;
2684 hull = isl_basic_set_alloc_space(isl_set_get_space(set), 0, 0, n_ineq);
2685 if (!hull)
2686 goto error;
2688 data = sh_data_alloc(set, 0);
2689 if (!data)
2690 goto error;
2692 dim = isl_set_dim(set, isl_dim_set);
2693 for (i = 0; i < n_ineq; ++i) {
2694 int hull_n_ineq = hull->n_ineq;
2695 int parallel;
2697 parallel = i > 0 && isl_seq_eq(ineq[i - 1] + 1, ineq[i] + 1,
2698 dim);
2699 if (parallel &&
2700 (last_added || isl_int_eq(ineq[i - 1][0], ineq[i][0])))
2701 continue;
2702 hull = add_bound_from_constraint(hull, data, set, ineq[i]);
2703 if (!hull)
2704 goto error;
2705 last_added = hull->n_ineq > hull_n_ineq;
2708 sh_data_free(data);
2709 isl_set_free(set);
2710 return hull;
2711 error:
2712 sh_data_free(data);
2713 isl_set_free(set);
2714 isl_basic_set_free(hull);
2715 return NULL;
2718 /* Collect pointers to all the inequalities in the elements of "list"
2719 * in "ineq". For equalities, store both a pointer to the equality and
2720 * a pointer to its opposite, which is first copied to "mat".
2721 * "ineq" and "mat" are assumed to have been preallocated to the right size
2722 * (the number of inequalities + 2 times the number of equalites and
2723 * the number of equalities, respectively).
2725 static __isl_give isl_mat *collect_inequalities(__isl_take isl_mat *mat,
2726 __isl_keep isl_basic_set_list *list, isl_int **ineq)
2728 int i, j, n, n_eq, n_ineq;
2730 if (!mat)
2731 return NULL;
2733 n_eq = 0;
2734 n_ineq = 0;
2735 n = isl_basic_set_list_n_basic_set(list);
2736 for (i = 0; i < n; ++i) {
2737 isl_basic_set *bset;
2738 bset = isl_basic_set_list_get_basic_set(list, i);
2739 if (!bset)
2740 return isl_mat_free(mat);
2741 for (j = 0; j < bset->n_eq; ++j) {
2742 ineq[n_ineq++] = mat->row[n_eq];
2743 ineq[n_ineq++] = bset->eq[j];
2744 isl_seq_neg(mat->row[n_eq++], bset->eq[j], mat->n_col);
2746 for (j = 0; j < bset->n_ineq; ++j)
2747 ineq[n_ineq++] = bset->ineq[j];
2748 isl_basic_set_free(bset);
2751 return mat;
2754 /* Comparison routine for use as an isl_sort callback.
2756 * Constraints with the same linear part are sorted together and
2757 * among constraints with the same linear part, those with smaller
2758 * constant term are sorted first.
2760 static int cmp_ineq(const void *a, const void *b, void *arg)
2762 unsigned dim = *(unsigned *) arg;
2763 isl_int * const *ineq1 = a;
2764 isl_int * const *ineq2 = b;
2765 int cmp;
2767 cmp = isl_seq_cmp((*ineq1) + 1, (*ineq2) + 1, dim);
2768 if (cmp != 0)
2769 return cmp;
2770 return isl_int_cmp((*ineq1)[0], (*ineq2)[0]);
2773 /* Compute a superset of the convex hull of "set" that is described
2774 * by only constraints in the elements of "list", where "set" has
2775 * no parameters or integer divisions.
2777 * We collect all the constraints in those elements and then
2778 * sort the constraints such that constraints with the same linear part
2779 * are sorted together and that those with smaller constant term are
2780 * sorted first.
2782 static __isl_give isl_basic_set *uset_unshifted_simple_hull_from_basic_set_list(
2783 __isl_take isl_set *set, __isl_take isl_basic_set_list *list)
2785 int i, n, n_eq, n_ineq;
2786 unsigned dim;
2787 isl_ctx *ctx;
2788 isl_mat *mat = NULL;
2789 isl_int **ineq = NULL;
2790 isl_basic_set *hull;
2792 if (!set)
2793 goto error;
2794 ctx = isl_set_get_ctx(set);
2796 n_eq = 0;
2797 n_ineq = 0;
2798 n = isl_basic_set_list_n_basic_set(list);
2799 for (i = 0; i < n; ++i) {
2800 isl_basic_set *bset;
2801 bset = isl_basic_set_list_get_basic_set(list, i);
2802 if (!bset)
2803 goto error;
2804 n_eq += bset->n_eq;
2805 n_ineq += 2 * bset->n_eq + bset->n_ineq;
2806 isl_basic_set_free(bset);
2809 ineq = isl_alloc_array(ctx, isl_int *, n_ineq);
2810 if (n_ineq > 0 && !ineq)
2811 goto error;
2813 dim = isl_set_dim(set, isl_dim_set);
2814 mat = isl_mat_alloc(ctx, n_eq, 1 + dim);
2815 mat = collect_inequalities(mat, list, ineq);
2816 if (!mat)
2817 goto error;
2819 if (isl_sort(ineq, n_ineq, sizeof(ineq[0]), &cmp_ineq, &dim) < 0)
2820 goto error;
2822 hull = uset_unshifted_simple_hull_from_constraints(set, n_ineq, ineq);
2824 isl_mat_free(mat);
2825 free(ineq);
2826 isl_basic_set_list_free(list);
2827 return hull;
2828 error:
2829 isl_mat_free(mat);
2830 free(ineq);
2831 isl_set_free(set);
2832 isl_basic_set_list_free(list);
2833 return NULL;
2836 /* Compute a superset of the convex hull of "map" that is described
2837 * by only constraints in the elements of "list".
2839 * If the list is empty, then we can only describe the universe set.
2840 * If the input map is empty, then all constraints are valid, so
2841 * we return the intersection of the elements in "list".
2843 * Otherwise, we align all divs and temporarily treat them
2844 * as regular variables, computing the unshifted simple hull in
2845 * uset_unshifted_simple_hull_from_basic_set_list.
2847 static __isl_give isl_basic_map *map_unshifted_simple_hull_from_basic_map_list(
2848 __isl_take isl_map *map, __isl_take isl_basic_map_list *list)
2850 isl_basic_map *model;
2851 isl_basic_map *hull;
2852 isl_set *set;
2853 isl_basic_set_list *bset_list;
2855 if (!map || !list)
2856 goto error;
2858 if (isl_basic_map_list_n_basic_map(list) == 0) {
2859 isl_space *space;
2861 space = isl_map_get_space(map);
2862 isl_map_free(map);
2863 isl_basic_map_list_free(list);
2864 return isl_basic_map_universe(space);
2866 if (isl_map_plain_is_empty(map)) {
2867 isl_map_free(map);
2868 return isl_basic_map_list_intersect(list);
2871 map = isl_map_align_divs_to_basic_map_list(map, list);
2872 if (!map)
2873 goto error;
2874 list = isl_basic_map_list_align_divs_to_basic_map(list, map->p[0]);
2876 model = isl_basic_map_list_get_basic_map(list, 0);
2878 set = isl_map_underlying_set(map);
2879 bset_list = isl_basic_map_list_underlying_set(list);
2881 hull = uset_unshifted_simple_hull_from_basic_set_list(set, bset_list);
2882 hull = isl_basic_map_overlying_set(hull, model);
2884 return hull;
2885 error:
2886 isl_map_free(map);
2887 isl_basic_map_list_free(list);
2888 return NULL;
2891 /* Return a sequence of the basic maps that make up the maps in "list".
2893 static __isl_give isl_basic_map_list *collect_basic_maps(
2894 __isl_take isl_map_list *list)
2896 int i, n;
2897 isl_ctx *ctx;
2898 isl_basic_map_list *bmap_list;
2900 if (!list)
2901 return NULL;
2902 n = isl_map_list_n_map(list);
2903 ctx = isl_map_list_get_ctx(list);
2904 bmap_list = isl_basic_map_list_alloc(ctx, 0);
2906 for (i = 0; i < n; ++i) {
2907 isl_map *map;
2908 isl_basic_map_list *list_i;
2910 map = isl_map_list_get_map(list, i);
2911 map = isl_map_compute_divs(map);
2912 list_i = isl_map_get_basic_map_list(map);
2913 isl_map_free(map);
2914 bmap_list = isl_basic_map_list_concat(bmap_list, list_i);
2917 isl_map_list_free(list);
2918 return bmap_list;
2921 /* Compute a superset of the convex hull of "map" that is described
2922 * by only constraints in the elements of "list".
2924 * If "map" is the universe, then the convex hull (and therefore
2925 * any superset of the convexhull) is the universe as well.
2927 * Otherwise, we collect all the basic maps in the map list and
2928 * continue with map_unshifted_simple_hull_from_basic_map_list.
2930 __isl_give isl_basic_map *isl_map_unshifted_simple_hull_from_map_list(
2931 __isl_take isl_map *map, __isl_take isl_map_list *list)
2933 isl_basic_map_list *bmap_list;
2934 int is_universe;
2936 is_universe = isl_map_plain_is_universe(map);
2937 if (is_universe < 0)
2938 map = isl_map_free(map);
2939 if (is_universe < 0 || is_universe) {
2940 isl_map_list_free(list);
2941 return isl_map_unshifted_simple_hull(map);
2944 bmap_list = collect_basic_maps(list);
2945 return map_unshifted_simple_hull_from_basic_map_list(map, bmap_list);
2948 /* Compute a superset of the convex hull of "set" that is described
2949 * by only constraints in the elements of "list".
2951 __isl_give isl_basic_set *isl_set_unshifted_simple_hull_from_set_list(
2952 __isl_take isl_set *set, __isl_take isl_set_list *list)
2954 return isl_map_unshifted_simple_hull_from_map_list(set, list);
2957 /* Given a set "set", return parametric bounds on the dimension "dim".
2959 static __isl_give isl_basic_set *set_bounds(__isl_keep isl_set *set, int dim)
2961 unsigned set_dim = isl_set_dim(set, isl_dim_set);
2962 set = isl_set_copy(set);
2963 set = isl_set_eliminate_dims(set, dim + 1, set_dim - (dim + 1));
2964 set = isl_set_eliminate_dims(set, 0, dim);
2965 return isl_set_convex_hull(set);
2968 /* Computes a "simple hull" and then check if each dimension in the
2969 * resulting hull is bounded by a symbolic constant. If not, the
2970 * hull is intersected with the corresponding bounds on the whole set.
2972 __isl_give isl_basic_set *isl_set_bounded_simple_hull(__isl_take isl_set *set)
2974 int i, j;
2975 struct isl_basic_set *hull;
2976 unsigned nparam, left;
2977 int removed_divs = 0;
2979 hull = isl_set_simple_hull(isl_set_copy(set));
2980 if (!hull)
2981 goto error;
2983 nparam = isl_basic_set_dim(hull, isl_dim_param);
2984 for (i = 0; i < isl_basic_set_dim(hull, isl_dim_set); ++i) {
2985 int lower = 0, upper = 0;
2986 struct isl_basic_set *bounds;
2988 left = isl_basic_set_total_dim(hull) - nparam - i - 1;
2989 for (j = 0; j < hull->n_eq; ++j) {
2990 if (isl_int_is_zero(hull->eq[j][1 + nparam + i]))
2991 continue;
2992 if (isl_seq_first_non_zero(hull->eq[j]+1+nparam+i+1,
2993 left) == -1)
2994 break;
2996 if (j < hull->n_eq)
2997 continue;
2999 for (j = 0; j < hull->n_ineq; ++j) {
3000 if (isl_int_is_zero(hull->ineq[j][1 + nparam + i]))
3001 continue;
3002 if (isl_seq_first_non_zero(hull->ineq[j]+1+nparam+i+1,
3003 left) != -1 ||
3004 isl_seq_first_non_zero(hull->ineq[j]+1+nparam,
3005 i) != -1)
3006 continue;
3007 if (isl_int_is_pos(hull->ineq[j][1 + nparam + i]))
3008 lower = 1;
3009 else
3010 upper = 1;
3011 if (lower && upper)
3012 break;
3015 if (lower && upper)
3016 continue;
3018 if (!removed_divs) {
3019 set = isl_set_remove_divs(set);
3020 if (!set)
3021 goto error;
3022 removed_divs = 1;
3024 bounds = set_bounds(set, i);
3025 hull = isl_basic_set_intersect(hull, bounds);
3026 if (!hull)
3027 goto error;
3030 isl_set_free(set);
3031 return hull;
3032 error:
3033 isl_set_free(set);
3034 isl_basic_set_free(hull);
3035 return NULL;