isl_ast_codegen.c: create_node_scaled: inline call to isl_ast_graft_alloc_level
[isl.git] / isl_convex_hull.c
blobb160ec9279558217dd7a610724635969e243ee64
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 static struct isl_basic_set *uset_convex_hull_wrap_bounded(struct isl_set *set);
28 /* Return 1 if constraint c is redundant with respect to the constraints
29 * in bmap. If c is a lower [upper] bound in some variable and bmap
30 * does not have a lower [upper] bound in that variable, then c cannot
31 * be redundant and we do not need solve any lp.
33 int isl_basic_map_constraint_is_redundant(struct isl_basic_map **bmap,
34 isl_int *c, isl_int *opt_n, isl_int *opt_d)
36 enum isl_lp_result res;
37 unsigned total;
38 int i, j;
40 if (!bmap)
41 return -1;
43 total = isl_basic_map_total_dim(*bmap);
44 for (i = 0; i < total; ++i) {
45 int sign;
46 if (isl_int_is_zero(c[1+i]))
47 continue;
48 sign = isl_int_sgn(c[1+i]);
49 for (j = 0; j < (*bmap)->n_ineq; ++j)
50 if (sign == isl_int_sgn((*bmap)->ineq[j][1+i]))
51 break;
52 if (j == (*bmap)->n_ineq)
53 break;
55 if (i < total)
56 return 0;
58 res = isl_basic_map_solve_lp(*bmap, 0, c, (*bmap)->ctx->one,
59 opt_n, opt_d, NULL);
60 if (res == isl_lp_unbounded)
61 return 0;
62 if (res == isl_lp_error)
63 return -1;
64 if (res == isl_lp_empty) {
65 *bmap = isl_basic_map_set_to_empty(*bmap);
66 return 0;
68 return !isl_int_is_neg(*opt_n);
71 int isl_basic_set_constraint_is_redundant(struct isl_basic_set **bset,
72 isl_int *c, isl_int *opt_n, isl_int *opt_d)
74 return isl_basic_map_constraint_is_redundant(
75 (struct isl_basic_map **)bset, c, opt_n, opt_d);
78 /* Remove redundant
79 * constraints. If the minimal value along the normal of a constraint
80 * is the same if the constraint is removed, then the constraint is redundant.
82 * Alternatively, we could have intersected the basic map with the
83 * corresponding equality and the checked if the dimension was that
84 * of a facet.
86 __isl_give isl_basic_map *isl_basic_map_remove_redundancies(
87 __isl_take isl_basic_map *bmap)
89 struct isl_tab *tab;
91 if (!bmap)
92 return NULL;
94 bmap = isl_basic_map_gauss(bmap, NULL);
95 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
96 return bmap;
97 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_NO_REDUNDANT))
98 return bmap;
99 if (bmap->n_ineq <= 1)
100 return bmap;
102 tab = isl_tab_from_basic_map(bmap, 0);
103 if (isl_tab_detect_implicit_equalities(tab) < 0)
104 goto error;
105 if (isl_tab_detect_redundant(tab) < 0)
106 goto error;
107 bmap = isl_basic_map_update_from_tab(bmap, tab);
108 isl_tab_free(tab);
109 ISL_F_SET(bmap, ISL_BASIC_MAP_NO_IMPLICIT);
110 ISL_F_SET(bmap, ISL_BASIC_MAP_NO_REDUNDANT);
111 return bmap;
112 error:
113 isl_tab_free(tab);
114 isl_basic_map_free(bmap);
115 return NULL;
118 __isl_give isl_basic_set *isl_basic_set_remove_redundancies(
119 __isl_take isl_basic_set *bset)
121 return (struct isl_basic_set *)
122 isl_basic_map_remove_redundancies((struct isl_basic_map *)bset);
125 /* Remove redundant constraints in each of the basic maps.
127 __isl_give isl_map *isl_map_remove_redundancies(__isl_take isl_map *map)
129 return isl_map_inline_foreach_basic_map(map,
130 &isl_basic_map_remove_redundancies);
133 __isl_give isl_set *isl_set_remove_redundancies(__isl_take isl_set *set)
135 return isl_map_remove_redundancies(set);
138 /* Check if the set set is bound in the direction of the affine
139 * constraint c and if so, set the constant term such that the
140 * resulting constraint is a bounding constraint for the set.
142 static int uset_is_bound(struct isl_set *set, isl_int *c, unsigned len)
144 int first;
145 int j;
146 isl_int opt;
147 isl_int opt_denom;
149 isl_int_init(opt);
150 isl_int_init(opt_denom);
151 first = 1;
152 for (j = 0; j < set->n; ++j) {
153 enum isl_lp_result res;
155 if (ISL_F_ISSET(set->p[j], ISL_BASIC_SET_EMPTY))
156 continue;
158 res = isl_basic_set_solve_lp(set->p[j],
159 0, c, set->ctx->one, &opt, &opt_denom, NULL);
160 if (res == isl_lp_unbounded)
161 break;
162 if (res == isl_lp_error)
163 goto error;
164 if (res == isl_lp_empty) {
165 set->p[j] = isl_basic_set_set_to_empty(set->p[j]);
166 if (!set->p[j])
167 goto error;
168 continue;
170 if (first || isl_int_is_neg(opt)) {
171 if (!isl_int_is_one(opt_denom))
172 isl_seq_scale(c, c, opt_denom, len);
173 isl_int_sub(c[0], c[0], opt);
175 first = 0;
177 isl_int_clear(opt);
178 isl_int_clear(opt_denom);
179 return j >= set->n;
180 error:
181 isl_int_clear(opt);
182 isl_int_clear(opt_denom);
183 return -1;
186 __isl_give isl_basic_map *isl_basic_map_set_rational(
187 __isl_take isl_basic_set *bmap)
189 if (!bmap)
190 return NULL;
192 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
193 return bmap;
195 bmap = isl_basic_map_cow(bmap);
196 if (!bmap)
197 return NULL;
199 ISL_F_SET(bmap, ISL_BASIC_MAP_RATIONAL);
201 return isl_basic_map_finalize(bmap);
204 __isl_give isl_basic_set *isl_basic_set_set_rational(
205 __isl_take isl_basic_set *bset)
207 return isl_basic_map_set_rational(bset);
210 __isl_give isl_map *isl_map_set_rational(__isl_take isl_map *map)
212 int i;
214 map = isl_map_cow(map);
215 if (!map)
216 return NULL;
217 for (i = 0; i < map->n; ++i) {
218 map->p[i] = isl_basic_map_set_rational(map->p[i]);
219 if (!map->p[i])
220 goto error;
222 return map;
223 error:
224 isl_map_free(map);
225 return NULL;
228 __isl_give isl_set *isl_set_set_rational(__isl_take isl_set *set)
230 return isl_map_set_rational(set);
233 static struct isl_basic_set *isl_basic_set_add_equality(
234 struct isl_basic_set *bset, isl_int *c)
236 int i;
237 unsigned dim;
239 if (!bset)
240 return NULL;
242 if (ISL_F_ISSET(bset, ISL_BASIC_SET_EMPTY))
243 return bset;
245 isl_assert(bset->ctx, isl_basic_set_n_param(bset) == 0, goto error);
246 isl_assert(bset->ctx, bset->n_div == 0, goto error);
247 dim = isl_basic_set_n_dim(bset);
248 bset = isl_basic_set_cow(bset);
249 bset = isl_basic_set_extend(bset, 0, dim, 0, 1, 0);
250 i = isl_basic_set_alloc_equality(bset);
251 if (i < 0)
252 goto error;
253 isl_seq_cpy(bset->eq[i], c, 1 + dim);
254 return bset;
255 error:
256 isl_basic_set_free(bset);
257 return NULL;
260 static struct isl_set *isl_set_add_basic_set_equality(struct isl_set *set, isl_int *c)
262 int i;
264 set = isl_set_cow(set);
265 if (!set)
266 return NULL;
267 for (i = 0; i < set->n; ++i) {
268 set->p[i] = isl_basic_set_add_equality(set->p[i], c);
269 if (!set->p[i])
270 goto error;
272 return set;
273 error:
274 isl_set_free(set);
275 return NULL;
278 /* Given a union of basic sets, construct the constraints for wrapping
279 * a facet around one of its ridges.
280 * In particular, if each of n the d-dimensional basic sets i in "set"
281 * contains the origin, satisfies the constraints x_1 >= 0 and x_2 >= 0
282 * and is defined by the constraints
283 * [ 1 ]
284 * A_i [ x ] >= 0
286 * then the resulting set is of dimension n*(1+d) and has as constraints
288 * [ a_i ]
289 * A_i [ x_i ] >= 0
291 * a_i >= 0
293 * \sum_i x_{i,1} = 1
295 static struct isl_basic_set *wrap_constraints(struct isl_set *set)
297 struct isl_basic_set *lp;
298 unsigned n_eq;
299 unsigned n_ineq;
300 int i, j, k;
301 unsigned dim, lp_dim;
303 if (!set)
304 return NULL;
306 dim = 1 + isl_set_n_dim(set);
307 n_eq = 1;
308 n_ineq = set->n;
309 for (i = 0; i < set->n; ++i) {
310 n_eq += set->p[i]->n_eq;
311 n_ineq += set->p[i]->n_ineq;
313 lp = isl_basic_set_alloc(set->ctx, 0, dim * set->n, 0, n_eq, n_ineq);
314 lp = isl_basic_set_set_rational(lp);
315 if (!lp)
316 return NULL;
317 lp_dim = isl_basic_set_n_dim(lp);
318 k = isl_basic_set_alloc_equality(lp);
319 isl_int_set_si(lp->eq[k][0], -1);
320 for (i = 0; i < set->n; ++i) {
321 isl_int_set_si(lp->eq[k][1+dim*i], 0);
322 isl_int_set_si(lp->eq[k][1+dim*i+1], 1);
323 isl_seq_clr(lp->eq[k]+1+dim*i+2, dim-2);
325 for (i = 0; i < set->n; ++i) {
326 k = isl_basic_set_alloc_inequality(lp);
327 isl_seq_clr(lp->ineq[k], 1+lp_dim);
328 isl_int_set_si(lp->ineq[k][1+dim*i], 1);
330 for (j = 0; j < set->p[i]->n_eq; ++j) {
331 k = isl_basic_set_alloc_equality(lp);
332 isl_seq_clr(lp->eq[k], 1+dim*i);
333 isl_seq_cpy(lp->eq[k]+1+dim*i, set->p[i]->eq[j], dim);
334 isl_seq_clr(lp->eq[k]+1+dim*(i+1), dim*(set->n-i-1));
337 for (j = 0; j < set->p[i]->n_ineq; ++j) {
338 k = isl_basic_set_alloc_inequality(lp);
339 isl_seq_clr(lp->ineq[k], 1+dim*i);
340 isl_seq_cpy(lp->ineq[k]+1+dim*i, set->p[i]->ineq[j], dim);
341 isl_seq_clr(lp->ineq[k]+1+dim*(i+1), dim*(set->n-i-1));
344 return lp;
347 /* Given a facet "facet" of the convex hull of "set" and a facet "ridge"
348 * of that facet, compute the other facet of the convex hull that contains
349 * the ridge.
351 * We first transform the set such that the facet constraint becomes
353 * x_1 >= 0
355 * I.e., the facet lies in
357 * x_1 = 0
359 * and on that facet, the constraint that defines the ridge is
361 * x_2 >= 0
363 * (This transformation is not strictly needed, all that is needed is
364 * that the ridge contains the origin.)
366 * Since the ridge contains the origin, the cone of the convex hull
367 * will be of the form
369 * x_1 >= 0
370 * x_2 >= a x_1
372 * with this second constraint defining the new facet.
373 * The constant a is obtained by settting x_1 in the cone of the
374 * convex hull to 1 and minimizing x_2.
375 * Now, each element in the cone of the convex hull is the sum
376 * of elements in the cones of the basic sets.
377 * If a_i is the dilation factor of basic set i, then the problem
378 * we need to solve is
380 * min \sum_i x_{i,2}
381 * st
382 * \sum_i x_{i,1} = 1
383 * a_i >= 0
384 * [ a_i ]
385 * A [ x_i ] >= 0
387 * with
388 * [ 1 ]
389 * A_i [ x_i ] >= 0
391 * the constraints of each (transformed) basic set.
392 * If a = n/d, then the constraint defining the new facet (in the transformed
393 * space) is
395 * -n x_1 + d x_2 >= 0
397 * In the original space, we need to take the same combination of the
398 * corresponding constraints "facet" and "ridge".
400 * If a = -infty = "-1/0", then we just return the original facet constraint.
401 * This means that the facet is unbounded, but has a bounded intersection
402 * with the union of sets.
404 isl_int *isl_set_wrap_facet(__isl_keep isl_set *set,
405 isl_int *facet, isl_int *ridge)
407 int i;
408 isl_ctx *ctx;
409 struct isl_mat *T = NULL;
410 struct isl_basic_set *lp = NULL;
411 struct isl_vec *obj;
412 enum isl_lp_result res;
413 isl_int num, den;
414 unsigned dim;
416 if (!set)
417 return NULL;
418 ctx = set->ctx;
419 set = isl_set_copy(set);
420 set = isl_set_set_rational(set);
422 dim = 1 + isl_set_n_dim(set);
423 T = isl_mat_alloc(ctx, 3, dim);
424 if (!T)
425 goto error;
426 isl_int_set_si(T->row[0][0], 1);
427 isl_seq_clr(T->row[0]+1, dim - 1);
428 isl_seq_cpy(T->row[1], facet, dim);
429 isl_seq_cpy(T->row[2], ridge, dim);
430 T = isl_mat_right_inverse(T);
431 set = isl_set_preimage(set, T);
432 T = NULL;
433 if (!set)
434 goto error;
435 lp = wrap_constraints(set);
436 obj = isl_vec_alloc(ctx, 1 + dim*set->n);
437 if (!obj)
438 goto error;
439 isl_int_set_si(obj->block.data[0], 0);
440 for (i = 0; i < set->n; ++i) {
441 isl_seq_clr(obj->block.data + 1 + dim*i, 2);
442 isl_int_set_si(obj->block.data[1 + dim*i+2], 1);
443 isl_seq_clr(obj->block.data + 1 + dim*i+3, dim-3);
445 isl_int_init(num);
446 isl_int_init(den);
447 res = isl_basic_set_solve_lp(lp, 0,
448 obj->block.data, ctx->one, &num, &den, NULL);
449 if (res == isl_lp_ok) {
450 isl_int_neg(num, num);
451 isl_seq_combine(facet, num, facet, den, ridge, dim);
452 isl_seq_normalize(ctx, facet, dim);
454 isl_int_clear(num);
455 isl_int_clear(den);
456 isl_vec_free(obj);
457 isl_basic_set_free(lp);
458 isl_set_free(set);
459 if (res == isl_lp_error)
460 return NULL;
461 isl_assert(ctx, res == isl_lp_ok || res == isl_lp_unbounded,
462 return NULL);
463 return facet;
464 error:
465 isl_basic_set_free(lp);
466 isl_mat_free(T);
467 isl_set_free(set);
468 return NULL;
471 /* Compute the constraint of a facet of "set".
473 * We first compute the intersection with a bounding constraint
474 * that is orthogonal to one of the coordinate axes.
475 * If the affine hull of this intersection has only one equality,
476 * we have found a facet.
477 * Otherwise, we wrap the current bounding constraint around
478 * one of the equalities of the face (one that is not equal to
479 * the current bounding constraint).
480 * This process continues until we have found a facet.
481 * The dimension of the intersection increases by at least
482 * one on each iteration, so termination is guaranteed.
484 static __isl_give isl_mat *initial_facet_constraint(__isl_keep isl_set *set)
486 struct isl_set *slice = NULL;
487 struct isl_basic_set *face = NULL;
488 int i;
489 unsigned dim = isl_set_n_dim(set);
490 int is_bound;
491 isl_mat *bounds = NULL;
493 isl_assert(set->ctx, set->n > 0, goto error);
494 bounds = isl_mat_alloc(set->ctx, 1, 1 + dim);
495 if (!bounds)
496 return NULL;
498 isl_seq_clr(bounds->row[0], dim);
499 isl_int_set_si(bounds->row[0][1 + dim - 1], 1);
500 is_bound = uset_is_bound(set, bounds->row[0], 1 + dim);
501 if (is_bound < 0)
502 goto error;
503 isl_assert(set->ctx, is_bound, goto error);
504 isl_seq_normalize(set->ctx, bounds->row[0], 1 + dim);
505 bounds->n_row = 1;
507 for (;;) {
508 slice = isl_set_copy(set);
509 slice = isl_set_add_basic_set_equality(slice, bounds->row[0]);
510 face = isl_set_affine_hull(slice);
511 if (!face)
512 goto error;
513 if (face->n_eq == 1) {
514 isl_basic_set_free(face);
515 break;
517 for (i = 0; i < face->n_eq; ++i)
518 if (!isl_seq_eq(bounds->row[0], face->eq[i], 1 + dim) &&
519 !isl_seq_is_neg(bounds->row[0],
520 face->eq[i], 1 + dim))
521 break;
522 isl_assert(set->ctx, i < face->n_eq, goto error);
523 if (!isl_set_wrap_facet(set, bounds->row[0], face->eq[i]))
524 goto error;
525 isl_seq_normalize(set->ctx, bounds->row[0], bounds->n_col);
526 isl_basic_set_free(face);
529 return bounds;
530 error:
531 isl_basic_set_free(face);
532 isl_mat_free(bounds);
533 return NULL;
536 /* Given the bounding constraint "c" of a facet of the convex hull of "set",
537 * compute a hyperplane description of the facet, i.e., compute the facets
538 * of the facet.
540 * We compute an affine transformation that transforms the constraint
542 * [ 1 ]
543 * c [ x ] = 0
545 * to the constraint
547 * z_1 = 0
549 * by computing the right inverse U of a matrix that starts with the rows
551 * [ 1 0 ]
552 * [ c ]
554 * Then
555 * [ 1 ] [ 1 ]
556 * [ x ] = U [ z ]
557 * and
558 * [ 1 ] [ 1 ]
559 * [ z ] = Q [ x ]
561 * with Q = U^{-1}
562 * Since z_1 is zero, we can drop this variable as well as the corresponding
563 * column of U to obtain
565 * [ 1 ] [ 1 ]
566 * [ x ] = U' [ z' ]
567 * and
568 * [ 1 ] [ 1 ]
569 * [ z' ] = Q' [ x ]
571 * with Q' equal to Q, but without the corresponding row.
572 * After computing the facets of the facet in the z' space,
573 * we convert them back to the x space through Q.
575 static struct isl_basic_set *compute_facet(struct isl_set *set, isl_int *c)
577 struct isl_mat *m, *U, *Q;
578 struct isl_basic_set *facet = NULL;
579 struct isl_ctx *ctx;
580 unsigned dim;
582 ctx = set->ctx;
583 set = isl_set_copy(set);
584 dim = isl_set_n_dim(set);
585 m = isl_mat_alloc(set->ctx, 2, 1 + dim);
586 if (!m)
587 goto error;
588 isl_int_set_si(m->row[0][0], 1);
589 isl_seq_clr(m->row[0]+1, dim);
590 isl_seq_cpy(m->row[1], c, 1+dim);
591 U = isl_mat_right_inverse(m);
592 Q = isl_mat_right_inverse(isl_mat_copy(U));
593 U = isl_mat_drop_cols(U, 1, 1);
594 Q = isl_mat_drop_rows(Q, 1, 1);
595 set = isl_set_preimage(set, U);
596 facet = uset_convex_hull_wrap_bounded(set);
597 facet = isl_basic_set_preimage(facet, Q);
598 if (facet)
599 isl_assert(ctx, facet->n_eq == 0, goto error);
600 return facet;
601 error:
602 isl_basic_set_free(facet);
603 isl_set_free(set);
604 return NULL;
607 /* Given an initial facet constraint, compute the remaining facets.
608 * We do this by running through all facets found so far and computing
609 * the adjacent facets through wrapping, adding those facets that we
610 * hadn't already found before.
612 * For each facet we have found so far, we first compute its facets
613 * in the resulting convex hull. That is, we compute the ridges
614 * of the resulting convex hull contained in the facet.
615 * We also compute the corresponding facet in the current approximation
616 * of the convex hull. There is no need to wrap around the ridges
617 * in this facet since that would result in a facet that is already
618 * present in the current approximation.
620 * This function can still be significantly optimized by checking which of
621 * the facets of the basic sets are also facets of the convex hull and
622 * using all the facets so far to help in constructing the facets of the
623 * facets
624 * and/or
625 * using the technique in section "3.1 Ridge Generation" of
626 * "Extended Convex Hull" by Fukuda et al.
628 static struct isl_basic_set *extend(struct isl_basic_set *hull,
629 struct isl_set *set)
631 int i, j, f;
632 int k;
633 struct isl_basic_set *facet = NULL;
634 struct isl_basic_set *hull_facet = NULL;
635 unsigned dim;
637 if (!hull)
638 return NULL;
640 isl_assert(set->ctx, set->n > 0, goto error);
642 dim = isl_set_n_dim(set);
644 for (i = 0; i < hull->n_ineq; ++i) {
645 facet = compute_facet(set, hull->ineq[i]);
646 facet = isl_basic_set_add_equality(facet, hull->ineq[i]);
647 facet = isl_basic_set_gauss(facet, NULL);
648 facet = isl_basic_set_normalize_constraints(facet);
649 hull_facet = isl_basic_set_copy(hull);
650 hull_facet = isl_basic_set_add_equality(hull_facet, hull->ineq[i]);
651 hull_facet = isl_basic_set_gauss(hull_facet, NULL);
652 hull_facet = isl_basic_set_normalize_constraints(hull_facet);
653 if (!facet || !hull_facet)
654 goto error;
655 hull = isl_basic_set_cow(hull);
656 hull = isl_basic_set_extend_space(hull,
657 isl_space_copy(hull->dim), 0, 0, facet->n_ineq);
658 if (!hull)
659 goto error;
660 for (j = 0; j < facet->n_ineq; ++j) {
661 for (f = 0; f < hull_facet->n_ineq; ++f)
662 if (isl_seq_eq(facet->ineq[j],
663 hull_facet->ineq[f], 1 + dim))
664 break;
665 if (f < hull_facet->n_ineq)
666 continue;
667 k = isl_basic_set_alloc_inequality(hull);
668 if (k < 0)
669 goto error;
670 isl_seq_cpy(hull->ineq[k], hull->ineq[i], 1+dim);
671 if (!isl_set_wrap_facet(set, hull->ineq[k], facet->ineq[j]))
672 goto error;
674 isl_basic_set_free(hull_facet);
675 isl_basic_set_free(facet);
677 hull = isl_basic_set_simplify(hull);
678 hull = isl_basic_set_finalize(hull);
679 return hull;
680 error:
681 isl_basic_set_free(hull_facet);
682 isl_basic_set_free(facet);
683 isl_basic_set_free(hull);
684 return NULL;
687 /* Special case for computing the convex hull of a one dimensional set.
688 * We simply collect the lower and upper bounds of each basic set
689 * and the biggest of those.
691 static struct isl_basic_set *convex_hull_1d(struct isl_set *set)
693 struct isl_mat *c = NULL;
694 isl_int *lower = NULL;
695 isl_int *upper = NULL;
696 int i, j, k;
697 isl_int a, b;
698 struct isl_basic_set *hull;
700 for (i = 0; i < set->n; ++i) {
701 set->p[i] = isl_basic_set_simplify(set->p[i]);
702 if (!set->p[i])
703 goto error;
705 set = isl_set_remove_empty_parts(set);
706 if (!set)
707 goto error;
708 isl_assert(set->ctx, set->n > 0, goto error);
709 c = isl_mat_alloc(set->ctx, 2, 2);
710 if (!c)
711 goto error;
713 if (set->p[0]->n_eq > 0) {
714 isl_assert(set->ctx, set->p[0]->n_eq == 1, goto error);
715 lower = c->row[0];
716 upper = c->row[1];
717 if (isl_int_is_pos(set->p[0]->eq[0][1])) {
718 isl_seq_cpy(lower, set->p[0]->eq[0], 2);
719 isl_seq_neg(upper, set->p[0]->eq[0], 2);
720 } else {
721 isl_seq_neg(lower, set->p[0]->eq[0], 2);
722 isl_seq_cpy(upper, set->p[0]->eq[0], 2);
724 } else {
725 for (j = 0; j < set->p[0]->n_ineq; ++j) {
726 if (isl_int_is_pos(set->p[0]->ineq[j][1])) {
727 lower = c->row[0];
728 isl_seq_cpy(lower, set->p[0]->ineq[j], 2);
729 } else {
730 upper = c->row[1];
731 isl_seq_cpy(upper, set->p[0]->ineq[j], 2);
736 isl_int_init(a);
737 isl_int_init(b);
738 for (i = 0; i < set->n; ++i) {
739 struct isl_basic_set *bset = set->p[i];
740 int has_lower = 0;
741 int has_upper = 0;
743 for (j = 0; j < bset->n_eq; ++j) {
744 has_lower = 1;
745 has_upper = 1;
746 if (lower) {
747 isl_int_mul(a, lower[0], bset->eq[j][1]);
748 isl_int_mul(b, lower[1], bset->eq[j][0]);
749 if (isl_int_lt(a, b) && isl_int_is_pos(bset->eq[j][1]))
750 isl_seq_cpy(lower, bset->eq[j], 2);
751 if (isl_int_gt(a, b) && isl_int_is_neg(bset->eq[j][1]))
752 isl_seq_neg(lower, bset->eq[j], 2);
754 if (upper) {
755 isl_int_mul(a, upper[0], bset->eq[j][1]);
756 isl_int_mul(b, upper[1], bset->eq[j][0]);
757 if (isl_int_lt(a, b) && isl_int_is_pos(bset->eq[j][1]))
758 isl_seq_neg(upper, bset->eq[j], 2);
759 if (isl_int_gt(a, b) && isl_int_is_neg(bset->eq[j][1]))
760 isl_seq_cpy(upper, bset->eq[j], 2);
763 for (j = 0; j < bset->n_ineq; ++j) {
764 if (isl_int_is_pos(bset->ineq[j][1]))
765 has_lower = 1;
766 if (isl_int_is_neg(bset->ineq[j][1]))
767 has_upper = 1;
768 if (lower && isl_int_is_pos(bset->ineq[j][1])) {
769 isl_int_mul(a, lower[0], bset->ineq[j][1]);
770 isl_int_mul(b, lower[1], bset->ineq[j][0]);
771 if (isl_int_lt(a, b))
772 isl_seq_cpy(lower, bset->ineq[j], 2);
774 if (upper && isl_int_is_neg(bset->ineq[j][1])) {
775 isl_int_mul(a, upper[0], bset->ineq[j][1]);
776 isl_int_mul(b, upper[1], bset->ineq[j][0]);
777 if (isl_int_gt(a, b))
778 isl_seq_cpy(upper, bset->ineq[j], 2);
781 if (!has_lower)
782 lower = NULL;
783 if (!has_upper)
784 upper = NULL;
786 isl_int_clear(a);
787 isl_int_clear(b);
789 hull = isl_basic_set_alloc(set->ctx, 0, 1, 0, 0, 2);
790 hull = isl_basic_set_set_rational(hull);
791 if (!hull)
792 goto error;
793 if (lower) {
794 k = isl_basic_set_alloc_inequality(hull);
795 isl_seq_cpy(hull->ineq[k], lower, 2);
797 if (upper) {
798 k = isl_basic_set_alloc_inequality(hull);
799 isl_seq_cpy(hull->ineq[k], upper, 2);
801 hull = isl_basic_set_finalize(hull);
802 isl_set_free(set);
803 isl_mat_free(c);
804 return hull;
805 error:
806 isl_set_free(set);
807 isl_mat_free(c);
808 return NULL;
811 static struct isl_basic_set *convex_hull_0d(struct isl_set *set)
813 struct isl_basic_set *convex_hull;
815 if (!set)
816 return NULL;
818 if (isl_set_is_empty(set))
819 convex_hull = isl_basic_set_empty(isl_space_copy(set->dim));
820 else
821 convex_hull = isl_basic_set_universe(isl_space_copy(set->dim));
822 isl_set_free(set);
823 return convex_hull;
826 /* Compute the convex hull of a pair of basic sets without any parameters or
827 * integer divisions using Fourier-Motzkin elimination.
828 * The convex hull is the set of all points that can be written as
829 * the sum of points from both basic sets (in homogeneous coordinates).
830 * We set up the constraints in a space with dimensions for each of
831 * the three sets and then project out the dimensions corresponding
832 * to the two original basic sets, retaining only those corresponding
833 * to the convex hull.
835 static struct isl_basic_set *convex_hull_pair_elim(struct isl_basic_set *bset1,
836 struct isl_basic_set *bset2)
838 int i, j, k;
839 struct isl_basic_set *bset[2];
840 struct isl_basic_set *hull = NULL;
841 unsigned dim;
843 if (!bset1 || !bset2)
844 goto error;
846 dim = isl_basic_set_n_dim(bset1);
847 hull = isl_basic_set_alloc(bset1->ctx, 0, 2 + 3 * dim, 0,
848 1 + dim + bset1->n_eq + bset2->n_eq,
849 2 + bset1->n_ineq + bset2->n_ineq);
850 bset[0] = bset1;
851 bset[1] = bset2;
852 for (i = 0; i < 2; ++i) {
853 for (j = 0; j < bset[i]->n_eq; ++j) {
854 k = isl_basic_set_alloc_equality(hull);
855 if (k < 0)
856 goto error;
857 isl_seq_clr(hull->eq[k], (i+1) * (1+dim));
858 isl_seq_clr(hull->eq[k]+(i+2)*(1+dim), (1-i)*(1+dim));
859 isl_seq_cpy(hull->eq[k]+(i+1)*(1+dim), bset[i]->eq[j],
860 1+dim);
862 for (j = 0; j < bset[i]->n_ineq; ++j) {
863 k = isl_basic_set_alloc_inequality(hull);
864 if (k < 0)
865 goto error;
866 isl_seq_clr(hull->ineq[k], (i+1) * (1+dim));
867 isl_seq_clr(hull->ineq[k]+(i+2)*(1+dim), (1-i)*(1+dim));
868 isl_seq_cpy(hull->ineq[k]+(i+1)*(1+dim),
869 bset[i]->ineq[j], 1+dim);
871 k = isl_basic_set_alloc_inequality(hull);
872 if (k < 0)
873 goto error;
874 isl_seq_clr(hull->ineq[k], 1+2+3*dim);
875 isl_int_set_si(hull->ineq[k][(i+1)*(1+dim)], 1);
877 for (j = 0; j < 1+dim; ++j) {
878 k = isl_basic_set_alloc_equality(hull);
879 if (k < 0)
880 goto error;
881 isl_seq_clr(hull->eq[k], 1+2+3*dim);
882 isl_int_set_si(hull->eq[k][j], -1);
883 isl_int_set_si(hull->eq[k][1+dim+j], 1);
884 isl_int_set_si(hull->eq[k][2*(1+dim)+j], 1);
886 hull = isl_basic_set_set_rational(hull);
887 hull = isl_basic_set_remove_dims(hull, isl_dim_set, dim, 2*(1+dim));
888 hull = isl_basic_set_remove_redundancies(hull);
889 isl_basic_set_free(bset1);
890 isl_basic_set_free(bset2);
891 return hull;
892 error:
893 isl_basic_set_free(bset1);
894 isl_basic_set_free(bset2);
895 isl_basic_set_free(hull);
896 return NULL;
899 /* Is the set bounded for each value of the parameters?
901 int isl_basic_set_is_bounded(__isl_keep isl_basic_set *bset)
903 struct isl_tab *tab;
904 int bounded;
906 if (!bset)
907 return -1;
908 if (isl_basic_set_plain_is_empty(bset))
909 return 1;
911 tab = isl_tab_from_recession_cone(bset, 1);
912 bounded = isl_tab_cone_is_bounded(tab);
913 isl_tab_free(tab);
914 return bounded;
917 /* Is the image bounded for each value of the parameters and
918 * the domain variables?
920 int isl_basic_map_image_is_bounded(__isl_keep isl_basic_map *bmap)
922 unsigned nparam = isl_basic_map_dim(bmap, isl_dim_param);
923 unsigned n_in = isl_basic_map_dim(bmap, isl_dim_in);
924 int bounded;
926 bmap = isl_basic_map_copy(bmap);
927 bmap = isl_basic_map_cow(bmap);
928 bmap = isl_basic_map_move_dims(bmap, isl_dim_param, nparam,
929 isl_dim_in, 0, n_in);
930 bounded = isl_basic_set_is_bounded((isl_basic_set *)bmap);
931 isl_basic_map_free(bmap);
933 return bounded;
936 /* Is the set bounded for each value of the parameters?
938 int isl_set_is_bounded(__isl_keep isl_set *set)
940 int i;
942 if (!set)
943 return -1;
945 for (i = 0; i < set->n; ++i) {
946 int bounded = isl_basic_set_is_bounded(set->p[i]);
947 if (!bounded || bounded < 0)
948 return bounded;
950 return 1;
953 /* Compute the lineality space of the convex hull of bset1 and bset2.
955 * We first compute the intersection of the recession cone of bset1
956 * with the negative of the recession cone of bset2 and then compute
957 * the linear hull of the resulting cone.
959 static struct isl_basic_set *induced_lineality_space(
960 struct isl_basic_set *bset1, struct isl_basic_set *bset2)
962 int i, k;
963 struct isl_basic_set *lin = NULL;
964 unsigned dim;
966 if (!bset1 || !bset2)
967 goto error;
969 dim = isl_basic_set_total_dim(bset1);
970 lin = isl_basic_set_alloc_space(isl_basic_set_get_space(bset1), 0,
971 bset1->n_eq + bset2->n_eq,
972 bset1->n_ineq + bset2->n_ineq);
973 lin = isl_basic_set_set_rational(lin);
974 if (!lin)
975 goto error;
976 for (i = 0; i < bset1->n_eq; ++i) {
977 k = isl_basic_set_alloc_equality(lin);
978 if (k < 0)
979 goto error;
980 isl_int_set_si(lin->eq[k][0], 0);
981 isl_seq_cpy(lin->eq[k] + 1, bset1->eq[i] + 1, dim);
983 for (i = 0; i < bset1->n_ineq; ++i) {
984 k = isl_basic_set_alloc_inequality(lin);
985 if (k < 0)
986 goto error;
987 isl_int_set_si(lin->ineq[k][0], 0);
988 isl_seq_cpy(lin->ineq[k] + 1, bset1->ineq[i] + 1, dim);
990 for (i = 0; i < bset2->n_eq; ++i) {
991 k = isl_basic_set_alloc_equality(lin);
992 if (k < 0)
993 goto error;
994 isl_int_set_si(lin->eq[k][0], 0);
995 isl_seq_neg(lin->eq[k] + 1, bset2->eq[i] + 1, dim);
997 for (i = 0; i < bset2->n_ineq; ++i) {
998 k = isl_basic_set_alloc_inequality(lin);
999 if (k < 0)
1000 goto error;
1001 isl_int_set_si(lin->ineq[k][0], 0);
1002 isl_seq_neg(lin->ineq[k] + 1, bset2->ineq[i] + 1, dim);
1005 isl_basic_set_free(bset1);
1006 isl_basic_set_free(bset2);
1007 return isl_basic_set_affine_hull(lin);
1008 error:
1009 isl_basic_set_free(lin);
1010 isl_basic_set_free(bset1);
1011 isl_basic_set_free(bset2);
1012 return NULL;
1015 static struct isl_basic_set *uset_convex_hull(struct isl_set *set);
1017 /* Given a set and a linear space "lin" of dimension n > 0,
1018 * project the linear space from the set, compute the convex hull
1019 * and then map the set back to the original space.
1021 * Let
1023 * M x = 0
1025 * describe the linear space. We first compute the Hermite normal
1026 * form H = M U of M = H Q, to obtain
1028 * H Q x = 0
1030 * The last n rows of H will be zero, so the last n variables of x' = Q x
1031 * are the one we want to project out. We do this by transforming each
1032 * basic set A x >= b to A U x' >= b and then removing the last n dimensions.
1033 * After computing the convex hull in x'_1, i.e., A' x'_1 >= b',
1034 * we transform the hull back to the original space as A' Q_1 x >= b',
1035 * with Q_1 all but the last n rows of Q.
1037 static struct isl_basic_set *modulo_lineality(struct isl_set *set,
1038 struct isl_basic_set *lin)
1040 unsigned total = isl_basic_set_total_dim(lin);
1041 unsigned lin_dim;
1042 struct isl_basic_set *hull;
1043 struct isl_mat *M, *U, *Q;
1045 if (!set || !lin)
1046 goto error;
1047 lin_dim = total - lin->n_eq;
1048 M = isl_mat_sub_alloc6(set->ctx, lin->eq, 0, lin->n_eq, 1, total);
1049 M = isl_mat_left_hermite(M, 0, &U, &Q);
1050 if (!M)
1051 goto error;
1052 isl_mat_free(M);
1053 isl_basic_set_free(lin);
1055 Q = isl_mat_drop_rows(Q, Q->n_row - lin_dim, lin_dim);
1057 U = isl_mat_lin_to_aff(U);
1058 Q = isl_mat_lin_to_aff(Q);
1060 set = isl_set_preimage(set, U);
1061 set = isl_set_remove_dims(set, isl_dim_set, total - lin_dim, lin_dim);
1062 hull = uset_convex_hull(set);
1063 hull = isl_basic_set_preimage(hull, Q);
1065 return hull;
1066 error:
1067 isl_basic_set_free(lin);
1068 isl_set_free(set);
1069 return NULL;
1072 /* Given two polyhedra with as constraints h_{ij} x >= 0 in homegeneous space,
1073 * set up an LP for solving
1075 * \sum_j \alpha_{1j} h_{1j} = \sum_j \alpha_{2j} h_{2j}
1077 * \alpha{i0} corresponds to the (implicit) positivity constraint 1 >= 0
1078 * The next \alpha{ij} correspond to the equalities and come in pairs.
1079 * The final \alpha{ij} correspond to the inequalities.
1081 static struct isl_basic_set *valid_direction_lp(
1082 struct isl_basic_set *bset1, struct isl_basic_set *bset2)
1084 isl_space *dim;
1085 struct isl_basic_set *lp;
1086 unsigned d;
1087 int n;
1088 int i, j, k;
1090 if (!bset1 || !bset2)
1091 goto error;
1092 d = 1 + isl_basic_set_total_dim(bset1);
1093 n = 2 +
1094 2 * bset1->n_eq + bset1->n_ineq + 2 * bset2->n_eq + bset2->n_ineq;
1095 dim = isl_space_set_alloc(bset1->ctx, 0, n);
1096 lp = isl_basic_set_alloc_space(dim, 0, d, n);
1097 if (!lp)
1098 goto error;
1099 for (i = 0; i < n; ++i) {
1100 k = isl_basic_set_alloc_inequality(lp);
1101 if (k < 0)
1102 goto error;
1103 isl_seq_clr(lp->ineq[k] + 1, n);
1104 isl_int_set_si(lp->ineq[k][0], -1);
1105 isl_int_set_si(lp->ineq[k][1 + i], 1);
1107 for (i = 0; i < d; ++i) {
1108 k = isl_basic_set_alloc_equality(lp);
1109 if (k < 0)
1110 goto error;
1111 n = 0;
1112 isl_int_set_si(lp->eq[k][n], 0); n++;
1113 /* positivity constraint 1 >= 0 */
1114 isl_int_set_si(lp->eq[k][n], i == 0); n++;
1115 for (j = 0; j < bset1->n_eq; ++j) {
1116 isl_int_set(lp->eq[k][n], bset1->eq[j][i]); n++;
1117 isl_int_neg(lp->eq[k][n], bset1->eq[j][i]); n++;
1119 for (j = 0; j < bset1->n_ineq; ++j) {
1120 isl_int_set(lp->eq[k][n], bset1->ineq[j][i]); n++;
1122 /* positivity constraint 1 >= 0 */
1123 isl_int_set_si(lp->eq[k][n], -(i == 0)); n++;
1124 for (j = 0; j < bset2->n_eq; ++j) {
1125 isl_int_neg(lp->eq[k][n], bset2->eq[j][i]); n++;
1126 isl_int_set(lp->eq[k][n], bset2->eq[j][i]); n++;
1128 for (j = 0; j < bset2->n_ineq; ++j) {
1129 isl_int_neg(lp->eq[k][n], bset2->ineq[j][i]); n++;
1132 lp = isl_basic_set_gauss(lp, NULL);
1133 isl_basic_set_free(bset1);
1134 isl_basic_set_free(bset2);
1135 return lp;
1136 error:
1137 isl_basic_set_free(bset1);
1138 isl_basic_set_free(bset2);
1139 return NULL;
1142 /* Compute a vector s in the homogeneous space such that <s, r> > 0
1143 * for all rays in the homogeneous space of the two cones that correspond
1144 * to the input polyhedra bset1 and bset2.
1146 * We compute s as a vector that satisfies
1148 * s = \sum_j \alpha_{ij} h_{ij} for i = 1,2 (*)
1150 * with h_{ij} the normals of the facets of polyhedron i
1151 * (including the "positivity constraint" 1 >= 0) and \alpha_{ij}
1152 * strictly positive numbers. For simplicity we impose \alpha_{ij} >= 1.
1153 * We first set up an LP with as variables the \alpha{ij}.
1154 * In this formulation, for each polyhedron i,
1155 * the first constraint is the positivity constraint, followed by pairs
1156 * of variables for the equalities, followed by variables for the inequalities.
1157 * We then simply pick a feasible solution and compute s using (*).
1159 * Note that we simply pick any valid direction and make no attempt
1160 * to pick a "good" or even the "best" valid direction.
1162 static struct isl_vec *valid_direction(
1163 struct isl_basic_set *bset1, struct isl_basic_set *bset2)
1165 struct isl_basic_set *lp;
1166 struct isl_tab *tab;
1167 struct isl_vec *sample = NULL;
1168 struct isl_vec *dir;
1169 unsigned d;
1170 int i;
1171 int n;
1173 if (!bset1 || !bset2)
1174 goto error;
1175 lp = valid_direction_lp(isl_basic_set_copy(bset1),
1176 isl_basic_set_copy(bset2));
1177 tab = isl_tab_from_basic_set(lp, 0);
1178 sample = isl_tab_get_sample_value(tab);
1179 isl_tab_free(tab);
1180 isl_basic_set_free(lp);
1181 if (!sample)
1182 goto error;
1183 d = isl_basic_set_total_dim(bset1);
1184 dir = isl_vec_alloc(bset1->ctx, 1 + d);
1185 if (!dir)
1186 goto error;
1187 isl_seq_clr(dir->block.data + 1, dir->size - 1);
1188 n = 1;
1189 /* positivity constraint 1 >= 0 */
1190 isl_int_set(dir->block.data[0], sample->block.data[n]); n++;
1191 for (i = 0; i < bset1->n_eq; ++i) {
1192 isl_int_sub(sample->block.data[n],
1193 sample->block.data[n], sample->block.data[n+1]);
1194 isl_seq_combine(dir->block.data,
1195 bset1->ctx->one, dir->block.data,
1196 sample->block.data[n], bset1->eq[i], 1 + d);
1198 n += 2;
1200 for (i = 0; i < bset1->n_ineq; ++i)
1201 isl_seq_combine(dir->block.data,
1202 bset1->ctx->one, dir->block.data,
1203 sample->block.data[n++], bset1->ineq[i], 1 + d);
1204 isl_vec_free(sample);
1205 isl_seq_normalize(bset1->ctx, dir->el, dir->size);
1206 isl_basic_set_free(bset1);
1207 isl_basic_set_free(bset2);
1208 return dir;
1209 error:
1210 isl_vec_free(sample);
1211 isl_basic_set_free(bset1);
1212 isl_basic_set_free(bset2);
1213 return NULL;
1216 /* Given a polyhedron b_i + A_i x >= 0 and a map T = S^{-1},
1217 * compute b_i' + A_i' x' >= 0, with
1219 * [ b_i A_i ] [ y' ] [ y' ]
1220 * [ 1 0 ] S^{-1} [ x' ] >= 0 or [ b_i' A_i' ] [ x' ] >= 0
1222 * In particular, add the "positivity constraint" and then perform
1223 * the mapping.
1225 static struct isl_basic_set *homogeneous_map(struct isl_basic_set *bset,
1226 struct isl_mat *T)
1228 int k;
1230 if (!bset)
1231 goto error;
1232 bset = isl_basic_set_extend_constraints(bset, 0, 1);
1233 k = isl_basic_set_alloc_inequality(bset);
1234 if (k < 0)
1235 goto error;
1236 isl_seq_clr(bset->ineq[k] + 1, isl_basic_set_total_dim(bset));
1237 isl_int_set_si(bset->ineq[k][0], 1);
1238 bset = isl_basic_set_preimage(bset, T);
1239 return bset;
1240 error:
1241 isl_mat_free(T);
1242 isl_basic_set_free(bset);
1243 return NULL;
1246 /* Compute the convex hull of a pair of basic sets without any parameters or
1247 * integer divisions, where the convex hull is known to be pointed,
1248 * but the basic sets may be unbounded.
1250 * We turn this problem into the computation of a convex hull of a pair
1251 * _bounded_ polyhedra by "changing the direction of the homogeneous
1252 * dimension". This idea is due to Matthias Koeppe.
1254 * Consider the cones in homogeneous space that correspond to the
1255 * input polyhedra. The rays of these cones are also rays of the
1256 * polyhedra if the coordinate that corresponds to the homogeneous
1257 * dimension is zero. That is, if the inner product of the rays
1258 * with the homogeneous direction is zero.
1259 * The cones in the homogeneous space can also be considered to
1260 * correspond to other pairs of polyhedra by chosing a different
1261 * homogeneous direction. To ensure that both of these polyhedra
1262 * are bounded, we need to make sure that all rays of the cones
1263 * correspond to vertices and not to rays.
1264 * Let s be a direction such that <s, r> > 0 for all rays r of both cones.
1265 * Then using s as a homogeneous direction, we obtain a pair of polytopes.
1266 * The vector s is computed in valid_direction.
1268 * Note that we need to consider _all_ rays of the cones and not just
1269 * the rays that correspond to rays in the polyhedra. If we were to
1270 * only consider those rays and turn them into vertices, then we
1271 * may inadvertently turn some vertices into rays.
1273 * The standard homogeneous direction is the unit vector in the 0th coordinate.
1274 * We therefore transform the two polyhedra such that the selected
1275 * direction is mapped onto this standard direction and then proceed
1276 * with the normal computation.
1277 * Let S be a non-singular square matrix with s as its first row,
1278 * then we want to map the polyhedra to the space
1280 * [ y' ] [ y ] [ y ] [ y' ]
1281 * [ x' ] = S [ x ] i.e., [ x ] = S^{-1} [ x' ]
1283 * We take S to be the unimodular completion of s to limit the growth
1284 * of the coefficients in the following computations.
1286 * Let b_i + A_i x >= 0 be the constraints of polyhedron i.
1287 * We first move to the homogeneous dimension
1289 * b_i y + A_i x >= 0 [ b_i A_i ] [ y ] [ 0 ]
1290 * y >= 0 or [ 1 0 ] [ x ] >= [ 0 ]
1292 * Then we change directoin
1294 * [ b_i A_i ] [ y' ] [ y' ]
1295 * [ 1 0 ] S^{-1} [ x' ] >= 0 or [ b_i' A_i' ] [ x' ] >= 0
1297 * Then we compute the convex hull of the polytopes b_i' + A_i' x' >= 0
1298 * resulting in b' + A' x' >= 0, which we then convert back
1300 * [ y ] [ y ]
1301 * [ b' A' ] S [ x ] >= 0 or [ b A ] [ x ] >= 0
1303 * The polyhedron b + A x >= 0 is then the convex hull of the input polyhedra.
1305 static struct isl_basic_set *convex_hull_pair_pointed(
1306 struct isl_basic_set *bset1, struct isl_basic_set *bset2)
1308 struct isl_ctx *ctx = NULL;
1309 struct isl_vec *dir = NULL;
1310 struct isl_mat *T = NULL;
1311 struct isl_mat *T2 = NULL;
1312 struct isl_basic_set *hull;
1313 struct isl_set *set;
1315 if (!bset1 || !bset2)
1316 goto error;
1317 ctx = bset1->ctx;
1318 dir = valid_direction(isl_basic_set_copy(bset1),
1319 isl_basic_set_copy(bset2));
1320 if (!dir)
1321 goto error;
1322 T = isl_mat_alloc(bset1->ctx, dir->size, dir->size);
1323 if (!T)
1324 goto error;
1325 isl_seq_cpy(T->row[0], dir->block.data, dir->size);
1326 T = isl_mat_unimodular_complete(T, 1);
1327 T2 = isl_mat_right_inverse(isl_mat_copy(T));
1329 bset1 = homogeneous_map(bset1, isl_mat_copy(T2));
1330 bset2 = homogeneous_map(bset2, T2);
1331 set = isl_set_alloc_space(isl_basic_set_get_space(bset1), 2, 0);
1332 set = isl_set_add_basic_set(set, bset1);
1333 set = isl_set_add_basic_set(set, bset2);
1334 hull = uset_convex_hull(set);
1335 hull = isl_basic_set_preimage(hull, T);
1337 isl_vec_free(dir);
1339 return hull;
1340 error:
1341 isl_vec_free(dir);
1342 isl_basic_set_free(bset1);
1343 isl_basic_set_free(bset2);
1344 return NULL;
1347 static struct isl_basic_set *uset_convex_hull_wrap(struct isl_set *set);
1348 static struct isl_basic_set *modulo_affine_hull(
1349 struct isl_set *set, struct isl_basic_set *affine_hull);
1351 /* Compute the convex hull of a pair of basic sets without any parameters or
1352 * integer divisions.
1354 * This function is called from uset_convex_hull_unbounded, which
1355 * means that the complete convex hull is unbounded. Some pairs
1356 * of basic sets may still be bounded, though.
1357 * They may even lie inside a lower dimensional space, in which
1358 * case they need to be handled inside their affine hull since
1359 * the main algorithm assumes that the result is full-dimensional.
1361 * If the convex hull of the two basic sets would have a non-trivial
1362 * lineality space, we first project out this lineality space.
1364 static struct isl_basic_set *convex_hull_pair(struct isl_basic_set *bset1,
1365 struct isl_basic_set *bset2)
1367 isl_basic_set *lin, *aff;
1368 int bounded1, bounded2;
1370 if (bset1->ctx->opt->convex == ISL_CONVEX_HULL_FM)
1371 return convex_hull_pair_elim(bset1, bset2);
1373 aff = isl_set_affine_hull(isl_basic_set_union(isl_basic_set_copy(bset1),
1374 isl_basic_set_copy(bset2)));
1375 if (!aff)
1376 goto error;
1377 if (aff->n_eq != 0)
1378 return modulo_affine_hull(isl_basic_set_union(bset1, bset2), aff);
1379 isl_basic_set_free(aff);
1381 bounded1 = isl_basic_set_is_bounded(bset1);
1382 bounded2 = isl_basic_set_is_bounded(bset2);
1384 if (bounded1 < 0 || bounded2 < 0)
1385 goto error;
1387 if (bounded1 && bounded2)
1388 uset_convex_hull_wrap(isl_basic_set_union(bset1, bset2));
1390 if (bounded1 || bounded2)
1391 return convex_hull_pair_pointed(bset1, bset2);
1393 lin = induced_lineality_space(isl_basic_set_copy(bset1),
1394 isl_basic_set_copy(bset2));
1395 if (!lin)
1396 goto error;
1397 if (isl_basic_set_is_universe(lin)) {
1398 isl_basic_set_free(bset1);
1399 isl_basic_set_free(bset2);
1400 return lin;
1402 if (lin->n_eq < isl_basic_set_total_dim(lin)) {
1403 struct isl_set *set;
1404 set = isl_set_alloc_space(isl_basic_set_get_space(bset1), 2, 0);
1405 set = isl_set_add_basic_set(set, bset1);
1406 set = isl_set_add_basic_set(set, bset2);
1407 return modulo_lineality(set, lin);
1409 isl_basic_set_free(lin);
1411 return convex_hull_pair_pointed(bset1, bset2);
1412 error:
1413 isl_basic_set_free(bset1);
1414 isl_basic_set_free(bset2);
1415 return NULL;
1418 /* Compute the lineality space of a basic set.
1419 * We currently do not allow the basic set to have any divs.
1420 * We basically just drop the constants and turn every inequality
1421 * into an equality.
1423 struct isl_basic_set *isl_basic_set_lineality_space(struct isl_basic_set *bset)
1425 int i, k;
1426 struct isl_basic_set *lin = NULL;
1427 unsigned dim;
1429 if (!bset)
1430 goto error;
1431 isl_assert(bset->ctx, bset->n_div == 0, goto error);
1432 dim = isl_basic_set_total_dim(bset);
1434 lin = isl_basic_set_alloc_space(isl_basic_set_get_space(bset), 0, dim, 0);
1435 if (!lin)
1436 goto error;
1437 for (i = 0; i < bset->n_eq; ++i) {
1438 k = isl_basic_set_alloc_equality(lin);
1439 if (k < 0)
1440 goto error;
1441 isl_int_set_si(lin->eq[k][0], 0);
1442 isl_seq_cpy(lin->eq[k] + 1, bset->eq[i] + 1, dim);
1444 lin = isl_basic_set_gauss(lin, NULL);
1445 if (!lin)
1446 goto error;
1447 for (i = 0; i < bset->n_ineq && lin->n_eq < dim; ++i) {
1448 k = isl_basic_set_alloc_equality(lin);
1449 if (k < 0)
1450 goto error;
1451 isl_int_set_si(lin->eq[k][0], 0);
1452 isl_seq_cpy(lin->eq[k] + 1, bset->ineq[i] + 1, dim);
1453 lin = isl_basic_set_gauss(lin, NULL);
1454 if (!lin)
1455 goto error;
1457 isl_basic_set_free(bset);
1458 return lin;
1459 error:
1460 isl_basic_set_free(lin);
1461 isl_basic_set_free(bset);
1462 return NULL;
1465 /* Compute the (linear) hull of the lineality spaces of the basic sets in the
1466 * "underlying" set "set".
1468 static struct isl_basic_set *uset_combined_lineality_space(struct isl_set *set)
1470 int i;
1471 struct isl_set *lin = NULL;
1473 if (!set)
1474 return NULL;
1475 if (set->n == 0) {
1476 isl_space *dim = isl_set_get_space(set);
1477 isl_set_free(set);
1478 return isl_basic_set_empty(dim);
1481 lin = isl_set_alloc_space(isl_set_get_space(set), set->n, 0);
1482 for (i = 0; i < set->n; ++i)
1483 lin = isl_set_add_basic_set(lin,
1484 isl_basic_set_lineality_space(isl_basic_set_copy(set->p[i])));
1485 isl_set_free(set);
1486 return isl_set_affine_hull(lin);
1489 /* Compute the convex hull of a set without any parameters or
1490 * integer divisions.
1491 * In each step, we combined two basic sets until only one
1492 * basic set is left.
1493 * The input basic sets are assumed not to have a non-trivial
1494 * lineality space. If any of the intermediate results has
1495 * a non-trivial lineality space, it is projected out.
1497 static struct isl_basic_set *uset_convex_hull_unbounded(struct isl_set *set)
1499 struct isl_basic_set *convex_hull = NULL;
1501 convex_hull = isl_set_copy_basic_set(set);
1502 set = isl_set_drop_basic_set(set, convex_hull);
1503 if (!set)
1504 goto error;
1505 while (set->n > 0) {
1506 struct isl_basic_set *t;
1507 t = isl_set_copy_basic_set(set);
1508 if (!t)
1509 goto error;
1510 set = isl_set_drop_basic_set(set, t);
1511 if (!set)
1512 goto error;
1513 convex_hull = convex_hull_pair(convex_hull, t);
1514 if (set->n == 0)
1515 break;
1516 t = isl_basic_set_lineality_space(isl_basic_set_copy(convex_hull));
1517 if (!t)
1518 goto error;
1519 if (isl_basic_set_is_universe(t)) {
1520 isl_basic_set_free(convex_hull);
1521 convex_hull = t;
1522 break;
1524 if (t->n_eq < isl_basic_set_total_dim(t)) {
1525 set = isl_set_add_basic_set(set, convex_hull);
1526 return modulo_lineality(set, t);
1528 isl_basic_set_free(t);
1530 isl_set_free(set);
1531 return convex_hull;
1532 error:
1533 isl_set_free(set);
1534 isl_basic_set_free(convex_hull);
1535 return NULL;
1538 /* Compute an initial hull for wrapping containing a single initial
1539 * facet.
1540 * This function assumes that the given set is bounded.
1542 static struct isl_basic_set *initial_hull(struct isl_basic_set *hull,
1543 struct isl_set *set)
1545 struct isl_mat *bounds = NULL;
1546 unsigned dim;
1547 int k;
1549 if (!hull)
1550 goto error;
1551 bounds = initial_facet_constraint(set);
1552 if (!bounds)
1553 goto error;
1554 k = isl_basic_set_alloc_inequality(hull);
1555 if (k < 0)
1556 goto error;
1557 dim = isl_set_n_dim(set);
1558 isl_assert(set->ctx, 1 + dim == bounds->n_col, goto error);
1559 isl_seq_cpy(hull->ineq[k], bounds->row[0], bounds->n_col);
1560 isl_mat_free(bounds);
1562 return hull;
1563 error:
1564 isl_basic_set_free(hull);
1565 isl_mat_free(bounds);
1566 return NULL;
1569 struct max_constraint {
1570 struct isl_mat *c;
1571 int count;
1572 int ineq;
1575 static int max_constraint_equal(const void *entry, const void *val)
1577 struct max_constraint *a = (struct max_constraint *)entry;
1578 isl_int *b = (isl_int *)val;
1580 return isl_seq_eq(a->c->row[0] + 1, b, a->c->n_col - 1);
1583 static void update_constraint(struct isl_ctx *ctx, struct isl_hash_table *table,
1584 isl_int *con, unsigned len, int n, int ineq)
1586 struct isl_hash_table_entry *entry;
1587 struct max_constraint *c;
1588 uint32_t c_hash;
1590 c_hash = isl_seq_get_hash(con + 1, len);
1591 entry = isl_hash_table_find(ctx, table, c_hash, max_constraint_equal,
1592 con + 1, 0);
1593 if (!entry)
1594 return;
1595 c = entry->data;
1596 if (c->count < n) {
1597 isl_hash_table_remove(ctx, table, entry);
1598 return;
1600 c->count++;
1601 if (isl_int_gt(c->c->row[0][0], con[0]))
1602 return;
1603 if (isl_int_eq(c->c->row[0][0], con[0])) {
1604 if (ineq)
1605 c->ineq = ineq;
1606 return;
1608 c->c = isl_mat_cow(c->c);
1609 isl_int_set(c->c->row[0][0], con[0]);
1610 c->ineq = ineq;
1613 /* Check whether the constraint hash table "table" constains the constraint
1614 * "con".
1616 static int has_constraint(struct isl_ctx *ctx, struct isl_hash_table *table,
1617 isl_int *con, unsigned len, int n)
1619 struct isl_hash_table_entry *entry;
1620 struct max_constraint *c;
1621 uint32_t c_hash;
1623 c_hash = isl_seq_get_hash(con + 1, len);
1624 entry = isl_hash_table_find(ctx, table, c_hash, max_constraint_equal,
1625 con + 1, 0);
1626 if (!entry)
1627 return 0;
1628 c = entry->data;
1629 if (c->count < n)
1630 return 0;
1631 return isl_int_eq(c->c->row[0][0], con[0]);
1634 /* Check for inequality constraints of a basic set without equalities
1635 * such that the same or more stringent copies of the constraint appear
1636 * in all of the basic sets. Such constraints are necessarily facet
1637 * constraints of the convex hull.
1639 * If the resulting basic set is by chance identical to one of
1640 * the basic sets in "set", then we know that this basic set contains
1641 * all other basic sets and is therefore the convex hull of set.
1642 * In this case we set *is_hull to 1.
1644 static struct isl_basic_set *common_constraints(struct isl_basic_set *hull,
1645 struct isl_set *set, int *is_hull)
1647 int i, j, s, n;
1648 int min_constraints;
1649 int best;
1650 struct max_constraint *constraints = NULL;
1651 struct isl_hash_table *table = NULL;
1652 unsigned total;
1654 *is_hull = 0;
1656 for (i = 0; i < set->n; ++i)
1657 if (set->p[i]->n_eq == 0)
1658 break;
1659 if (i >= set->n)
1660 return hull;
1661 min_constraints = set->p[i]->n_ineq;
1662 best = i;
1663 for (i = best + 1; i < set->n; ++i) {
1664 if (set->p[i]->n_eq != 0)
1665 continue;
1666 if (set->p[i]->n_ineq >= min_constraints)
1667 continue;
1668 min_constraints = set->p[i]->n_ineq;
1669 best = i;
1671 constraints = isl_calloc_array(hull->ctx, struct max_constraint,
1672 min_constraints);
1673 if (!constraints)
1674 return hull;
1675 table = isl_alloc_type(hull->ctx, struct isl_hash_table);
1676 if (isl_hash_table_init(hull->ctx, table, min_constraints))
1677 goto error;
1679 total = isl_space_dim(set->dim, isl_dim_all);
1680 for (i = 0; i < set->p[best]->n_ineq; ++i) {
1681 constraints[i].c = isl_mat_sub_alloc6(hull->ctx,
1682 set->p[best]->ineq + i, 0, 1, 0, 1 + total);
1683 if (!constraints[i].c)
1684 goto error;
1685 constraints[i].ineq = 1;
1687 for (i = 0; i < min_constraints; ++i) {
1688 struct isl_hash_table_entry *entry;
1689 uint32_t c_hash;
1690 c_hash = isl_seq_get_hash(constraints[i].c->row[0] + 1, total);
1691 entry = isl_hash_table_find(hull->ctx, table, c_hash,
1692 max_constraint_equal, constraints[i].c->row[0] + 1, 1);
1693 if (!entry)
1694 goto error;
1695 isl_assert(hull->ctx, !entry->data, goto error);
1696 entry->data = &constraints[i];
1699 n = 0;
1700 for (s = 0; s < set->n; ++s) {
1701 if (s == best)
1702 continue;
1704 for (i = 0; i < set->p[s]->n_eq; ++i) {
1705 isl_int *eq = set->p[s]->eq[i];
1706 for (j = 0; j < 2; ++j) {
1707 isl_seq_neg(eq, eq, 1 + total);
1708 update_constraint(hull->ctx, table,
1709 eq, total, n, 0);
1712 for (i = 0; i < set->p[s]->n_ineq; ++i) {
1713 isl_int *ineq = set->p[s]->ineq[i];
1714 update_constraint(hull->ctx, table, ineq, total, n,
1715 set->p[s]->n_eq == 0);
1717 ++n;
1720 for (i = 0; i < min_constraints; ++i) {
1721 if (constraints[i].count < n)
1722 continue;
1723 if (!constraints[i].ineq)
1724 continue;
1725 j = isl_basic_set_alloc_inequality(hull);
1726 if (j < 0)
1727 goto error;
1728 isl_seq_cpy(hull->ineq[j], constraints[i].c->row[0], 1 + total);
1731 for (s = 0; s < set->n; ++s) {
1732 if (set->p[s]->n_eq)
1733 continue;
1734 if (set->p[s]->n_ineq != hull->n_ineq)
1735 continue;
1736 for (i = 0; i < set->p[s]->n_ineq; ++i) {
1737 isl_int *ineq = set->p[s]->ineq[i];
1738 if (!has_constraint(hull->ctx, table, ineq, total, n))
1739 break;
1741 if (i == set->p[s]->n_ineq)
1742 *is_hull = 1;
1745 isl_hash_table_clear(table);
1746 for (i = 0; i < min_constraints; ++i)
1747 isl_mat_free(constraints[i].c);
1748 free(constraints);
1749 free(table);
1750 return hull;
1751 error:
1752 isl_hash_table_clear(table);
1753 free(table);
1754 if (constraints)
1755 for (i = 0; i < min_constraints; ++i)
1756 isl_mat_free(constraints[i].c);
1757 free(constraints);
1758 return hull;
1761 /* Create a template for the convex hull of "set" and fill it up
1762 * obvious facet constraints, if any. If the result happens to
1763 * be the convex hull of "set" then *is_hull is set to 1.
1765 static struct isl_basic_set *proto_hull(struct isl_set *set, int *is_hull)
1767 struct isl_basic_set *hull;
1768 unsigned n_ineq;
1769 int i;
1771 n_ineq = 1;
1772 for (i = 0; i < set->n; ++i) {
1773 n_ineq += set->p[i]->n_eq;
1774 n_ineq += set->p[i]->n_ineq;
1776 hull = isl_basic_set_alloc_space(isl_space_copy(set->dim), 0, 0, n_ineq);
1777 hull = isl_basic_set_set_rational(hull);
1778 if (!hull)
1779 return NULL;
1780 return common_constraints(hull, set, is_hull);
1783 static struct isl_basic_set *uset_convex_hull_wrap(struct isl_set *set)
1785 struct isl_basic_set *hull;
1786 int is_hull;
1788 hull = proto_hull(set, &is_hull);
1789 if (hull && !is_hull) {
1790 if (hull->n_ineq == 0)
1791 hull = initial_hull(hull, set);
1792 hull = extend(hull, set);
1794 isl_set_free(set);
1796 return hull;
1799 /* Compute the convex hull of a set without any parameters or
1800 * integer divisions. Depending on whether the set is bounded,
1801 * we pass control to the wrapping based convex hull or
1802 * the Fourier-Motzkin elimination based convex hull.
1803 * We also handle a few special cases before checking the boundedness.
1805 static struct isl_basic_set *uset_convex_hull(struct isl_set *set)
1807 struct isl_basic_set *convex_hull = NULL;
1808 struct isl_basic_set *lin;
1810 if (isl_set_n_dim(set) == 0)
1811 return convex_hull_0d(set);
1813 set = isl_set_coalesce(set);
1814 set = isl_set_set_rational(set);
1816 if (!set)
1817 goto error;
1818 if (!set)
1819 return NULL;
1820 if (set->n == 1) {
1821 convex_hull = isl_basic_set_copy(set->p[0]);
1822 isl_set_free(set);
1823 return convex_hull;
1825 if (isl_set_n_dim(set) == 1)
1826 return convex_hull_1d(set);
1828 if (isl_set_is_bounded(set) &&
1829 set->ctx->opt->convex == ISL_CONVEX_HULL_WRAP)
1830 return uset_convex_hull_wrap(set);
1832 lin = uset_combined_lineality_space(isl_set_copy(set));
1833 if (!lin)
1834 goto error;
1835 if (isl_basic_set_is_universe(lin)) {
1836 isl_set_free(set);
1837 return lin;
1839 if (lin->n_eq < isl_basic_set_total_dim(lin))
1840 return modulo_lineality(set, lin);
1841 isl_basic_set_free(lin);
1843 return uset_convex_hull_unbounded(set);
1844 error:
1845 isl_set_free(set);
1846 isl_basic_set_free(convex_hull);
1847 return NULL;
1850 /* This is the core procedure, where "set" is a "pure" set, i.e.,
1851 * without parameters or divs and where the convex hull of set is
1852 * known to be full-dimensional.
1854 static struct isl_basic_set *uset_convex_hull_wrap_bounded(struct isl_set *set)
1856 struct isl_basic_set *convex_hull = NULL;
1858 if (!set)
1859 goto error;
1861 if (isl_set_n_dim(set) == 0) {
1862 convex_hull = isl_basic_set_universe(isl_space_copy(set->dim));
1863 isl_set_free(set);
1864 convex_hull = isl_basic_set_set_rational(convex_hull);
1865 return convex_hull;
1868 set = isl_set_set_rational(set);
1869 set = isl_set_coalesce(set);
1870 if (!set)
1871 goto error;
1872 if (set->n == 1) {
1873 convex_hull = isl_basic_set_copy(set->p[0]);
1874 isl_set_free(set);
1875 convex_hull = isl_basic_map_remove_redundancies(convex_hull);
1876 return convex_hull;
1878 if (isl_set_n_dim(set) == 1)
1879 return convex_hull_1d(set);
1881 return uset_convex_hull_wrap(set);
1882 error:
1883 isl_set_free(set);
1884 return NULL;
1887 /* Compute the convex hull of set "set" with affine hull "affine_hull",
1888 * We first remove the equalities (transforming the set), compute the
1889 * convex hull of the transformed set and then add the equalities back
1890 * (after performing the inverse transformation.
1892 static struct isl_basic_set *modulo_affine_hull(
1893 struct isl_set *set, struct isl_basic_set *affine_hull)
1895 struct isl_mat *T;
1896 struct isl_mat *T2;
1897 struct isl_basic_set *dummy;
1898 struct isl_basic_set *convex_hull;
1900 dummy = isl_basic_set_remove_equalities(
1901 isl_basic_set_copy(affine_hull), &T, &T2);
1902 if (!dummy)
1903 goto error;
1904 isl_basic_set_free(dummy);
1905 set = isl_set_preimage(set, T);
1906 convex_hull = uset_convex_hull(set);
1907 convex_hull = isl_basic_set_preimage(convex_hull, T2);
1908 convex_hull = isl_basic_set_intersect(convex_hull, affine_hull);
1909 return convex_hull;
1910 error:
1911 isl_basic_set_free(affine_hull);
1912 isl_set_free(set);
1913 return NULL;
1916 /* Compute the convex hull of a map.
1918 * The implementation was inspired by "Extended Convex Hull" by Fukuda et al.,
1919 * specifically, the wrapping of facets to obtain new facets.
1921 struct isl_basic_map *isl_map_convex_hull(struct isl_map *map)
1923 struct isl_basic_set *bset;
1924 struct isl_basic_map *model = NULL;
1925 struct isl_basic_set *affine_hull = NULL;
1926 struct isl_basic_map *convex_hull = NULL;
1927 struct isl_set *set = NULL;
1928 struct isl_ctx *ctx;
1930 if (!map)
1931 goto error;
1933 ctx = map->ctx;
1934 if (map->n == 0) {
1935 convex_hull = isl_basic_map_empty_like_map(map);
1936 isl_map_free(map);
1937 return convex_hull;
1940 map = isl_map_detect_equalities(map);
1941 map = isl_map_align_divs(map);
1942 if (!map)
1943 goto error;
1944 model = isl_basic_map_copy(map->p[0]);
1945 set = isl_map_underlying_set(map);
1946 if (!set)
1947 goto error;
1949 affine_hull = isl_set_affine_hull(isl_set_copy(set));
1950 if (!affine_hull)
1951 goto error;
1952 if (affine_hull->n_eq != 0)
1953 bset = modulo_affine_hull(set, affine_hull);
1954 else {
1955 isl_basic_set_free(affine_hull);
1956 bset = uset_convex_hull(set);
1959 convex_hull = isl_basic_map_overlying_set(bset, model);
1960 if (!convex_hull)
1961 return NULL;
1963 ISL_F_SET(convex_hull, ISL_BASIC_MAP_NO_IMPLICIT);
1964 ISL_F_SET(convex_hull, ISL_BASIC_MAP_ALL_EQUALITIES);
1965 ISL_F_CLR(convex_hull, ISL_BASIC_MAP_RATIONAL);
1966 return convex_hull;
1967 error:
1968 isl_set_free(set);
1969 isl_basic_map_free(model);
1970 return NULL;
1973 struct isl_basic_set *isl_set_convex_hull(struct isl_set *set)
1975 return (struct isl_basic_set *)
1976 isl_map_convex_hull((struct isl_map *)set);
1979 __isl_give isl_basic_map *isl_map_polyhedral_hull(__isl_take isl_map *map)
1981 isl_basic_map *hull;
1983 hull = isl_map_convex_hull(map);
1984 return isl_basic_map_remove_divs(hull);
1987 __isl_give isl_basic_set *isl_set_polyhedral_hull(__isl_take isl_set *set)
1989 return (isl_basic_set *)isl_map_polyhedral_hull((isl_map *)set);
1992 struct sh_data_entry {
1993 struct isl_hash_table *table;
1994 struct isl_tab *tab;
1997 /* Holds the data needed during the simple hull computation.
1998 * In particular,
1999 * n the number of basic sets in the original set
2000 * hull_table a hash table of already computed constraints
2001 * in the simple hull
2002 * p for each basic set,
2003 * table a hash table of the constraints
2004 * tab the tableau corresponding to the basic set
2006 struct sh_data {
2007 struct isl_ctx *ctx;
2008 unsigned n;
2009 struct isl_hash_table *hull_table;
2010 struct sh_data_entry p[1];
2013 static void sh_data_free(struct sh_data *data)
2015 int i;
2017 if (!data)
2018 return;
2019 isl_hash_table_free(data->ctx, data->hull_table);
2020 for (i = 0; i < data->n; ++i) {
2021 isl_hash_table_free(data->ctx, data->p[i].table);
2022 isl_tab_free(data->p[i].tab);
2024 free(data);
2027 struct ineq_cmp_data {
2028 unsigned len;
2029 isl_int *p;
2032 static int has_ineq(const void *entry, const void *val)
2034 isl_int *row = (isl_int *)entry;
2035 struct ineq_cmp_data *v = (struct ineq_cmp_data *)val;
2037 return isl_seq_eq(row + 1, v->p + 1, v->len) ||
2038 isl_seq_is_neg(row + 1, v->p + 1, v->len);
2041 static int hash_ineq(struct isl_ctx *ctx, struct isl_hash_table *table,
2042 isl_int *ineq, unsigned len)
2044 uint32_t c_hash;
2045 struct ineq_cmp_data v;
2046 struct isl_hash_table_entry *entry;
2048 v.len = len;
2049 v.p = ineq;
2050 c_hash = isl_seq_get_hash(ineq + 1, len);
2051 entry = isl_hash_table_find(ctx, table, c_hash, has_ineq, &v, 1);
2052 if (!entry)
2053 return - 1;
2054 entry->data = ineq;
2055 return 0;
2058 /* Fill hash table "table" with the constraints of "bset".
2059 * Equalities are added as two inequalities.
2060 * The value in the hash table is a pointer to the (in)equality of "bset".
2062 static int hash_basic_set(struct isl_hash_table *table,
2063 struct isl_basic_set *bset)
2065 int i, j;
2066 unsigned dim = isl_basic_set_total_dim(bset);
2068 for (i = 0; i < bset->n_eq; ++i) {
2069 for (j = 0; j < 2; ++j) {
2070 isl_seq_neg(bset->eq[i], bset->eq[i], 1 + dim);
2071 if (hash_ineq(bset->ctx, table, bset->eq[i], dim) < 0)
2072 return -1;
2075 for (i = 0; i < bset->n_ineq; ++i) {
2076 if (hash_ineq(bset->ctx, table, bset->ineq[i], dim) < 0)
2077 return -1;
2079 return 0;
2082 static struct sh_data *sh_data_alloc(struct isl_set *set, unsigned n_ineq)
2084 struct sh_data *data;
2085 int i;
2087 data = isl_calloc(set->ctx, struct sh_data,
2088 sizeof(struct sh_data) +
2089 (set->n - 1) * sizeof(struct sh_data_entry));
2090 if (!data)
2091 return NULL;
2092 data->ctx = set->ctx;
2093 data->n = set->n;
2094 data->hull_table = isl_hash_table_alloc(set->ctx, n_ineq);
2095 if (!data->hull_table)
2096 goto error;
2097 for (i = 0; i < set->n; ++i) {
2098 data->p[i].table = isl_hash_table_alloc(set->ctx,
2099 2 * set->p[i]->n_eq + set->p[i]->n_ineq);
2100 if (!data->p[i].table)
2101 goto error;
2102 if (hash_basic_set(data->p[i].table, set->p[i]) < 0)
2103 goto error;
2105 return data;
2106 error:
2107 sh_data_free(data);
2108 return NULL;
2111 /* Check if inequality "ineq" is a bound for basic set "j" or if
2112 * it can be relaxed (by increasing the constant term) to become
2113 * a bound for that basic set. In the latter case, the constant
2114 * term is updated.
2115 * Relaxation of the constant term is only allowed if "shift" is set.
2117 * Return 1 if "ineq" is a bound
2118 * 0 if "ineq" may attain arbitrarily small values on basic set "j"
2119 * -1 if some error occurred
2121 static int is_bound(struct sh_data *data, struct isl_set *set, int j,
2122 isl_int *ineq, int shift)
2124 enum isl_lp_result res;
2125 isl_int opt;
2127 if (!data->p[j].tab) {
2128 data->p[j].tab = isl_tab_from_basic_set(set->p[j], 0);
2129 if (!data->p[j].tab)
2130 return -1;
2133 isl_int_init(opt);
2135 res = isl_tab_min(data->p[j].tab, ineq, data->ctx->one,
2136 &opt, NULL, 0);
2137 if (res == isl_lp_ok && isl_int_is_neg(opt)) {
2138 if (shift)
2139 isl_int_sub(ineq[0], ineq[0], opt);
2140 else
2141 res = isl_lp_unbounded;
2144 isl_int_clear(opt);
2146 return (res == isl_lp_ok || res == isl_lp_empty) ? 1 :
2147 res == isl_lp_unbounded ? 0 : -1;
2150 /* Check if inequality "ineq" from basic set "i" is or can be relaxed to
2151 * become a bound on the whole set. If so, add the (relaxed) inequality
2152 * to "hull". Relaxation is only allowed if "shift" is set.
2154 * We first check if "hull" already contains a translate of the inequality.
2155 * If so, we are done.
2156 * Then, we check if any of the previous basic sets contains a translate
2157 * of the inequality. If so, then we have already considered this
2158 * inequality and we are done.
2159 * Otherwise, for each basic set other than "i", we check if the inequality
2160 * is a bound on the basic set.
2161 * For previous basic sets, we know that they do not contain a translate
2162 * of the inequality, so we directly call is_bound.
2163 * For following basic sets, we first check if a translate of the
2164 * inequality appears in its description and if so directly update
2165 * the inequality accordingly.
2167 static struct isl_basic_set *add_bound(struct isl_basic_set *hull,
2168 struct sh_data *data, struct isl_set *set, int i, isl_int *ineq,
2169 int shift)
2171 uint32_t c_hash;
2172 struct ineq_cmp_data v;
2173 struct isl_hash_table_entry *entry;
2174 int j, k;
2176 if (!hull)
2177 return NULL;
2179 v.len = isl_basic_set_total_dim(hull);
2180 v.p = ineq;
2181 c_hash = isl_seq_get_hash(ineq + 1, v.len);
2183 entry = isl_hash_table_find(hull->ctx, data->hull_table, c_hash,
2184 has_ineq, &v, 0);
2185 if (entry)
2186 return hull;
2188 for (j = 0; j < i; ++j) {
2189 entry = isl_hash_table_find(hull->ctx, data->p[j].table,
2190 c_hash, has_ineq, &v, 0);
2191 if (entry)
2192 break;
2194 if (j < i)
2195 return hull;
2197 k = isl_basic_set_alloc_inequality(hull);
2198 if (k < 0)
2199 goto error;
2200 isl_seq_cpy(hull->ineq[k], ineq, 1 + v.len);
2202 for (j = 0; j < i; ++j) {
2203 int bound;
2204 bound = is_bound(data, set, j, hull->ineq[k], shift);
2205 if (bound < 0)
2206 goto error;
2207 if (!bound)
2208 break;
2210 if (j < i) {
2211 isl_basic_set_free_inequality(hull, 1);
2212 return hull;
2215 for (j = i + 1; j < set->n; ++j) {
2216 int bound, neg;
2217 isl_int *ineq_j;
2218 entry = isl_hash_table_find(hull->ctx, data->p[j].table,
2219 c_hash, has_ineq, &v, 0);
2220 if (entry) {
2221 ineq_j = entry->data;
2222 neg = isl_seq_is_neg(ineq_j + 1,
2223 hull->ineq[k] + 1, v.len);
2224 if (neg)
2225 isl_int_neg(ineq_j[0], ineq_j[0]);
2226 if (isl_int_gt(ineq_j[0], hull->ineq[k][0]))
2227 isl_int_set(hull->ineq[k][0], ineq_j[0]);
2228 if (neg)
2229 isl_int_neg(ineq_j[0], ineq_j[0]);
2230 continue;
2232 bound = is_bound(data, set, j, hull->ineq[k], shift);
2233 if (bound < 0)
2234 goto error;
2235 if (!bound)
2236 break;
2238 if (j < set->n) {
2239 isl_basic_set_free_inequality(hull, 1);
2240 return hull;
2243 entry = isl_hash_table_find(hull->ctx, data->hull_table, c_hash,
2244 has_ineq, &v, 1);
2245 if (!entry)
2246 goto error;
2247 entry->data = hull->ineq[k];
2249 return hull;
2250 error:
2251 isl_basic_set_free(hull);
2252 return NULL;
2255 /* Check if any inequality from basic set "i" is or can be relaxed to
2256 * become a bound on the whole set. If so, add the (relaxed) inequality
2257 * to "hull". Relaxation is only allowed if "shift" is set.
2259 static struct isl_basic_set *add_bounds(struct isl_basic_set *bset,
2260 struct sh_data *data, struct isl_set *set, int i, int shift)
2262 int j, k;
2263 unsigned dim = isl_basic_set_total_dim(bset);
2265 for (j = 0; j < set->p[i]->n_eq; ++j) {
2266 for (k = 0; k < 2; ++k) {
2267 isl_seq_neg(set->p[i]->eq[j], set->p[i]->eq[j], 1+dim);
2268 bset = add_bound(bset, data, set, i, set->p[i]->eq[j],
2269 shift);
2272 for (j = 0; j < set->p[i]->n_ineq; ++j)
2273 bset = add_bound(bset, data, set, i, set->p[i]->ineq[j], shift);
2274 return bset;
2277 /* Compute a superset of the convex hull of set that is described
2278 * by only (translates of) the constraints in the constituents of set.
2279 * Translation is only allowed if "shift" is set.
2281 static __isl_give isl_basic_set *uset_simple_hull(__isl_take isl_set *set,
2282 int shift)
2284 struct sh_data *data = NULL;
2285 struct isl_basic_set *hull = NULL;
2286 unsigned n_ineq;
2287 int i;
2289 if (!set)
2290 return NULL;
2292 n_ineq = 0;
2293 for (i = 0; i < set->n; ++i) {
2294 if (!set->p[i])
2295 goto error;
2296 n_ineq += 2 * set->p[i]->n_eq + set->p[i]->n_ineq;
2299 hull = isl_basic_set_alloc_space(isl_space_copy(set->dim), 0, 0, n_ineq);
2300 if (!hull)
2301 goto error;
2303 data = sh_data_alloc(set, n_ineq);
2304 if (!data)
2305 goto error;
2307 for (i = 0; i < set->n; ++i)
2308 hull = add_bounds(hull, data, set, i, shift);
2310 sh_data_free(data);
2311 isl_set_free(set);
2313 return hull;
2314 error:
2315 sh_data_free(data);
2316 isl_basic_set_free(hull);
2317 isl_set_free(set);
2318 return NULL;
2321 /* Compute a superset of the convex hull of map that is described
2322 * by only (translates of) the constraints in the constituents of map.
2323 * Translation is only allowed if "shift" is set.
2325 static __isl_give isl_basic_map *map_simple_hull(__isl_take isl_map *map,
2326 int shift)
2328 struct isl_set *set = NULL;
2329 struct isl_basic_map *model = NULL;
2330 struct isl_basic_map *hull;
2331 struct isl_basic_map *affine_hull;
2332 struct isl_basic_set *bset = NULL;
2334 if (!map)
2335 return NULL;
2336 if (map->n == 0) {
2337 hull = isl_basic_map_empty_like_map(map);
2338 isl_map_free(map);
2339 return hull;
2341 if (map->n == 1) {
2342 hull = isl_basic_map_copy(map->p[0]);
2343 isl_map_free(map);
2344 return hull;
2347 map = isl_map_detect_equalities(map);
2348 affine_hull = isl_map_affine_hull(isl_map_copy(map));
2349 map = isl_map_align_divs(map);
2350 model = map ? isl_basic_map_copy(map->p[0]) : NULL;
2352 set = isl_map_underlying_set(map);
2354 bset = uset_simple_hull(set, shift);
2356 hull = isl_basic_map_overlying_set(bset, model);
2358 hull = isl_basic_map_intersect(hull, affine_hull);
2359 hull = isl_basic_map_remove_redundancies(hull);
2361 if (!hull)
2362 return NULL;
2363 ISL_F_SET(hull, ISL_BASIC_MAP_NO_IMPLICIT);
2364 ISL_F_SET(hull, ISL_BASIC_MAP_ALL_EQUALITIES);
2366 hull = isl_basic_map_finalize(hull);
2368 return hull;
2371 /* Compute a superset of the convex hull of map that is described
2372 * by only translates of the constraints in the constituents of map.
2374 __isl_give isl_basic_map *isl_map_simple_hull(__isl_take isl_map *map)
2376 return map_simple_hull(map, 1);
2379 struct isl_basic_set *isl_set_simple_hull(struct isl_set *set)
2381 return (struct isl_basic_set *)
2382 isl_map_simple_hull((struct isl_map *)set);
2385 /* Compute a superset of the convex hull of map that is described
2386 * by only the constraints in the constituents of map.
2388 __isl_give isl_basic_map *isl_map_unshifted_simple_hull(
2389 __isl_take isl_map *map)
2391 return map_simple_hull(map, 0);
2394 __isl_give isl_basic_set *isl_set_unshifted_simple_hull(
2395 __isl_take isl_set *set)
2397 return isl_map_unshifted_simple_hull(set);
2400 /* Check if "ineq" is a bound on "set" and, if so, add it to "hull".
2402 * For each basic set in "set", we first check if the basic set
2403 * contains a translate of "ineq". If this translate is more relaxed,
2404 * then we assume that "ineq" is not a bound on this basic set.
2405 * Otherwise, we know that it is a bound.
2406 * If the basic set does not contain a translate of "ineq", then
2407 * we call is_bound to perform the test.
2409 static __isl_give isl_basic_set *add_bound_from_constraint(
2410 __isl_take isl_basic_set *hull, struct sh_data *data,
2411 __isl_keep isl_set *set, isl_int *ineq)
2413 int i, k;
2414 isl_ctx *ctx;
2415 uint32_t c_hash;
2416 struct ineq_cmp_data v;
2418 if (!hull || !set)
2419 return isl_basic_set_free(hull);
2421 v.len = isl_basic_set_total_dim(hull);
2422 v.p = ineq;
2423 c_hash = isl_seq_get_hash(ineq + 1, v.len);
2425 ctx = isl_basic_set_get_ctx(hull);
2426 for (i = 0; i < set->n; ++i) {
2427 int bound;
2428 struct isl_hash_table_entry *entry;
2430 entry = isl_hash_table_find(ctx, data->p[i].table,
2431 c_hash, &has_ineq, &v, 0);
2432 if (entry) {
2433 isl_int *ineq_i = entry->data;
2434 int neg, more_relaxed;
2436 neg = isl_seq_is_neg(ineq_i + 1, ineq + 1, v.len);
2437 if (neg)
2438 isl_int_neg(ineq_i[0], ineq_i[0]);
2439 more_relaxed = isl_int_gt(ineq_i[0], ineq[0]);
2440 if (neg)
2441 isl_int_neg(ineq_i[0], ineq_i[0]);
2442 if (more_relaxed)
2443 break;
2444 else
2445 continue;
2447 bound = is_bound(data, set, i, ineq, 0);
2448 if (bound < 0)
2449 return isl_basic_set_free(hull);
2450 if (!bound)
2451 break;
2453 if (i < set->n)
2454 return hull;
2456 k = isl_basic_set_alloc_inequality(hull);
2457 if (k < 0)
2458 return isl_basic_set_free(hull);
2459 isl_seq_cpy(hull->ineq[k], ineq, 1 + v.len);
2461 return hull;
2464 /* Compute a superset of the convex hull of "set" that is described
2465 * by only some of the "n_ineq" constraints in the list "ineq", where "set"
2466 * has no parameters or integer divisions.
2468 * The inequalities in "ineq" are assumed to have been sorted such
2469 * that constraints with the same linear part appear together and
2470 * that among constraints with the same linear part, those with
2471 * smaller constant term appear first.
2473 * We reuse the same data structure that is used by uset_simple_hull,
2474 * but we do not need the hull table since we will not consider the
2475 * same constraint more than once. We therefore allocate it with zero size.
2477 * We run through the constraints and try to add them one by one,
2478 * skipping identical constraints. If we have added a constraint and
2479 * the next constraint is a more relaxed translate, then we skip this
2480 * next constraint as well.
2482 static __isl_give isl_basic_set *uset_unshifted_simple_hull_from_constraints(
2483 __isl_take isl_set *set, int n_ineq, isl_int **ineq)
2485 int i;
2486 int last_added = 0;
2487 struct sh_data *data = NULL;
2488 isl_basic_set *hull = NULL;
2489 unsigned dim;
2491 hull = isl_basic_set_alloc_space(isl_set_get_space(set), 0, 0, n_ineq);
2492 if (!hull)
2493 goto error;
2495 data = sh_data_alloc(set, 0);
2496 if (!data)
2497 goto error;
2499 dim = isl_set_dim(set, isl_dim_set);
2500 for (i = 0; i < n_ineq; ++i) {
2501 int hull_n_ineq = hull->n_ineq;
2502 int parallel;
2504 parallel = i > 0 && isl_seq_eq(ineq[i - 1] + 1, ineq[i] + 1,
2505 dim);
2506 if (parallel &&
2507 (last_added || isl_int_eq(ineq[i - 1][0], ineq[i][0])))
2508 continue;
2509 hull = add_bound_from_constraint(hull, data, set, ineq[i]);
2510 if (!hull)
2511 goto error;
2512 last_added = hull->n_ineq > hull_n_ineq;
2515 sh_data_free(data);
2516 isl_set_free(set);
2517 return hull;
2518 error:
2519 sh_data_free(data);
2520 isl_set_free(set);
2521 isl_basic_set_free(hull);
2522 return NULL;
2525 /* Collect pointers to all the inequalities in the elements of "list"
2526 * in "ineq". For equalities, store both a pointer to the equality and
2527 * a pointer to its opposite, which is first copied to "mat".
2528 * "ineq" and "mat" are assumed to have been preallocated to the right size
2529 * (the number of inequalities + 2 times the number of equalites and
2530 * the number of equalities, respectively).
2532 static __isl_give isl_mat *collect_inequalities(__isl_take isl_mat *mat,
2533 __isl_keep isl_basic_set_list *list, isl_int **ineq)
2535 int i, j, n, n_eq, n_ineq;
2537 if (!mat)
2538 return NULL;
2540 n_eq = 0;
2541 n_ineq = 0;
2542 n = isl_basic_set_list_n_basic_set(list);
2543 for (i = 0; i < n; ++i) {
2544 isl_basic_set *bset;
2545 bset = isl_basic_set_list_get_basic_set(list, i);
2546 if (!bset)
2547 return isl_mat_free(mat);
2548 for (j = 0; j < bset->n_eq; ++j) {
2549 ineq[n_ineq++] = mat->row[n_eq];
2550 ineq[n_ineq++] = bset->eq[j];
2551 isl_seq_neg(mat->row[n_eq++], bset->eq[j], mat->n_col);
2553 for (j = 0; j < bset->n_ineq; ++j)
2554 ineq[n_ineq++] = bset->ineq[j];
2555 isl_basic_set_free(bset);
2558 return mat;
2561 /* Comparison routine for use as an isl_sort callback.
2563 * Constraints with the same linear part are sorted together and
2564 * among constraints with the same linear part, those with smaller
2565 * constant term are sorted first.
2567 static int cmp_ineq(const void *a, const void *b, void *arg)
2569 unsigned dim = *(unsigned *) arg;
2570 isl_int * const *ineq1 = a;
2571 isl_int * const *ineq2 = b;
2572 int cmp;
2574 cmp = isl_seq_cmp((*ineq1) + 1, (*ineq2) + 1, dim);
2575 if (cmp != 0)
2576 return cmp;
2577 return isl_int_cmp((*ineq1)[0], (*ineq2)[0]);
2580 /* Compute a superset of the convex hull of "set" that is described
2581 * by only constraints in the elements of "list", where "set" has
2582 * no parameters or integer divisions.
2584 * We collect all the constraints in those elements and then
2585 * sort the constraints such that constraints with the same linear part
2586 * are sorted together and that those with smaller constant term are
2587 * sorted first.
2589 static __isl_give isl_basic_set *uset_unshifted_simple_hull_from_basic_set_list(
2590 __isl_take isl_set *set, __isl_take isl_basic_set_list *list)
2592 int i, n, n_eq, n_ineq;
2593 unsigned dim;
2594 isl_ctx *ctx;
2595 isl_mat *mat = NULL;
2596 isl_int **ineq = NULL;
2597 isl_basic_set *hull;
2599 if (!set)
2600 goto error;
2601 ctx = isl_set_get_ctx(set);
2603 n_eq = 0;
2604 n_ineq = 0;
2605 n = isl_basic_set_list_n_basic_set(list);
2606 for (i = 0; i < n; ++i) {
2607 isl_basic_set *bset;
2608 bset = isl_basic_set_list_get_basic_set(list, i);
2609 if (!bset)
2610 goto error;
2611 n_eq += bset->n_eq;
2612 n_ineq += 2 * bset->n_eq + bset->n_ineq;
2613 isl_basic_set_free(bset);
2616 ineq = isl_alloc_array(ctx, isl_int *, n_ineq);
2617 if (n_ineq > 0 && !ineq)
2618 goto error;
2620 dim = isl_set_dim(set, isl_dim_set);
2621 mat = isl_mat_alloc(ctx, n_eq, 1 + dim);
2622 mat = collect_inequalities(mat, list, ineq);
2623 if (!mat)
2624 goto error;
2626 if (isl_sort(ineq, n_ineq, sizeof(ineq[0]), &cmp_ineq, &dim) < 0)
2627 goto error;
2629 hull = uset_unshifted_simple_hull_from_constraints(set, n_ineq, ineq);
2631 isl_mat_free(mat);
2632 free(ineq);
2633 isl_basic_set_list_free(list);
2634 return hull;
2635 error:
2636 isl_mat_free(mat);
2637 free(ineq);
2638 isl_set_free(set);
2639 isl_basic_set_list_free(list);
2640 return NULL;
2643 /* Compute a superset of the convex hull of "set" that is described
2644 * by only constraints in the elements of "list".
2646 * If the list is empty, then we can only describe the universe set.
2647 * If the input set is empty, then all constraints are valid, so
2648 * we return the intersection of the elements in "list".
2650 * Otherwise, we align all divs and temporarily treat them
2651 * as regular variables, computing the unshifted simple hull in
2652 * uset_unshifted_simple_hull_from_basic_set_list.
2654 static __isl_give isl_basic_set *set_unshifted_simple_hull_from_basic_set_list(
2655 __isl_take isl_set *set, __isl_take isl_basic_set_list *list)
2657 isl_basic_set *model;
2658 isl_basic_set *hull;
2660 if (!set || !list)
2661 goto error;
2663 if (isl_basic_set_list_n_basic_set(list) == 0) {
2664 isl_space *space;
2666 space = isl_set_get_space(set);
2667 isl_set_free(set);
2668 isl_basic_set_list_free(list);
2669 return isl_basic_set_universe(space);
2671 if (isl_set_plain_is_empty(set)) {
2672 isl_set_free(set);
2673 return isl_basic_set_list_intersect(list);
2676 set = isl_set_align_divs_to_basic_set_list(set, list);
2677 if (!set)
2678 goto error;
2679 list = isl_basic_set_list_align_divs_to_basic_set(list, set->p[0]);
2681 model = isl_basic_set_list_get_basic_set(list, 0);
2683 set = isl_set_to_underlying_set(set);
2684 list = isl_basic_set_list_underlying_set(list);
2686 hull = uset_unshifted_simple_hull_from_basic_set_list(set, list);
2687 hull = isl_basic_map_overlying_set(hull, model);
2689 return hull;
2690 error:
2691 isl_set_free(set);
2692 isl_basic_set_list_free(list);
2693 return NULL;
2696 /* Return a sequence of the basic sets that make up the sets in "list".
2698 static __isl_give isl_basic_set_list *collect_basic_sets(
2699 __isl_take isl_set_list *list)
2701 int i, n;
2702 isl_ctx *ctx;
2703 isl_basic_set_list *bset_list;
2705 if (!list)
2706 return NULL;
2707 n = isl_set_list_n_set(list);
2708 ctx = isl_set_list_get_ctx(list);
2709 bset_list = isl_basic_set_list_alloc(ctx, 0);
2711 for (i = 0; i < n; ++i) {
2712 isl_set *set;
2713 isl_basic_set_list *list_i;
2715 set = isl_set_list_get_set(list, i);
2716 set = isl_set_compute_divs(set);
2717 list_i = isl_set_get_basic_set_list(set);
2718 isl_set_free(set);
2719 bset_list = isl_basic_set_list_concat(bset_list, list_i);
2722 isl_set_list_free(list);
2723 return bset_list;
2726 /* Compute a superset of the convex hull of "set" that is described
2727 * by only constraints in the elements of "list".
2729 * If "set" is the universe, then the convex hull (and therefore
2730 * any superset of the convexhull) is the universe as well.
2732 * Otherwise, we collect all the basic sets in the set list and
2733 * continue with set_unshifted_simple_hull_from_basic_set_list.
2735 __isl_give isl_basic_set *isl_set_unshifted_simple_hull_from_set_list(
2736 __isl_take isl_set *set, __isl_take isl_set_list *list)
2738 isl_basic_set_list *bset_list;
2739 int is_universe;
2741 is_universe = isl_set_plain_is_universe(set);
2742 if (is_universe < 0)
2743 set = isl_set_free(set);
2744 if (is_universe < 0 || is_universe) {
2745 isl_set_list_free(list);
2746 return isl_set_unshifted_simple_hull(set);
2749 bset_list = collect_basic_sets(list);
2750 return set_unshifted_simple_hull_from_basic_set_list(set, bset_list);
2753 /* Given a set "set", return parametric bounds on the dimension "dim".
2755 static struct isl_basic_set *set_bounds(struct isl_set *set, int dim)
2757 unsigned set_dim = isl_set_dim(set, isl_dim_set);
2758 set = isl_set_copy(set);
2759 set = isl_set_eliminate_dims(set, dim + 1, set_dim - (dim + 1));
2760 set = isl_set_eliminate_dims(set, 0, dim);
2761 return isl_set_convex_hull(set);
2764 /* Computes a "simple hull" and then check if each dimension in the
2765 * resulting hull is bounded by a symbolic constant. If not, the
2766 * hull is intersected with the corresponding bounds on the whole set.
2768 struct isl_basic_set *isl_set_bounded_simple_hull(struct isl_set *set)
2770 int i, j;
2771 struct isl_basic_set *hull;
2772 unsigned nparam, left;
2773 int removed_divs = 0;
2775 hull = isl_set_simple_hull(isl_set_copy(set));
2776 if (!hull)
2777 goto error;
2779 nparam = isl_basic_set_dim(hull, isl_dim_param);
2780 for (i = 0; i < isl_basic_set_dim(hull, isl_dim_set); ++i) {
2781 int lower = 0, upper = 0;
2782 struct isl_basic_set *bounds;
2784 left = isl_basic_set_total_dim(hull) - nparam - i - 1;
2785 for (j = 0; j < hull->n_eq; ++j) {
2786 if (isl_int_is_zero(hull->eq[j][1 + nparam + i]))
2787 continue;
2788 if (isl_seq_first_non_zero(hull->eq[j]+1+nparam+i+1,
2789 left) == -1)
2790 break;
2792 if (j < hull->n_eq)
2793 continue;
2795 for (j = 0; j < hull->n_ineq; ++j) {
2796 if (isl_int_is_zero(hull->ineq[j][1 + nparam + i]))
2797 continue;
2798 if (isl_seq_first_non_zero(hull->ineq[j]+1+nparam+i+1,
2799 left) != -1 ||
2800 isl_seq_first_non_zero(hull->ineq[j]+1+nparam,
2801 i) != -1)
2802 continue;
2803 if (isl_int_is_pos(hull->ineq[j][1 + nparam + i]))
2804 lower = 1;
2805 else
2806 upper = 1;
2807 if (lower && upper)
2808 break;
2811 if (lower && upper)
2812 continue;
2814 if (!removed_divs) {
2815 set = isl_set_remove_divs(set);
2816 if (!set)
2817 goto error;
2818 removed_divs = 1;
2820 bounds = set_bounds(set, i);
2821 hull = isl_basic_set_intersect(hull, bounds);
2822 if (!hull)
2823 goto error;
2826 isl_set_free(set);
2827 return hull;
2828 error:
2829 isl_set_free(set);
2830 return NULL;