isl_map_union_disjoint: check spaces before handling special cases
[isl.git] / isl_convex_hull.c
blob316b8dca6355c59c011be0dceb7823355abd8c78
1 /*
2 * Copyright 2008-2009 Katholieke Universiteit Leuven
4 * Use of this software is governed by the MIT license
6 * Written by Sven Verdoolaege, K.U.Leuven, Departement
7 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
8 */
10 #include <isl_ctx_private.h>
11 #include <isl_map_private.h>
12 #include <isl_lp_private.h>
13 #include <isl/map.h>
14 #include <isl_mat_private.h>
15 #include <isl_vec_private.h>
16 #include <isl/set.h>
17 #include <isl_seq.h>
18 #include <isl_options_private.h>
19 #include "isl_equalities.h"
20 #include "isl_tab.h"
22 static struct isl_basic_set *uset_convex_hull_wrap_bounded(struct isl_set *set);
24 /* Return 1 if constraint c is redundant with respect to the constraints
25 * in bmap. If c is a lower [upper] bound in some variable and bmap
26 * does not have a lower [upper] bound in that variable, then c cannot
27 * be redundant and we do not need solve any lp.
29 int isl_basic_map_constraint_is_redundant(struct isl_basic_map **bmap,
30 isl_int *c, isl_int *opt_n, isl_int *opt_d)
32 enum isl_lp_result res;
33 unsigned total;
34 int i, j;
36 if (!bmap)
37 return -1;
39 total = isl_basic_map_total_dim(*bmap);
40 for (i = 0; i < total; ++i) {
41 int sign;
42 if (isl_int_is_zero(c[1+i]))
43 continue;
44 sign = isl_int_sgn(c[1+i]);
45 for (j = 0; j < (*bmap)->n_ineq; ++j)
46 if (sign == isl_int_sgn((*bmap)->ineq[j][1+i]))
47 break;
48 if (j == (*bmap)->n_ineq)
49 break;
51 if (i < total)
52 return 0;
54 res = isl_basic_map_solve_lp(*bmap, 0, c, (*bmap)->ctx->one,
55 opt_n, opt_d, NULL);
56 if (res == isl_lp_unbounded)
57 return 0;
58 if (res == isl_lp_error)
59 return -1;
60 if (res == isl_lp_empty) {
61 *bmap = isl_basic_map_set_to_empty(*bmap);
62 return 0;
64 return !isl_int_is_neg(*opt_n);
67 int isl_basic_set_constraint_is_redundant(struct isl_basic_set **bset,
68 isl_int *c, isl_int *opt_n, isl_int *opt_d)
70 return isl_basic_map_constraint_is_redundant(
71 (struct isl_basic_map **)bset, c, opt_n, opt_d);
74 /* Remove redundant
75 * constraints. If the minimal value along the normal of a constraint
76 * is the same if the constraint is removed, then the constraint is redundant.
78 * Alternatively, we could have intersected the basic map with the
79 * corresponding equality and the checked if the dimension was that
80 * of a facet.
82 __isl_give isl_basic_map *isl_basic_map_remove_redundancies(
83 __isl_take isl_basic_map *bmap)
85 struct isl_tab *tab;
87 if (!bmap)
88 return NULL;
90 bmap = isl_basic_map_gauss(bmap, NULL);
91 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
92 return bmap;
93 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_NO_REDUNDANT))
94 return bmap;
95 if (bmap->n_ineq <= 1)
96 return bmap;
98 tab = isl_tab_from_basic_map(bmap, 0);
99 if (isl_tab_detect_implicit_equalities(tab) < 0)
100 goto error;
101 if (isl_tab_detect_redundant(tab) < 0)
102 goto error;
103 bmap = isl_basic_map_update_from_tab(bmap, tab);
104 isl_tab_free(tab);
105 ISL_F_SET(bmap, ISL_BASIC_MAP_NO_IMPLICIT);
106 ISL_F_SET(bmap, ISL_BASIC_MAP_NO_REDUNDANT);
107 return bmap;
108 error:
109 isl_tab_free(tab);
110 isl_basic_map_free(bmap);
111 return NULL;
114 __isl_give isl_basic_set *isl_basic_set_remove_redundancies(
115 __isl_take isl_basic_set *bset)
117 return (struct isl_basic_set *)
118 isl_basic_map_remove_redundancies((struct isl_basic_map *)bset);
121 /* Remove redundant constraints in each of the basic maps.
123 __isl_give isl_map *isl_map_remove_redundancies(__isl_take isl_map *map)
125 return isl_map_inline_foreach_basic_map(map,
126 &isl_basic_map_remove_redundancies);
129 __isl_give isl_set *isl_set_remove_redundancies(__isl_take isl_set *set)
131 return isl_map_remove_redundancies(set);
134 /* Check if the set set is bound in the direction of the affine
135 * constraint c and if so, set the constant term such that the
136 * resulting constraint is a bounding constraint for the set.
138 static int uset_is_bound(struct isl_set *set, isl_int *c, unsigned len)
140 int first;
141 int j;
142 isl_int opt;
143 isl_int opt_denom;
145 isl_int_init(opt);
146 isl_int_init(opt_denom);
147 first = 1;
148 for (j = 0; j < set->n; ++j) {
149 enum isl_lp_result res;
151 if (ISL_F_ISSET(set->p[j], ISL_BASIC_SET_EMPTY))
152 continue;
154 res = isl_basic_set_solve_lp(set->p[j],
155 0, c, set->ctx->one, &opt, &opt_denom, NULL);
156 if (res == isl_lp_unbounded)
157 break;
158 if (res == isl_lp_error)
159 goto error;
160 if (res == isl_lp_empty) {
161 set->p[j] = isl_basic_set_set_to_empty(set->p[j]);
162 if (!set->p[j])
163 goto error;
164 continue;
166 if (first || isl_int_is_neg(opt)) {
167 if (!isl_int_is_one(opt_denom))
168 isl_seq_scale(c, c, opt_denom, len);
169 isl_int_sub(c[0], c[0], opt);
171 first = 0;
173 isl_int_clear(opt);
174 isl_int_clear(opt_denom);
175 return j >= set->n;
176 error:
177 isl_int_clear(opt);
178 isl_int_clear(opt_denom);
179 return -1;
182 __isl_give isl_basic_map *isl_basic_map_set_rational(
183 __isl_take isl_basic_set *bmap)
185 if (!bmap)
186 return NULL;
188 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
189 return bmap;
191 bmap = isl_basic_map_cow(bmap);
192 if (!bmap)
193 return NULL;
195 ISL_F_SET(bmap, ISL_BASIC_MAP_RATIONAL);
197 return isl_basic_map_finalize(bmap);
200 __isl_give isl_basic_set *isl_basic_set_set_rational(
201 __isl_take isl_basic_set *bset)
203 return isl_basic_map_set_rational(bset);
206 __isl_give isl_map *isl_map_set_rational(__isl_take isl_map *map)
208 int i;
210 map = isl_map_cow(map);
211 if (!map)
212 return NULL;
213 for (i = 0; i < map->n; ++i) {
214 map->p[i] = isl_basic_map_set_rational(map->p[i]);
215 if (!map->p[i])
216 goto error;
218 return map;
219 error:
220 isl_map_free(map);
221 return NULL;
224 __isl_give isl_set *isl_set_set_rational(__isl_take isl_set *set)
226 return isl_map_set_rational(set);
229 static struct isl_basic_set *isl_basic_set_add_equality(
230 struct isl_basic_set *bset, isl_int *c)
232 int i;
233 unsigned dim;
235 if (!bset)
236 return NULL;
238 if (ISL_F_ISSET(bset, ISL_BASIC_SET_EMPTY))
239 return bset;
241 isl_assert(bset->ctx, isl_basic_set_n_param(bset) == 0, goto error);
242 isl_assert(bset->ctx, bset->n_div == 0, goto error);
243 dim = isl_basic_set_n_dim(bset);
244 bset = isl_basic_set_cow(bset);
245 bset = isl_basic_set_extend(bset, 0, dim, 0, 1, 0);
246 i = isl_basic_set_alloc_equality(bset);
247 if (i < 0)
248 goto error;
249 isl_seq_cpy(bset->eq[i], c, 1 + dim);
250 return bset;
251 error:
252 isl_basic_set_free(bset);
253 return NULL;
256 static struct isl_set *isl_set_add_basic_set_equality(struct isl_set *set, isl_int *c)
258 int i;
260 set = isl_set_cow(set);
261 if (!set)
262 return NULL;
263 for (i = 0; i < set->n; ++i) {
264 set->p[i] = isl_basic_set_add_equality(set->p[i], c);
265 if (!set->p[i])
266 goto error;
268 return set;
269 error:
270 isl_set_free(set);
271 return NULL;
274 /* Given a union of basic sets, construct the constraints for wrapping
275 * a facet around one of its ridges.
276 * In particular, if each of n the d-dimensional basic sets i in "set"
277 * contains the origin, satisfies the constraints x_1 >= 0 and x_2 >= 0
278 * and is defined by the constraints
279 * [ 1 ]
280 * A_i [ x ] >= 0
282 * then the resulting set is of dimension n*(1+d) and has as constraints
284 * [ a_i ]
285 * A_i [ x_i ] >= 0
287 * a_i >= 0
289 * \sum_i x_{i,1} = 1
291 static struct isl_basic_set *wrap_constraints(struct isl_set *set)
293 struct isl_basic_set *lp;
294 unsigned n_eq;
295 unsigned n_ineq;
296 int i, j, k;
297 unsigned dim, lp_dim;
299 if (!set)
300 return NULL;
302 dim = 1 + isl_set_n_dim(set);
303 n_eq = 1;
304 n_ineq = set->n;
305 for (i = 0; i < set->n; ++i) {
306 n_eq += set->p[i]->n_eq;
307 n_ineq += set->p[i]->n_ineq;
309 lp = isl_basic_set_alloc(set->ctx, 0, dim * set->n, 0, n_eq, n_ineq);
310 lp = isl_basic_set_set_rational(lp);
311 if (!lp)
312 return NULL;
313 lp_dim = isl_basic_set_n_dim(lp);
314 k = isl_basic_set_alloc_equality(lp);
315 isl_int_set_si(lp->eq[k][0], -1);
316 for (i = 0; i < set->n; ++i) {
317 isl_int_set_si(lp->eq[k][1+dim*i], 0);
318 isl_int_set_si(lp->eq[k][1+dim*i+1], 1);
319 isl_seq_clr(lp->eq[k]+1+dim*i+2, dim-2);
321 for (i = 0; i < set->n; ++i) {
322 k = isl_basic_set_alloc_inequality(lp);
323 isl_seq_clr(lp->ineq[k], 1+lp_dim);
324 isl_int_set_si(lp->ineq[k][1+dim*i], 1);
326 for (j = 0; j < set->p[i]->n_eq; ++j) {
327 k = isl_basic_set_alloc_equality(lp);
328 isl_seq_clr(lp->eq[k], 1+dim*i);
329 isl_seq_cpy(lp->eq[k]+1+dim*i, set->p[i]->eq[j], dim);
330 isl_seq_clr(lp->eq[k]+1+dim*(i+1), dim*(set->n-i-1));
333 for (j = 0; j < set->p[i]->n_ineq; ++j) {
334 k = isl_basic_set_alloc_inequality(lp);
335 isl_seq_clr(lp->ineq[k], 1+dim*i);
336 isl_seq_cpy(lp->ineq[k]+1+dim*i, set->p[i]->ineq[j], dim);
337 isl_seq_clr(lp->ineq[k]+1+dim*(i+1), dim*(set->n-i-1));
340 return lp;
343 /* Given a facet "facet" of the convex hull of "set" and a facet "ridge"
344 * of that facet, compute the other facet of the convex hull that contains
345 * the ridge.
347 * We first transform the set such that the facet constraint becomes
349 * x_1 >= 0
351 * I.e., the facet lies in
353 * x_1 = 0
355 * and on that facet, the constraint that defines the ridge is
357 * x_2 >= 0
359 * (This transformation is not strictly needed, all that is needed is
360 * that the ridge contains the origin.)
362 * Since the ridge contains the origin, the cone of the convex hull
363 * will be of the form
365 * x_1 >= 0
366 * x_2 >= a x_1
368 * with this second constraint defining the new facet.
369 * The constant a is obtained by settting x_1 in the cone of the
370 * convex hull to 1 and minimizing x_2.
371 * Now, each element in the cone of the convex hull is the sum
372 * of elements in the cones of the basic sets.
373 * If a_i is the dilation factor of basic set i, then the problem
374 * we need to solve is
376 * min \sum_i x_{i,2}
377 * st
378 * \sum_i x_{i,1} = 1
379 * a_i >= 0
380 * [ a_i ]
381 * A [ x_i ] >= 0
383 * with
384 * [ 1 ]
385 * A_i [ x_i ] >= 0
387 * the constraints of each (transformed) basic set.
388 * If a = n/d, then the constraint defining the new facet (in the transformed
389 * space) is
391 * -n x_1 + d x_2 >= 0
393 * In the original space, we need to take the same combination of the
394 * corresponding constraints "facet" and "ridge".
396 * If a = -infty = "-1/0", then we just return the original facet constraint.
397 * This means that the facet is unbounded, but has a bounded intersection
398 * with the union of sets.
400 isl_int *isl_set_wrap_facet(__isl_keep isl_set *set,
401 isl_int *facet, isl_int *ridge)
403 int i;
404 isl_ctx *ctx;
405 struct isl_mat *T = NULL;
406 struct isl_basic_set *lp = NULL;
407 struct isl_vec *obj;
408 enum isl_lp_result res;
409 isl_int num, den;
410 unsigned dim;
412 if (!set)
413 return NULL;
414 ctx = set->ctx;
415 set = isl_set_copy(set);
416 set = isl_set_set_rational(set);
418 dim = 1 + isl_set_n_dim(set);
419 T = isl_mat_alloc(ctx, 3, dim);
420 if (!T)
421 goto error;
422 isl_int_set_si(T->row[0][0], 1);
423 isl_seq_clr(T->row[0]+1, dim - 1);
424 isl_seq_cpy(T->row[1], facet, dim);
425 isl_seq_cpy(T->row[2], ridge, dim);
426 T = isl_mat_right_inverse(T);
427 set = isl_set_preimage(set, T);
428 T = NULL;
429 if (!set)
430 goto error;
431 lp = wrap_constraints(set);
432 obj = isl_vec_alloc(ctx, 1 + dim*set->n);
433 if (!obj)
434 goto error;
435 isl_int_set_si(obj->block.data[0], 0);
436 for (i = 0; i < set->n; ++i) {
437 isl_seq_clr(obj->block.data + 1 + dim*i, 2);
438 isl_int_set_si(obj->block.data[1 + dim*i+2], 1);
439 isl_seq_clr(obj->block.data + 1 + dim*i+3, dim-3);
441 isl_int_init(num);
442 isl_int_init(den);
443 res = isl_basic_set_solve_lp(lp, 0,
444 obj->block.data, ctx->one, &num, &den, NULL);
445 if (res == isl_lp_ok) {
446 isl_int_neg(num, num);
447 isl_seq_combine(facet, num, facet, den, ridge, dim);
448 isl_seq_normalize(ctx, facet, dim);
450 isl_int_clear(num);
451 isl_int_clear(den);
452 isl_vec_free(obj);
453 isl_basic_set_free(lp);
454 isl_set_free(set);
455 if (res == isl_lp_error)
456 return NULL;
457 isl_assert(ctx, res == isl_lp_ok || res == isl_lp_unbounded,
458 return NULL);
459 return facet;
460 error:
461 isl_basic_set_free(lp);
462 isl_mat_free(T);
463 isl_set_free(set);
464 return NULL;
467 /* Compute the constraint of a facet of "set".
469 * We first compute the intersection with a bounding constraint
470 * that is orthogonal to one of the coordinate axes.
471 * If the affine hull of this intersection has only one equality,
472 * we have found a facet.
473 * Otherwise, we wrap the current bounding constraint around
474 * one of the equalities of the face (one that is not equal to
475 * the current bounding constraint).
476 * This process continues until we have found a facet.
477 * The dimension of the intersection increases by at least
478 * one on each iteration, so termination is guaranteed.
480 static __isl_give isl_mat *initial_facet_constraint(__isl_keep isl_set *set)
482 struct isl_set *slice = NULL;
483 struct isl_basic_set *face = NULL;
484 int i;
485 unsigned dim = isl_set_n_dim(set);
486 int is_bound;
487 isl_mat *bounds = NULL;
489 isl_assert(set->ctx, set->n > 0, goto error);
490 bounds = isl_mat_alloc(set->ctx, 1, 1 + dim);
491 if (!bounds)
492 return NULL;
494 isl_seq_clr(bounds->row[0], dim);
495 isl_int_set_si(bounds->row[0][1 + dim - 1], 1);
496 is_bound = uset_is_bound(set, bounds->row[0], 1 + dim);
497 if (is_bound < 0)
498 goto error;
499 isl_assert(set->ctx, is_bound, goto error);
500 isl_seq_normalize(set->ctx, bounds->row[0], 1 + dim);
501 bounds->n_row = 1;
503 for (;;) {
504 slice = isl_set_copy(set);
505 slice = isl_set_add_basic_set_equality(slice, bounds->row[0]);
506 face = isl_set_affine_hull(slice);
507 if (!face)
508 goto error;
509 if (face->n_eq == 1) {
510 isl_basic_set_free(face);
511 break;
513 for (i = 0; i < face->n_eq; ++i)
514 if (!isl_seq_eq(bounds->row[0], face->eq[i], 1 + dim) &&
515 !isl_seq_is_neg(bounds->row[0],
516 face->eq[i], 1 + dim))
517 break;
518 isl_assert(set->ctx, i < face->n_eq, goto error);
519 if (!isl_set_wrap_facet(set, bounds->row[0], face->eq[i]))
520 goto error;
521 isl_seq_normalize(set->ctx, bounds->row[0], bounds->n_col);
522 isl_basic_set_free(face);
525 return bounds;
526 error:
527 isl_basic_set_free(face);
528 isl_mat_free(bounds);
529 return NULL;
532 /* Given the bounding constraint "c" of a facet of the convex hull of "set",
533 * compute a hyperplane description of the facet, i.e., compute the facets
534 * of the facet.
536 * We compute an affine transformation that transforms the constraint
538 * [ 1 ]
539 * c [ x ] = 0
541 * to the constraint
543 * z_1 = 0
545 * by computing the right inverse U of a matrix that starts with the rows
547 * [ 1 0 ]
548 * [ c ]
550 * Then
551 * [ 1 ] [ 1 ]
552 * [ x ] = U [ z ]
553 * and
554 * [ 1 ] [ 1 ]
555 * [ z ] = Q [ x ]
557 * with Q = U^{-1}
558 * Since z_1 is zero, we can drop this variable as well as the corresponding
559 * column of U to obtain
561 * [ 1 ] [ 1 ]
562 * [ x ] = U' [ z' ]
563 * and
564 * [ 1 ] [ 1 ]
565 * [ z' ] = Q' [ x ]
567 * with Q' equal to Q, but without the corresponding row.
568 * After computing the facets of the facet in the z' space,
569 * we convert them back to the x space through Q.
571 static struct isl_basic_set *compute_facet(struct isl_set *set, isl_int *c)
573 struct isl_mat *m, *U, *Q;
574 struct isl_basic_set *facet = NULL;
575 struct isl_ctx *ctx;
576 unsigned dim;
578 ctx = set->ctx;
579 set = isl_set_copy(set);
580 dim = isl_set_n_dim(set);
581 m = isl_mat_alloc(set->ctx, 2, 1 + dim);
582 if (!m)
583 goto error;
584 isl_int_set_si(m->row[0][0], 1);
585 isl_seq_clr(m->row[0]+1, dim);
586 isl_seq_cpy(m->row[1], c, 1+dim);
587 U = isl_mat_right_inverse(m);
588 Q = isl_mat_right_inverse(isl_mat_copy(U));
589 U = isl_mat_drop_cols(U, 1, 1);
590 Q = isl_mat_drop_rows(Q, 1, 1);
591 set = isl_set_preimage(set, U);
592 facet = uset_convex_hull_wrap_bounded(set);
593 facet = isl_basic_set_preimage(facet, Q);
594 if (facet)
595 isl_assert(ctx, facet->n_eq == 0, goto error);
596 return facet;
597 error:
598 isl_basic_set_free(facet);
599 isl_set_free(set);
600 return NULL;
603 /* Given an initial facet constraint, compute the remaining facets.
604 * We do this by running through all facets found so far and computing
605 * the adjacent facets through wrapping, adding those facets that we
606 * hadn't already found before.
608 * For each facet we have found so far, we first compute its facets
609 * in the resulting convex hull. That is, we compute the ridges
610 * of the resulting convex hull contained in the facet.
611 * We also compute the corresponding facet in the current approximation
612 * of the convex hull. There is no need to wrap around the ridges
613 * in this facet since that would result in a facet that is already
614 * present in the current approximation.
616 * This function can still be significantly optimized by checking which of
617 * the facets of the basic sets are also facets of the convex hull and
618 * using all the facets so far to help in constructing the facets of the
619 * facets
620 * and/or
621 * using the technique in section "3.1 Ridge Generation" of
622 * "Extended Convex Hull" by Fukuda et al.
624 static struct isl_basic_set *extend(struct isl_basic_set *hull,
625 struct isl_set *set)
627 int i, j, f;
628 int k;
629 struct isl_basic_set *facet = NULL;
630 struct isl_basic_set *hull_facet = NULL;
631 unsigned dim;
633 if (!hull)
634 return NULL;
636 isl_assert(set->ctx, set->n > 0, goto error);
638 dim = isl_set_n_dim(set);
640 for (i = 0; i < hull->n_ineq; ++i) {
641 facet = compute_facet(set, hull->ineq[i]);
642 facet = isl_basic_set_add_equality(facet, hull->ineq[i]);
643 facet = isl_basic_set_gauss(facet, NULL);
644 facet = isl_basic_set_normalize_constraints(facet);
645 hull_facet = isl_basic_set_copy(hull);
646 hull_facet = isl_basic_set_add_equality(hull_facet, hull->ineq[i]);
647 hull_facet = isl_basic_set_gauss(hull_facet, NULL);
648 hull_facet = isl_basic_set_normalize_constraints(hull_facet);
649 if (!facet || !hull_facet)
650 goto error;
651 hull = isl_basic_set_cow(hull);
652 hull = isl_basic_set_extend_space(hull,
653 isl_space_copy(hull->dim), 0, 0, facet->n_ineq);
654 if (!hull)
655 goto error;
656 for (j = 0; j < facet->n_ineq; ++j) {
657 for (f = 0; f < hull_facet->n_ineq; ++f)
658 if (isl_seq_eq(facet->ineq[j],
659 hull_facet->ineq[f], 1 + dim))
660 break;
661 if (f < hull_facet->n_ineq)
662 continue;
663 k = isl_basic_set_alloc_inequality(hull);
664 if (k < 0)
665 goto error;
666 isl_seq_cpy(hull->ineq[k], hull->ineq[i], 1+dim);
667 if (!isl_set_wrap_facet(set, hull->ineq[k], facet->ineq[j]))
668 goto error;
670 isl_basic_set_free(hull_facet);
671 isl_basic_set_free(facet);
673 hull = isl_basic_set_simplify(hull);
674 hull = isl_basic_set_finalize(hull);
675 return hull;
676 error:
677 isl_basic_set_free(hull_facet);
678 isl_basic_set_free(facet);
679 isl_basic_set_free(hull);
680 return NULL;
683 /* Special case for computing the convex hull of a one dimensional set.
684 * We simply collect the lower and upper bounds of each basic set
685 * and the biggest of those.
687 static struct isl_basic_set *convex_hull_1d(struct isl_set *set)
689 struct isl_mat *c = NULL;
690 isl_int *lower = NULL;
691 isl_int *upper = NULL;
692 int i, j, k;
693 isl_int a, b;
694 struct isl_basic_set *hull;
696 for (i = 0; i < set->n; ++i) {
697 set->p[i] = isl_basic_set_simplify(set->p[i]);
698 if (!set->p[i])
699 goto error;
701 set = isl_set_remove_empty_parts(set);
702 if (!set)
703 goto error;
704 isl_assert(set->ctx, set->n > 0, goto error);
705 c = isl_mat_alloc(set->ctx, 2, 2);
706 if (!c)
707 goto error;
709 if (set->p[0]->n_eq > 0) {
710 isl_assert(set->ctx, set->p[0]->n_eq == 1, goto error);
711 lower = c->row[0];
712 upper = c->row[1];
713 if (isl_int_is_pos(set->p[0]->eq[0][1])) {
714 isl_seq_cpy(lower, set->p[0]->eq[0], 2);
715 isl_seq_neg(upper, set->p[0]->eq[0], 2);
716 } else {
717 isl_seq_neg(lower, set->p[0]->eq[0], 2);
718 isl_seq_cpy(upper, set->p[0]->eq[0], 2);
720 } else {
721 for (j = 0; j < set->p[0]->n_ineq; ++j) {
722 if (isl_int_is_pos(set->p[0]->ineq[j][1])) {
723 lower = c->row[0];
724 isl_seq_cpy(lower, set->p[0]->ineq[j], 2);
725 } else {
726 upper = c->row[1];
727 isl_seq_cpy(upper, set->p[0]->ineq[j], 2);
732 isl_int_init(a);
733 isl_int_init(b);
734 for (i = 0; i < set->n; ++i) {
735 struct isl_basic_set *bset = set->p[i];
736 int has_lower = 0;
737 int has_upper = 0;
739 for (j = 0; j < bset->n_eq; ++j) {
740 has_lower = 1;
741 has_upper = 1;
742 if (lower) {
743 isl_int_mul(a, lower[0], bset->eq[j][1]);
744 isl_int_mul(b, lower[1], bset->eq[j][0]);
745 if (isl_int_lt(a, b) && isl_int_is_pos(bset->eq[j][1]))
746 isl_seq_cpy(lower, bset->eq[j], 2);
747 if (isl_int_gt(a, b) && isl_int_is_neg(bset->eq[j][1]))
748 isl_seq_neg(lower, bset->eq[j], 2);
750 if (upper) {
751 isl_int_mul(a, upper[0], bset->eq[j][1]);
752 isl_int_mul(b, upper[1], bset->eq[j][0]);
753 if (isl_int_lt(a, b) && isl_int_is_pos(bset->eq[j][1]))
754 isl_seq_neg(upper, bset->eq[j], 2);
755 if (isl_int_gt(a, b) && isl_int_is_neg(bset->eq[j][1]))
756 isl_seq_cpy(upper, bset->eq[j], 2);
759 for (j = 0; j < bset->n_ineq; ++j) {
760 if (isl_int_is_pos(bset->ineq[j][1]))
761 has_lower = 1;
762 if (isl_int_is_neg(bset->ineq[j][1]))
763 has_upper = 1;
764 if (lower && isl_int_is_pos(bset->ineq[j][1])) {
765 isl_int_mul(a, lower[0], bset->ineq[j][1]);
766 isl_int_mul(b, lower[1], bset->ineq[j][0]);
767 if (isl_int_lt(a, b))
768 isl_seq_cpy(lower, bset->ineq[j], 2);
770 if (upper && isl_int_is_neg(bset->ineq[j][1])) {
771 isl_int_mul(a, upper[0], bset->ineq[j][1]);
772 isl_int_mul(b, upper[1], bset->ineq[j][0]);
773 if (isl_int_gt(a, b))
774 isl_seq_cpy(upper, bset->ineq[j], 2);
777 if (!has_lower)
778 lower = NULL;
779 if (!has_upper)
780 upper = NULL;
782 isl_int_clear(a);
783 isl_int_clear(b);
785 hull = isl_basic_set_alloc(set->ctx, 0, 1, 0, 0, 2);
786 hull = isl_basic_set_set_rational(hull);
787 if (!hull)
788 goto error;
789 if (lower) {
790 k = isl_basic_set_alloc_inequality(hull);
791 isl_seq_cpy(hull->ineq[k], lower, 2);
793 if (upper) {
794 k = isl_basic_set_alloc_inequality(hull);
795 isl_seq_cpy(hull->ineq[k], upper, 2);
797 hull = isl_basic_set_finalize(hull);
798 isl_set_free(set);
799 isl_mat_free(c);
800 return hull;
801 error:
802 isl_set_free(set);
803 isl_mat_free(c);
804 return NULL;
807 static struct isl_basic_set *convex_hull_0d(struct isl_set *set)
809 struct isl_basic_set *convex_hull;
811 if (!set)
812 return NULL;
814 if (isl_set_is_empty(set))
815 convex_hull = isl_basic_set_empty(isl_space_copy(set->dim));
816 else
817 convex_hull = isl_basic_set_universe(isl_space_copy(set->dim));
818 isl_set_free(set);
819 return convex_hull;
822 /* Compute the convex hull of a pair of basic sets without any parameters or
823 * integer divisions using Fourier-Motzkin elimination.
824 * The convex hull is the set of all points that can be written as
825 * the sum of points from both basic sets (in homogeneous coordinates).
826 * We set up the constraints in a space with dimensions for each of
827 * the three sets and then project out the dimensions corresponding
828 * to the two original basic sets, retaining only those corresponding
829 * to the convex hull.
831 static struct isl_basic_set *convex_hull_pair_elim(struct isl_basic_set *bset1,
832 struct isl_basic_set *bset2)
834 int i, j, k;
835 struct isl_basic_set *bset[2];
836 struct isl_basic_set *hull = NULL;
837 unsigned dim;
839 if (!bset1 || !bset2)
840 goto error;
842 dim = isl_basic_set_n_dim(bset1);
843 hull = isl_basic_set_alloc(bset1->ctx, 0, 2 + 3 * dim, 0,
844 1 + dim + bset1->n_eq + bset2->n_eq,
845 2 + bset1->n_ineq + bset2->n_ineq);
846 bset[0] = bset1;
847 bset[1] = bset2;
848 for (i = 0; i < 2; ++i) {
849 for (j = 0; j < bset[i]->n_eq; ++j) {
850 k = isl_basic_set_alloc_equality(hull);
851 if (k < 0)
852 goto error;
853 isl_seq_clr(hull->eq[k], (i+1) * (1+dim));
854 isl_seq_clr(hull->eq[k]+(i+2)*(1+dim), (1-i)*(1+dim));
855 isl_seq_cpy(hull->eq[k]+(i+1)*(1+dim), bset[i]->eq[j],
856 1+dim);
858 for (j = 0; j < bset[i]->n_ineq; ++j) {
859 k = isl_basic_set_alloc_inequality(hull);
860 if (k < 0)
861 goto error;
862 isl_seq_clr(hull->ineq[k], (i+1) * (1+dim));
863 isl_seq_clr(hull->ineq[k]+(i+2)*(1+dim), (1-i)*(1+dim));
864 isl_seq_cpy(hull->ineq[k]+(i+1)*(1+dim),
865 bset[i]->ineq[j], 1+dim);
867 k = isl_basic_set_alloc_inequality(hull);
868 if (k < 0)
869 goto error;
870 isl_seq_clr(hull->ineq[k], 1+2+3*dim);
871 isl_int_set_si(hull->ineq[k][(i+1)*(1+dim)], 1);
873 for (j = 0; j < 1+dim; ++j) {
874 k = isl_basic_set_alloc_equality(hull);
875 if (k < 0)
876 goto error;
877 isl_seq_clr(hull->eq[k], 1+2+3*dim);
878 isl_int_set_si(hull->eq[k][j], -1);
879 isl_int_set_si(hull->eq[k][1+dim+j], 1);
880 isl_int_set_si(hull->eq[k][2*(1+dim)+j], 1);
882 hull = isl_basic_set_set_rational(hull);
883 hull = isl_basic_set_remove_dims(hull, isl_dim_set, dim, 2*(1+dim));
884 hull = isl_basic_set_remove_redundancies(hull);
885 isl_basic_set_free(bset1);
886 isl_basic_set_free(bset2);
887 return hull;
888 error:
889 isl_basic_set_free(bset1);
890 isl_basic_set_free(bset2);
891 isl_basic_set_free(hull);
892 return NULL;
895 /* Is the set bounded for each value of the parameters?
897 int isl_basic_set_is_bounded(__isl_keep isl_basic_set *bset)
899 struct isl_tab *tab;
900 int bounded;
902 if (!bset)
903 return -1;
904 if (isl_basic_set_plain_is_empty(bset))
905 return 1;
907 tab = isl_tab_from_recession_cone(bset, 1);
908 bounded = isl_tab_cone_is_bounded(tab);
909 isl_tab_free(tab);
910 return bounded;
913 /* Is the image bounded for each value of the parameters and
914 * the domain variables?
916 int isl_basic_map_image_is_bounded(__isl_keep isl_basic_map *bmap)
918 unsigned nparam = isl_basic_map_dim(bmap, isl_dim_param);
919 unsigned n_in = isl_basic_map_dim(bmap, isl_dim_in);
920 int bounded;
922 bmap = isl_basic_map_copy(bmap);
923 bmap = isl_basic_map_cow(bmap);
924 bmap = isl_basic_map_move_dims(bmap, isl_dim_param, nparam,
925 isl_dim_in, 0, n_in);
926 bounded = isl_basic_set_is_bounded((isl_basic_set *)bmap);
927 isl_basic_map_free(bmap);
929 return bounded;
932 /* Is the set bounded for each value of the parameters?
934 int isl_set_is_bounded(__isl_keep isl_set *set)
936 int i;
938 if (!set)
939 return -1;
941 for (i = 0; i < set->n; ++i) {
942 int bounded = isl_basic_set_is_bounded(set->p[i]);
943 if (!bounded || bounded < 0)
944 return bounded;
946 return 1;
949 /* Compute the lineality space of the convex hull of bset1 and bset2.
951 * We first compute the intersection of the recession cone of bset1
952 * with the negative of the recession cone of bset2 and then compute
953 * the linear hull of the resulting cone.
955 static struct isl_basic_set *induced_lineality_space(
956 struct isl_basic_set *bset1, struct isl_basic_set *bset2)
958 int i, k;
959 struct isl_basic_set *lin = NULL;
960 unsigned dim;
962 if (!bset1 || !bset2)
963 goto error;
965 dim = isl_basic_set_total_dim(bset1);
966 lin = isl_basic_set_alloc_space(isl_basic_set_get_space(bset1), 0,
967 bset1->n_eq + bset2->n_eq,
968 bset1->n_ineq + bset2->n_ineq);
969 lin = isl_basic_set_set_rational(lin);
970 if (!lin)
971 goto error;
972 for (i = 0; i < bset1->n_eq; ++i) {
973 k = isl_basic_set_alloc_equality(lin);
974 if (k < 0)
975 goto error;
976 isl_int_set_si(lin->eq[k][0], 0);
977 isl_seq_cpy(lin->eq[k] + 1, bset1->eq[i] + 1, dim);
979 for (i = 0; i < bset1->n_ineq; ++i) {
980 k = isl_basic_set_alloc_inequality(lin);
981 if (k < 0)
982 goto error;
983 isl_int_set_si(lin->ineq[k][0], 0);
984 isl_seq_cpy(lin->ineq[k] + 1, bset1->ineq[i] + 1, dim);
986 for (i = 0; i < bset2->n_eq; ++i) {
987 k = isl_basic_set_alloc_equality(lin);
988 if (k < 0)
989 goto error;
990 isl_int_set_si(lin->eq[k][0], 0);
991 isl_seq_neg(lin->eq[k] + 1, bset2->eq[i] + 1, dim);
993 for (i = 0; i < bset2->n_ineq; ++i) {
994 k = isl_basic_set_alloc_inequality(lin);
995 if (k < 0)
996 goto error;
997 isl_int_set_si(lin->ineq[k][0], 0);
998 isl_seq_neg(lin->ineq[k] + 1, bset2->ineq[i] + 1, dim);
1001 isl_basic_set_free(bset1);
1002 isl_basic_set_free(bset2);
1003 return isl_basic_set_affine_hull(lin);
1004 error:
1005 isl_basic_set_free(lin);
1006 isl_basic_set_free(bset1);
1007 isl_basic_set_free(bset2);
1008 return NULL;
1011 static struct isl_basic_set *uset_convex_hull(struct isl_set *set);
1013 /* Given a set and a linear space "lin" of dimension n > 0,
1014 * project the linear space from the set, compute the convex hull
1015 * and then map the set back to the original space.
1017 * Let
1019 * M x = 0
1021 * describe the linear space. We first compute the Hermite normal
1022 * form H = M U of M = H Q, to obtain
1024 * H Q x = 0
1026 * The last n rows of H will be zero, so the last n variables of x' = Q x
1027 * are the one we want to project out. We do this by transforming each
1028 * basic set A x >= b to A U x' >= b and then removing the last n dimensions.
1029 * After computing the convex hull in x'_1, i.e., A' x'_1 >= b',
1030 * we transform the hull back to the original space as A' Q_1 x >= b',
1031 * with Q_1 all but the last n rows of Q.
1033 static struct isl_basic_set *modulo_lineality(struct isl_set *set,
1034 struct isl_basic_set *lin)
1036 unsigned total = isl_basic_set_total_dim(lin);
1037 unsigned lin_dim;
1038 struct isl_basic_set *hull;
1039 struct isl_mat *M, *U, *Q;
1041 if (!set || !lin)
1042 goto error;
1043 lin_dim = total - lin->n_eq;
1044 M = isl_mat_sub_alloc6(set->ctx, lin->eq, 0, lin->n_eq, 1, total);
1045 M = isl_mat_left_hermite(M, 0, &U, &Q);
1046 if (!M)
1047 goto error;
1048 isl_mat_free(M);
1049 isl_basic_set_free(lin);
1051 Q = isl_mat_drop_rows(Q, Q->n_row - lin_dim, lin_dim);
1053 U = isl_mat_lin_to_aff(U);
1054 Q = isl_mat_lin_to_aff(Q);
1056 set = isl_set_preimage(set, U);
1057 set = isl_set_remove_dims(set, isl_dim_set, total - lin_dim, lin_dim);
1058 hull = uset_convex_hull(set);
1059 hull = isl_basic_set_preimage(hull, Q);
1061 return hull;
1062 error:
1063 isl_basic_set_free(lin);
1064 isl_set_free(set);
1065 return NULL;
1068 /* Given two polyhedra with as constraints h_{ij} x >= 0 in homegeneous space,
1069 * set up an LP for solving
1071 * \sum_j \alpha_{1j} h_{1j} = \sum_j \alpha_{2j} h_{2j}
1073 * \alpha{i0} corresponds to the (implicit) positivity constraint 1 >= 0
1074 * The next \alpha{ij} correspond to the equalities and come in pairs.
1075 * The final \alpha{ij} correspond to the inequalities.
1077 static struct isl_basic_set *valid_direction_lp(
1078 struct isl_basic_set *bset1, struct isl_basic_set *bset2)
1080 isl_space *dim;
1081 struct isl_basic_set *lp;
1082 unsigned d;
1083 int n;
1084 int i, j, k;
1086 if (!bset1 || !bset2)
1087 goto error;
1088 d = 1 + isl_basic_set_total_dim(bset1);
1089 n = 2 +
1090 2 * bset1->n_eq + bset1->n_ineq + 2 * bset2->n_eq + bset2->n_ineq;
1091 dim = isl_space_set_alloc(bset1->ctx, 0, n);
1092 lp = isl_basic_set_alloc_space(dim, 0, d, n);
1093 if (!lp)
1094 goto error;
1095 for (i = 0; i < n; ++i) {
1096 k = isl_basic_set_alloc_inequality(lp);
1097 if (k < 0)
1098 goto error;
1099 isl_seq_clr(lp->ineq[k] + 1, n);
1100 isl_int_set_si(lp->ineq[k][0], -1);
1101 isl_int_set_si(lp->ineq[k][1 + i], 1);
1103 for (i = 0; i < d; ++i) {
1104 k = isl_basic_set_alloc_equality(lp);
1105 if (k < 0)
1106 goto error;
1107 n = 0;
1108 isl_int_set_si(lp->eq[k][n], 0); n++;
1109 /* positivity constraint 1 >= 0 */
1110 isl_int_set_si(lp->eq[k][n], i == 0); n++;
1111 for (j = 0; j < bset1->n_eq; ++j) {
1112 isl_int_set(lp->eq[k][n], bset1->eq[j][i]); n++;
1113 isl_int_neg(lp->eq[k][n], bset1->eq[j][i]); n++;
1115 for (j = 0; j < bset1->n_ineq; ++j) {
1116 isl_int_set(lp->eq[k][n], bset1->ineq[j][i]); n++;
1118 /* positivity constraint 1 >= 0 */
1119 isl_int_set_si(lp->eq[k][n], -(i == 0)); n++;
1120 for (j = 0; j < bset2->n_eq; ++j) {
1121 isl_int_neg(lp->eq[k][n], bset2->eq[j][i]); n++;
1122 isl_int_set(lp->eq[k][n], bset2->eq[j][i]); n++;
1124 for (j = 0; j < bset2->n_ineq; ++j) {
1125 isl_int_neg(lp->eq[k][n], bset2->ineq[j][i]); n++;
1128 lp = isl_basic_set_gauss(lp, NULL);
1129 isl_basic_set_free(bset1);
1130 isl_basic_set_free(bset2);
1131 return lp;
1132 error:
1133 isl_basic_set_free(bset1);
1134 isl_basic_set_free(bset2);
1135 return NULL;
1138 /* Compute a vector s in the homogeneous space such that <s, r> > 0
1139 * for all rays in the homogeneous space of the two cones that correspond
1140 * to the input polyhedra bset1 and bset2.
1142 * We compute s as a vector that satisfies
1144 * s = \sum_j \alpha_{ij} h_{ij} for i = 1,2 (*)
1146 * with h_{ij} the normals of the facets of polyhedron i
1147 * (including the "positivity constraint" 1 >= 0) and \alpha_{ij}
1148 * strictly positive numbers. For simplicity we impose \alpha_{ij} >= 1.
1149 * We first set up an LP with as variables the \alpha{ij}.
1150 * In this formulation, for each polyhedron i,
1151 * the first constraint is the positivity constraint, followed by pairs
1152 * of variables for the equalities, followed by variables for the inequalities.
1153 * We then simply pick a feasible solution and compute s using (*).
1155 * Note that we simply pick any valid direction and make no attempt
1156 * to pick a "good" or even the "best" valid direction.
1158 static struct isl_vec *valid_direction(
1159 struct isl_basic_set *bset1, struct isl_basic_set *bset2)
1161 struct isl_basic_set *lp;
1162 struct isl_tab *tab;
1163 struct isl_vec *sample = NULL;
1164 struct isl_vec *dir;
1165 unsigned d;
1166 int i;
1167 int n;
1169 if (!bset1 || !bset2)
1170 goto error;
1171 lp = valid_direction_lp(isl_basic_set_copy(bset1),
1172 isl_basic_set_copy(bset2));
1173 tab = isl_tab_from_basic_set(lp, 0);
1174 sample = isl_tab_get_sample_value(tab);
1175 isl_tab_free(tab);
1176 isl_basic_set_free(lp);
1177 if (!sample)
1178 goto error;
1179 d = isl_basic_set_total_dim(bset1);
1180 dir = isl_vec_alloc(bset1->ctx, 1 + d);
1181 if (!dir)
1182 goto error;
1183 isl_seq_clr(dir->block.data + 1, dir->size - 1);
1184 n = 1;
1185 /* positivity constraint 1 >= 0 */
1186 isl_int_set(dir->block.data[0], sample->block.data[n]); n++;
1187 for (i = 0; i < bset1->n_eq; ++i) {
1188 isl_int_sub(sample->block.data[n],
1189 sample->block.data[n], sample->block.data[n+1]);
1190 isl_seq_combine(dir->block.data,
1191 bset1->ctx->one, dir->block.data,
1192 sample->block.data[n], bset1->eq[i], 1 + d);
1194 n += 2;
1196 for (i = 0; i < bset1->n_ineq; ++i)
1197 isl_seq_combine(dir->block.data,
1198 bset1->ctx->one, dir->block.data,
1199 sample->block.data[n++], bset1->ineq[i], 1 + d);
1200 isl_vec_free(sample);
1201 isl_seq_normalize(bset1->ctx, dir->el, dir->size);
1202 isl_basic_set_free(bset1);
1203 isl_basic_set_free(bset2);
1204 return dir;
1205 error:
1206 isl_vec_free(sample);
1207 isl_basic_set_free(bset1);
1208 isl_basic_set_free(bset2);
1209 return NULL;
1212 /* Given a polyhedron b_i + A_i x >= 0 and a map T = S^{-1},
1213 * compute b_i' + A_i' x' >= 0, with
1215 * [ b_i A_i ] [ y' ] [ y' ]
1216 * [ 1 0 ] S^{-1} [ x' ] >= 0 or [ b_i' A_i' ] [ x' ] >= 0
1218 * In particular, add the "positivity constraint" and then perform
1219 * the mapping.
1221 static struct isl_basic_set *homogeneous_map(struct isl_basic_set *bset,
1222 struct isl_mat *T)
1224 int k;
1226 if (!bset)
1227 goto error;
1228 bset = isl_basic_set_extend_constraints(bset, 0, 1);
1229 k = isl_basic_set_alloc_inequality(bset);
1230 if (k < 0)
1231 goto error;
1232 isl_seq_clr(bset->ineq[k] + 1, isl_basic_set_total_dim(bset));
1233 isl_int_set_si(bset->ineq[k][0], 1);
1234 bset = isl_basic_set_preimage(bset, T);
1235 return bset;
1236 error:
1237 isl_mat_free(T);
1238 isl_basic_set_free(bset);
1239 return NULL;
1242 /* Compute the convex hull of a pair of basic sets without any parameters or
1243 * integer divisions, where the convex hull is known to be pointed,
1244 * but the basic sets may be unbounded.
1246 * We turn this problem into the computation of a convex hull of a pair
1247 * _bounded_ polyhedra by "changing the direction of the homogeneous
1248 * dimension". This idea is due to Matthias Koeppe.
1250 * Consider the cones in homogeneous space that correspond to the
1251 * input polyhedra. The rays of these cones are also rays of the
1252 * polyhedra if the coordinate that corresponds to the homogeneous
1253 * dimension is zero. That is, if the inner product of the rays
1254 * with the homogeneous direction is zero.
1255 * The cones in the homogeneous space can also be considered to
1256 * correspond to other pairs of polyhedra by chosing a different
1257 * homogeneous direction. To ensure that both of these polyhedra
1258 * are bounded, we need to make sure that all rays of the cones
1259 * correspond to vertices and not to rays.
1260 * Let s be a direction such that <s, r> > 0 for all rays r of both cones.
1261 * Then using s as a homogeneous direction, we obtain a pair of polytopes.
1262 * The vector s is computed in valid_direction.
1264 * Note that we need to consider _all_ rays of the cones and not just
1265 * the rays that correspond to rays in the polyhedra. If we were to
1266 * only consider those rays and turn them into vertices, then we
1267 * may inadvertently turn some vertices into rays.
1269 * The standard homogeneous direction is the unit vector in the 0th coordinate.
1270 * We therefore transform the two polyhedra such that the selected
1271 * direction is mapped onto this standard direction and then proceed
1272 * with the normal computation.
1273 * Let S be a non-singular square matrix with s as its first row,
1274 * then we want to map the polyhedra to the space
1276 * [ y' ] [ y ] [ y ] [ y' ]
1277 * [ x' ] = S [ x ] i.e., [ x ] = S^{-1} [ x' ]
1279 * We take S to be the unimodular completion of s to limit the growth
1280 * of the coefficients in the following computations.
1282 * Let b_i + A_i x >= 0 be the constraints of polyhedron i.
1283 * We first move to the homogeneous dimension
1285 * b_i y + A_i x >= 0 [ b_i A_i ] [ y ] [ 0 ]
1286 * y >= 0 or [ 1 0 ] [ x ] >= [ 0 ]
1288 * Then we change directoin
1290 * [ b_i A_i ] [ y' ] [ y' ]
1291 * [ 1 0 ] S^{-1} [ x' ] >= 0 or [ b_i' A_i' ] [ x' ] >= 0
1293 * Then we compute the convex hull of the polytopes b_i' + A_i' x' >= 0
1294 * resulting in b' + A' x' >= 0, which we then convert back
1296 * [ y ] [ y ]
1297 * [ b' A' ] S [ x ] >= 0 or [ b A ] [ x ] >= 0
1299 * The polyhedron b + A x >= 0 is then the convex hull of the input polyhedra.
1301 static struct isl_basic_set *convex_hull_pair_pointed(
1302 struct isl_basic_set *bset1, struct isl_basic_set *bset2)
1304 struct isl_ctx *ctx = NULL;
1305 struct isl_vec *dir = NULL;
1306 struct isl_mat *T = NULL;
1307 struct isl_mat *T2 = NULL;
1308 struct isl_basic_set *hull;
1309 struct isl_set *set;
1311 if (!bset1 || !bset2)
1312 goto error;
1313 ctx = bset1->ctx;
1314 dir = valid_direction(isl_basic_set_copy(bset1),
1315 isl_basic_set_copy(bset2));
1316 if (!dir)
1317 goto error;
1318 T = isl_mat_alloc(bset1->ctx, dir->size, dir->size);
1319 if (!T)
1320 goto error;
1321 isl_seq_cpy(T->row[0], dir->block.data, dir->size);
1322 T = isl_mat_unimodular_complete(T, 1);
1323 T2 = isl_mat_right_inverse(isl_mat_copy(T));
1325 bset1 = homogeneous_map(bset1, isl_mat_copy(T2));
1326 bset2 = homogeneous_map(bset2, T2);
1327 set = isl_set_alloc_space(isl_basic_set_get_space(bset1), 2, 0);
1328 set = isl_set_add_basic_set(set, bset1);
1329 set = isl_set_add_basic_set(set, bset2);
1330 hull = uset_convex_hull(set);
1331 hull = isl_basic_set_preimage(hull, T);
1333 isl_vec_free(dir);
1335 return hull;
1336 error:
1337 isl_vec_free(dir);
1338 isl_basic_set_free(bset1);
1339 isl_basic_set_free(bset2);
1340 return NULL;
1343 static struct isl_basic_set *uset_convex_hull_wrap(struct isl_set *set);
1344 static struct isl_basic_set *modulo_affine_hull(
1345 struct isl_set *set, struct isl_basic_set *affine_hull);
1347 /* Compute the convex hull of a pair of basic sets without any parameters or
1348 * integer divisions.
1350 * This function is called from uset_convex_hull_unbounded, which
1351 * means that the complete convex hull is unbounded. Some pairs
1352 * of basic sets may still be bounded, though.
1353 * They may even lie inside a lower dimensional space, in which
1354 * case they need to be handled inside their affine hull since
1355 * the main algorithm assumes that the result is full-dimensional.
1357 * If the convex hull of the two basic sets would have a non-trivial
1358 * lineality space, we first project out this lineality space.
1360 static struct isl_basic_set *convex_hull_pair(struct isl_basic_set *bset1,
1361 struct isl_basic_set *bset2)
1363 isl_basic_set *lin, *aff;
1364 int bounded1, bounded2;
1366 if (bset1->ctx->opt->convex == ISL_CONVEX_HULL_FM)
1367 return convex_hull_pair_elim(bset1, bset2);
1369 aff = isl_set_affine_hull(isl_basic_set_union(isl_basic_set_copy(bset1),
1370 isl_basic_set_copy(bset2)));
1371 if (!aff)
1372 goto error;
1373 if (aff->n_eq != 0)
1374 return modulo_affine_hull(isl_basic_set_union(bset1, bset2), aff);
1375 isl_basic_set_free(aff);
1377 bounded1 = isl_basic_set_is_bounded(bset1);
1378 bounded2 = isl_basic_set_is_bounded(bset2);
1380 if (bounded1 < 0 || bounded2 < 0)
1381 goto error;
1383 if (bounded1 && bounded2)
1384 uset_convex_hull_wrap(isl_basic_set_union(bset1, bset2));
1386 if (bounded1 || bounded2)
1387 return convex_hull_pair_pointed(bset1, bset2);
1389 lin = induced_lineality_space(isl_basic_set_copy(bset1),
1390 isl_basic_set_copy(bset2));
1391 if (!lin)
1392 goto error;
1393 if (isl_basic_set_is_universe(lin)) {
1394 isl_basic_set_free(bset1);
1395 isl_basic_set_free(bset2);
1396 return lin;
1398 if (lin->n_eq < isl_basic_set_total_dim(lin)) {
1399 struct isl_set *set;
1400 set = isl_set_alloc_space(isl_basic_set_get_space(bset1), 2, 0);
1401 set = isl_set_add_basic_set(set, bset1);
1402 set = isl_set_add_basic_set(set, bset2);
1403 return modulo_lineality(set, lin);
1405 isl_basic_set_free(lin);
1407 return convex_hull_pair_pointed(bset1, bset2);
1408 error:
1409 isl_basic_set_free(bset1);
1410 isl_basic_set_free(bset2);
1411 return NULL;
1414 /* Compute the lineality space of a basic set.
1415 * We currently do not allow the basic set to have any divs.
1416 * We basically just drop the constants and turn every inequality
1417 * into an equality.
1419 struct isl_basic_set *isl_basic_set_lineality_space(struct isl_basic_set *bset)
1421 int i, k;
1422 struct isl_basic_set *lin = NULL;
1423 unsigned dim;
1425 if (!bset)
1426 goto error;
1427 isl_assert(bset->ctx, bset->n_div == 0, goto error);
1428 dim = isl_basic_set_total_dim(bset);
1430 lin = isl_basic_set_alloc_space(isl_basic_set_get_space(bset), 0, dim, 0);
1431 if (!lin)
1432 goto error;
1433 for (i = 0; i < bset->n_eq; ++i) {
1434 k = isl_basic_set_alloc_equality(lin);
1435 if (k < 0)
1436 goto error;
1437 isl_int_set_si(lin->eq[k][0], 0);
1438 isl_seq_cpy(lin->eq[k] + 1, bset->eq[i] + 1, dim);
1440 lin = isl_basic_set_gauss(lin, NULL);
1441 if (!lin)
1442 goto error;
1443 for (i = 0; i < bset->n_ineq && lin->n_eq < dim; ++i) {
1444 k = isl_basic_set_alloc_equality(lin);
1445 if (k < 0)
1446 goto error;
1447 isl_int_set_si(lin->eq[k][0], 0);
1448 isl_seq_cpy(lin->eq[k] + 1, bset->ineq[i] + 1, dim);
1449 lin = isl_basic_set_gauss(lin, NULL);
1450 if (!lin)
1451 goto error;
1453 isl_basic_set_free(bset);
1454 return lin;
1455 error:
1456 isl_basic_set_free(lin);
1457 isl_basic_set_free(bset);
1458 return NULL;
1461 /* Compute the (linear) hull of the lineality spaces of the basic sets in the
1462 * "underlying" set "set".
1464 static struct isl_basic_set *uset_combined_lineality_space(struct isl_set *set)
1466 int i;
1467 struct isl_set *lin = NULL;
1469 if (!set)
1470 return NULL;
1471 if (set->n == 0) {
1472 isl_space *dim = isl_set_get_space(set);
1473 isl_set_free(set);
1474 return isl_basic_set_empty(dim);
1477 lin = isl_set_alloc_space(isl_set_get_space(set), set->n, 0);
1478 for (i = 0; i < set->n; ++i)
1479 lin = isl_set_add_basic_set(lin,
1480 isl_basic_set_lineality_space(isl_basic_set_copy(set->p[i])));
1481 isl_set_free(set);
1482 return isl_set_affine_hull(lin);
1485 /* Compute the convex hull of a set without any parameters or
1486 * integer divisions.
1487 * In each step, we combined two basic sets until only one
1488 * basic set is left.
1489 * The input basic sets are assumed not to have a non-trivial
1490 * lineality space. If any of the intermediate results has
1491 * a non-trivial lineality space, it is projected out.
1493 static struct isl_basic_set *uset_convex_hull_unbounded(struct isl_set *set)
1495 struct isl_basic_set *convex_hull = NULL;
1497 convex_hull = isl_set_copy_basic_set(set);
1498 set = isl_set_drop_basic_set(set, convex_hull);
1499 if (!set)
1500 goto error;
1501 while (set->n > 0) {
1502 struct isl_basic_set *t;
1503 t = isl_set_copy_basic_set(set);
1504 if (!t)
1505 goto error;
1506 set = isl_set_drop_basic_set(set, t);
1507 if (!set)
1508 goto error;
1509 convex_hull = convex_hull_pair(convex_hull, t);
1510 if (set->n == 0)
1511 break;
1512 t = isl_basic_set_lineality_space(isl_basic_set_copy(convex_hull));
1513 if (!t)
1514 goto error;
1515 if (isl_basic_set_is_universe(t)) {
1516 isl_basic_set_free(convex_hull);
1517 convex_hull = t;
1518 break;
1520 if (t->n_eq < isl_basic_set_total_dim(t)) {
1521 set = isl_set_add_basic_set(set, convex_hull);
1522 return modulo_lineality(set, t);
1524 isl_basic_set_free(t);
1526 isl_set_free(set);
1527 return convex_hull;
1528 error:
1529 isl_set_free(set);
1530 isl_basic_set_free(convex_hull);
1531 return NULL;
1534 /* Compute an initial hull for wrapping containing a single initial
1535 * facet.
1536 * This function assumes that the given set is bounded.
1538 static struct isl_basic_set *initial_hull(struct isl_basic_set *hull,
1539 struct isl_set *set)
1541 struct isl_mat *bounds = NULL;
1542 unsigned dim;
1543 int k;
1545 if (!hull)
1546 goto error;
1547 bounds = initial_facet_constraint(set);
1548 if (!bounds)
1549 goto error;
1550 k = isl_basic_set_alloc_inequality(hull);
1551 if (k < 0)
1552 goto error;
1553 dim = isl_set_n_dim(set);
1554 isl_assert(set->ctx, 1 + dim == bounds->n_col, goto error);
1555 isl_seq_cpy(hull->ineq[k], bounds->row[0], bounds->n_col);
1556 isl_mat_free(bounds);
1558 return hull;
1559 error:
1560 isl_basic_set_free(hull);
1561 isl_mat_free(bounds);
1562 return NULL;
1565 struct max_constraint {
1566 struct isl_mat *c;
1567 int count;
1568 int ineq;
1571 static int max_constraint_equal(const void *entry, const void *val)
1573 struct max_constraint *a = (struct max_constraint *)entry;
1574 isl_int *b = (isl_int *)val;
1576 return isl_seq_eq(a->c->row[0] + 1, b, a->c->n_col - 1);
1579 static void update_constraint(struct isl_ctx *ctx, struct isl_hash_table *table,
1580 isl_int *con, unsigned len, int n, int ineq)
1582 struct isl_hash_table_entry *entry;
1583 struct max_constraint *c;
1584 uint32_t c_hash;
1586 c_hash = isl_seq_get_hash(con + 1, len);
1587 entry = isl_hash_table_find(ctx, table, c_hash, max_constraint_equal,
1588 con + 1, 0);
1589 if (!entry)
1590 return;
1591 c = entry->data;
1592 if (c->count < n) {
1593 isl_hash_table_remove(ctx, table, entry);
1594 return;
1596 c->count++;
1597 if (isl_int_gt(c->c->row[0][0], con[0]))
1598 return;
1599 if (isl_int_eq(c->c->row[0][0], con[0])) {
1600 if (ineq)
1601 c->ineq = ineq;
1602 return;
1604 c->c = isl_mat_cow(c->c);
1605 isl_int_set(c->c->row[0][0], con[0]);
1606 c->ineq = ineq;
1609 /* Check whether the constraint hash table "table" constains the constraint
1610 * "con".
1612 static int has_constraint(struct isl_ctx *ctx, struct isl_hash_table *table,
1613 isl_int *con, unsigned len, int n)
1615 struct isl_hash_table_entry *entry;
1616 struct max_constraint *c;
1617 uint32_t c_hash;
1619 c_hash = isl_seq_get_hash(con + 1, len);
1620 entry = isl_hash_table_find(ctx, table, c_hash, max_constraint_equal,
1621 con + 1, 0);
1622 if (!entry)
1623 return 0;
1624 c = entry->data;
1625 if (c->count < n)
1626 return 0;
1627 return isl_int_eq(c->c->row[0][0], con[0]);
1630 /* Check for inequality constraints of a basic set without equalities
1631 * such that the same or more stringent copies of the constraint appear
1632 * in all of the basic sets. Such constraints are necessarily facet
1633 * constraints of the convex hull.
1635 * If the resulting basic set is by chance identical to one of
1636 * the basic sets in "set", then we know that this basic set contains
1637 * all other basic sets and is therefore the convex hull of set.
1638 * In this case we set *is_hull to 1.
1640 static struct isl_basic_set *common_constraints(struct isl_basic_set *hull,
1641 struct isl_set *set, int *is_hull)
1643 int i, j, s, n;
1644 int min_constraints;
1645 int best;
1646 struct max_constraint *constraints = NULL;
1647 struct isl_hash_table *table = NULL;
1648 unsigned total;
1650 *is_hull = 0;
1652 for (i = 0; i < set->n; ++i)
1653 if (set->p[i]->n_eq == 0)
1654 break;
1655 if (i >= set->n)
1656 return hull;
1657 min_constraints = set->p[i]->n_ineq;
1658 best = i;
1659 for (i = best + 1; i < set->n; ++i) {
1660 if (set->p[i]->n_eq != 0)
1661 continue;
1662 if (set->p[i]->n_ineq >= min_constraints)
1663 continue;
1664 min_constraints = set->p[i]->n_ineq;
1665 best = i;
1667 constraints = isl_calloc_array(hull->ctx, struct max_constraint,
1668 min_constraints);
1669 if (!constraints)
1670 return hull;
1671 table = isl_alloc_type(hull->ctx, struct isl_hash_table);
1672 if (isl_hash_table_init(hull->ctx, table, min_constraints))
1673 goto error;
1675 total = isl_space_dim(set->dim, isl_dim_all);
1676 for (i = 0; i < set->p[best]->n_ineq; ++i) {
1677 constraints[i].c = isl_mat_sub_alloc6(hull->ctx,
1678 set->p[best]->ineq + i, 0, 1, 0, 1 + total);
1679 if (!constraints[i].c)
1680 goto error;
1681 constraints[i].ineq = 1;
1683 for (i = 0; i < min_constraints; ++i) {
1684 struct isl_hash_table_entry *entry;
1685 uint32_t c_hash;
1686 c_hash = isl_seq_get_hash(constraints[i].c->row[0] + 1, total);
1687 entry = isl_hash_table_find(hull->ctx, table, c_hash,
1688 max_constraint_equal, constraints[i].c->row[0] + 1, 1);
1689 if (!entry)
1690 goto error;
1691 isl_assert(hull->ctx, !entry->data, goto error);
1692 entry->data = &constraints[i];
1695 n = 0;
1696 for (s = 0; s < set->n; ++s) {
1697 if (s == best)
1698 continue;
1700 for (i = 0; i < set->p[s]->n_eq; ++i) {
1701 isl_int *eq = set->p[s]->eq[i];
1702 for (j = 0; j < 2; ++j) {
1703 isl_seq_neg(eq, eq, 1 + total);
1704 update_constraint(hull->ctx, table,
1705 eq, total, n, 0);
1708 for (i = 0; i < set->p[s]->n_ineq; ++i) {
1709 isl_int *ineq = set->p[s]->ineq[i];
1710 update_constraint(hull->ctx, table, ineq, total, n,
1711 set->p[s]->n_eq == 0);
1713 ++n;
1716 for (i = 0; i < min_constraints; ++i) {
1717 if (constraints[i].count < n)
1718 continue;
1719 if (!constraints[i].ineq)
1720 continue;
1721 j = isl_basic_set_alloc_inequality(hull);
1722 if (j < 0)
1723 goto error;
1724 isl_seq_cpy(hull->ineq[j], constraints[i].c->row[0], 1 + total);
1727 for (s = 0; s < set->n; ++s) {
1728 if (set->p[s]->n_eq)
1729 continue;
1730 if (set->p[s]->n_ineq != hull->n_ineq)
1731 continue;
1732 for (i = 0; i < set->p[s]->n_ineq; ++i) {
1733 isl_int *ineq = set->p[s]->ineq[i];
1734 if (!has_constraint(hull->ctx, table, ineq, total, n))
1735 break;
1737 if (i == set->p[s]->n_ineq)
1738 *is_hull = 1;
1741 isl_hash_table_clear(table);
1742 for (i = 0; i < min_constraints; ++i)
1743 isl_mat_free(constraints[i].c);
1744 free(constraints);
1745 free(table);
1746 return hull;
1747 error:
1748 isl_hash_table_clear(table);
1749 free(table);
1750 if (constraints)
1751 for (i = 0; i < min_constraints; ++i)
1752 isl_mat_free(constraints[i].c);
1753 free(constraints);
1754 return hull;
1757 /* Create a template for the convex hull of "set" and fill it up
1758 * obvious facet constraints, if any. If the result happens to
1759 * be the convex hull of "set" then *is_hull is set to 1.
1761 static struct isl_basic_set *proto_hull(struct isl_set *set, int *is_hull)
1763 struct isl_basic_set *hull;
1764 unsigned n_ineq;
1765 int i;
1767 n_ineq = 1;
1768 for (i = 0; i < set->n; ++i) {
1769 n_ineq += set->p[i]->n_eq;
1770 n_ineq += set->p[i]->n_ineq;
1772 hull = isl_basic_set_alloc_space(isl_space_copy(set->dim), 0, 0, n_ineq);
1773 hull = isl_basic_set_set_rational(hull);
1774 if (!hull)
1775 return NULL;
1776 return common_constraints(hull, set, is_hull);
1779 static struct isl_basic_set *uset_convex_hull_wrap(struct isl_set *set)
1781 struct isl_basic_set *hull;
1782 int is_hull;
1784 hull = proto_hull(set, &is_hull);
1785 if (hull && !is_hull) {
1786 if (hull->n_ineq == 0)
1787 hull = initial_hull(hull, set);
1788 hull = extend(hull, set);
1790 isl_set_free(set);
1792 return hull;
1795 /* Compute the convex hull of a set without any parameters or
1796 * integer divisions. Depending on whether the set is bounded,
1797 * we pass control to the wrapping based convex hull or
1798 * the Fourier-Motzkin elimination based convex hull.
1799 * We also handle a few special cases before checking the boundedness.
1801 static struct isl_basic_set *uset_convex_hull(struct isl_set *set)
1803 struct isl_basic_set *convex_hull = NULL;
1804 struct isl_basic_set *lin;
1806 if (isl_set_n_dim(set) == 0)
1807 return convex_hull_0d(set);
1809 set = isl_set_coalesce(set);
1810 set = isl_set_set_rational(set);
1812 if (!set)
1813 goto error;
1814 if (!set)
1815 return NULL;
1816 if (set->n == 1) {
1817 convex_hull = isl_basic_set_copy(set->p[0]);
1818 isl_set_free(set);
1819 return convex_hull;
1821 if (isl_set_n_dim(set) == 1)
1822 return convex_hull_1d(set);
1824 if (isl_set_is_bounded(set) &&
1825 set->ctx->opt->convex == ISL_CONVEX_HULL_WRAP)
1826 return uset_convex_hull_wrap(set);
1828 lin = uset_combined_lineality_space(isl_set_copy(set));
1829 if (!lin)
1830 goto error;
1831 if (isl_basic_set_is_universe(lin)) {
1832 isl_set_free(set);
1833 return lin;
1835 if (lin->n_eq < isl_basic_set_total_dim(lin))
1836 return modulo_lineality(set, lin);
1837 isl_basic_set_free(lin);
1839 return uset_convex_hull_unbounded(set);
1840 error:
1841 isl_set_free(set);
1842 isl_basic_set_free(convex_hull);
1843 return NULL;
1846 /* This is the core procedure, where "set" is a "pure" set, i.e.,
1847 * without parameters or divs and where the convex hull of set is
1848 * known to be full-dimensional.
1850 static struct isl_basic_set *uset_convex_hull_wrap_bounded(struct isl_set *set)
1852 struct isl_basic_set *convex_hull = NULL;
1854 if (!set)
1855 goto error;
1857 if (isl_set_n_dim(set) == 0) {
1858 convex_hull = isl_basic_set_universe(isl_space_copy(set->dim));
1859 isl_set_free(set);
1860 convex_hull = isl_basic_set_set_rational(convex_hull);
1861 return convex_hull;
1864 set = isl_set_set_rational(set);
1865 set = isl_set_coalesce(set);
1866 if (!set)
1867 goto error;
1868 if (set->n == 1) {
1869 convex_hull = isl_basic_set_copy(set->p[0]);
1870 isl_set_free(set);
1871 convex_hull = isl_basic_map_remove_redundancies(convex_hull);
1872 return convex_hull;
1874 if (isl_set_n_dim(set) == 1)
1875 return convex_hull_1d(set);
1877 return uset_convex_hull_wrap(set);
1878 error:
1879 isl_set_free(set);
1880 return NULL;
1883 /* Compute the convex hull of set "set" with affine hull "affine_hull",
1884 * We first remove the equalities (transforming the set), compute the
1885 * convex hull of the transformed set and then add the equalities back
1886 * (after performing the inverse transformation.
1888 static struct isl_basic_set *modulo_affine_hull(
1889 struct isl_set *set, struct isl_basic_set *affine_hull)
1891 struct isl_mat *T;
1892 struct isl_mat *T2;
1893 struct isl_basic_set *dummy;
1894 struct isl_basic_set *convex_hull;
1896 dummy = isl_basic_set_remove_equalities(
1897 isl_basic_set_copy(affine_hull), &T, &T2);
1898 if (!dummy)
1899 goto error;
1900 isl_basic_set_free(dummy);
1901 set = isl_set_preimage(set, T);
1902 convex_hull = uset_convex_hull(set);
1903 convex_hull = isl_basic_set_preimage(convex_hull, T2);
1904 convex_hull = isl_basic_set_intersect(convex_hull, affine_hull);
1905 return convex_hull;
1906 error:
1907 isl_basic_set_free(affine_hull);
1908 isl_set_free(set);
1909 return NULL;
1912 /* Compute the convex hull of a map.
1914 * The implementation was inspired by "Extended Convex Hull" by Fukuda et al.,
1915 * specifically, the wrapping of facets to obtain new facets.
1917 struct isl_basic_map *isl_map_convex_hull(struct isl_map *map)
1919 struct isl_basic_set *bset;
1920 struct isl_basic_map *model = NULL;
1921 struct isl_basic_set *affine_hull = NULL;
1922 struct isl_basic_map *convex_hull = NULL;
1923 struct isl_set *set = NULL;
1924 struct isl_ctx *ctx;
1926 if (!map)
1927 goto error;
1929 ctx = map->ctx;
1930 if (map->n == 0) {
1931 convex_hull = isl_basic_map_empty_like_map(map);
1932 isl_map_free(map);
1933 return convex_hull;
1936 map = isl_map_detect_equalities(map);
1937 map = isl_map_align_divs(map);
1938 if (!map)
1939 goto error;
1940 model = isl_basic_map_copy(map->p[0]);
1941 set = isl_map_underlying_set(map);
1942 if (!set)
1943 goto error;
1945 affine_hull = isl_set_affine_hull(isl_set_copy(set));
1946 if (!affine_hull)
1947 goto error;
1948 if (affine_hull->n_eq != 0)
1949 bset = modulo_affine_hull(set, affine_hull);
1950 else {
1951 isl_basic_set_free(affine_hull);
1952 bset = uset_convex_hull(set);
1955 convex_hull = isl_basic_map_overlying_set(bset, model);
1956 if (!convex_hull)
1957 return NULL;
1959 ISL_F_SET(convex_hull, ISL_BASIC_MAP_NO_IMPLICIT);
1960 ISL_F_SET(convex_hull, ISL_BASIC_MAP_ALL_EQUALITIES);
1961 ISL_F_CLR(convex_hull, ISL_BASIC_MAP_RATIONAL);
1962 return convex_hull;
1963 error:
1964 isl_set_free(set);
1965 isl_basic_map_free(model);
1966 return NULL;
1969 struct isl_basic_set *isl_set_convex_hull(struct isl_set *set)
1971 return (struct isl_basic_set *)
1972 isl_map_convex_hull((struct isl_map *)set);
1975 __isl_give isl_basic_map *isl_map_polyhedral_hull(__isl_take isl_map *map)
1977 isl_basic_map *hull;
1979 hull = isl_map_convex_hull(map);
1980 return isl_basic_map_remove_divs(hull);
1983 __isl_give isl_basic_set *isl_set_polyhedral_hull(__isl_take isl_set *set)
1985 return (isl_basic_set *)isl_map_polyhedral_hull((isl_map *)set);
1988 struct sh_data_entry {
1989 struct isl_hash_table *table;
1990 struct isl_tab *tab;
1993 /* Holds the data needed during the simple hull computation.
1994 * In particular,
1995 * n the number of basic sets in the original set
1996 * hull_table a hash table of already computed constraints
1997 * in the simple hull
1998 * p for each basic set,
1999 * table a hash table of the constraints
2000 * tab the tableau corresponding to the basic set
2002 struct sh_data {
2003 struct isl_ctx *ctx;
2004 unsigned n;
2005 struct isl_hash_table *hull_table;
2006 struct sh_data_entry p[1];
2009 static void sh_data_free(struct sh_data *data)
2011 int i;
2013 if (!data)
2014 return;
2015 isl_hash_table_free(data->ctx, data->hull_table);
2016 for (i = 0; i < data->n; ++i) {
2017 isl_hash_table_free(data->ctx, data->p[i].table);
2018 isl_tab_free(data->p[i].tab);
2020 free(data);
2023 struct ineq_cmp_data {
2024 unsigned len;
2025 isl_int *p;
2028 static int has_ineq(const void *entry, const void *val)
2030 isl_int *row = (isl_int *)entry;
2031 struct ineq_cmp_data *v = (struct ineq_cmp_data *)val;
2033 return isl_seq_eq(row + 1, v->p + 1, v->len) ||
2034 isl_seq_is_neg(row + 1, v->p + 1, v->len);
2037 static int hash_ineq(struct isl_ctx *ctx, struct isl_hash_table *table,
2038 isl_int *ineq, unsigned len)
2040 uint32_t c_hash;
2041 struct ineq_cmp_data v;
2042 struct isl_hash_table_entry *entry;
2044 v.len = len;
2045 v.p = ineq;
2046 c_hash = isl_seq_get_hash(ineq + 1, len);
2047 entry = isl_hash_table_find(ctx, table, c_hash, has_ineq, &v, 1);
2048 if (!entry)
2049 return - 1;
2050 entry->data = ineq;
2051 return 0;
2054 /* Fill hash table "table" with the constraints of "bset".
2055 * Equalities are added as two inequalities.
2056 * The value in the hash table is a pointer to the (in)equality of "bset".
2058 static int hash_basic_set(struct isl_hash_table *table,
2059 struct isl_basic_set *bset)
2061 int i, j;
2062 unsigned dim = isl_basic_set_total_dim(bset);
2064 for (i = 0; i < bset->n_eq; ++i) {
2065 for (j = 0; j < 2; ++j) {
2066 isl_seq_neg(bset->eq[i], bset->eq[i], 1 + dim);
2067 if (hash_ineq(bset->ctx, table, bset->eq[i], dim) < 0)
2068 return -1;
2071 for (i = 0; i < bset->n_ineq; ++i) {
2072 if (hash_ineq(bset->ctx, table, bset->ineq[i], dim) < 0)
2073 return -1;
2075 return 0;
2078 static struct sh_data *sh_data_alloc(struct isl_set *set, unsigned n_ineq)
2080 struct sh_data *data;
2081 int i;
2083 data = isl_calloc(set->ctx, struct sh_data,
2084 sizeof(struct sh_data) +
2085 (set->n - 1) * sizeof(struct sh_data_entry));
2086 if (!data)
2087 return NULL;
2088 data->ctx = set->ctx;
2089 data->n = set->n;
2090 data->hull_table = isl_hash_table_alloc(set->ctx, n_ineq);
2091 if (!data->hull_table)
2092 goto error;
2093 for (i = 0; i < set->n; ++i) {
2094 data->p[i].table = isl_hash_table_alloc(set->ctx,
2095 2 * set->p[i]->n_eq + set->p[i]->n_ineq);
2096 if (!data->p[i].table)
2097 goto error;
2098 if (hash_basic_set(data->p[i].table, set->p[i]) < 0)
2099 goto error;
2101 return data;
2102 error:
2103 sh_data_free(data);
2104 return NULL;
2107 /* Check if inequality "ineq" is a bound for basic set "j" or if
2108 * it can be relaxed (by increasing the constant term) to become
2109 * a bound for that basic set. In the latter case, the constant
2110 * term is updated.
2111 * Relaxation of the constant term is only allowed if "shift" is set.
2113 * Return 1 if "ineq" is a bound
2114 * 0 if "ineq" may attain arbitrarily small values on basic set "j"
2115 * -1 if some error occurred
2117 static int is_bound(struct sh_data *data, struct isl_set *set, int j,
2118 isl_int *ineq, int shift)
2120 enum isl_lp_result res;
2121 isl_int opt;
2123 if (!data->p[j].tab) {
2124 data->p[j].tab = isl_tab_from_basic_set(set->p[j], 0);
2125 if (!data->p[j].tab)
2126 return -1;
2129 isl_int_init(opt);
2131 res = isl_tab_min(data->p[j].tab, ineq, data->ctx->one,
2132 &opt, NULL, 0);
2133 if (res == isl_lp_ok && isl_int_is_neg(opt)) {
2134 if (shift)
2135 isl_int_sub(ineq[0], ineq[0], opt);
2136 else
2137 res = isl_lp_unbounded;
2140 isl_int_clear(opt);
2142 return (res == isl_lp_ok || res == isl_lp_empty) ? 1 :
2143 res == isl_lp_unbounded ? 0 : -1;
2146 /* Check if inequality "ineq" from basic set "i" is or can be relaxed to
2147 * become a bound on the whole set. If so, add the (relaxed) inequality
2148 * to "hull". Relaxation is only allowed if "shift" is set.
2150 * We first check if "hull" already contains a translate of the inequality.
2151 * If so, we are done.
2152 * Then, we check if any of the previous basic sets contains a translate
2153 * of the inequality. If so, then we have already considered this
2154 * inequality and we are done.
2155 * Otherwise, for each basic set other than "i", we check if the inequality
2156 * is a bound on the basic set.
2157 * For previous basic sets, we know that they do not contain a translate
2158 * of the inequality, so we directly call is_bound.
2159 * For following basic sets, we first check if a translate of the
2160 * inequality appears in its description and if so directly update
2161 * the inequality accordingly.
2163 static struct isl_basic_set *add_bound(struct isl_basic_set *hull,
2164 struct sh_data *data, struct isl_set *set, int i, isl_int *ineq,
2165 int shift)
2167 uint32_t c_hash;
2168 struct ineq_cmp_data v;
2169 struct isl_hash_table_entry *entry;
2170 int j, k;
2172 if (!hull)
2173 return NULL;
2175 v.len = isl_basic_set_total_dim(hull);
2176 v.p = ineq;
2177 c_hash = isl_seq_get_hash(ineq + 1, v.len);
2179 entry = isl_hash_table_find(hull->ctx, data->hull_table, c_hash,
2180 has_ineq, &v, 0);
2181 if (entry)
2182 return hull;
2184 for (j = 0; j < i; ++j) {
2185 entry = isl_hash_table_find(hull->ctx, data->p[j].table,
2186 c_hash, has_ineq, &v, 0);
2187 if (entry)
2188 break;
2190 if (j < i)
2191 return hull;
2193 k = isl_basic_set_alloc_inequality(hull);
2194 isl_seq_cpy(hull->ineq[k], ineq, 1 + v.len);
2195 if (k < 0)
2196 goto error;
2198 for (j = 0; j < i; ++j) {
2199 int bound;
2200 bound = is_bound(data, set, j, hull->ineq[k], shift);
2201 if (bound < 0)
2202 goto error;
2203 if (!bound)
2204 break;
2206 if (j < i) {
2207 isl_basic_set_free_inequality(hull, 1);
2208 return hull;
2211 for (j = i + 1; j < set->n; ++j) {
2212 int bound, neg;
2213 isl_int *ineq_j;
2214 entry = isl_hash_table_find(hull->ctx, data->p[j].table,
2215 c_hash, has_ineq, &v, 0);
2216 if (entry) {
2217 ineq_j = entry->data;
2218 neg = isl_seq_is_neg(ineq_j + 1,
2219 hull->ineq[k] + 1, v.len);
2220 if (neg)
2221 isl_int_neg(ineq_j[0], ineq_j[0]);
2222 if (isl_int_gt(ineq_j[0], hull->ineq[k][0]))
2223 isl_int_set(hull->ineq[k][0], ineq_j[0]);
2224 if (neg)
2225 isl_int_neg(ineq_j[0], ineq_j[0]);
2226 continue;
2228 bound = is_bound(data, set, j, hull->ineq[k], shift);
2229 if (bound < 0)
2230 goto error;
2231 if (!bound)
2232 break;
2234 if (j < set->n) {
2235 isl_basic_set_free_inequality(hull, 1);
2236 return hull;
2239 entry = isl_hash_table_find(hull->ctx, data->hull_table, c_hash,
2240 has_ineq, &v, 1);
2241 if (!entry)
2242 goto error;
2243 entry->data = hull->ineq[k];
2245 return hull;
2246 error:
2247 isl_basic_set_free(hull);
2248 return NULL;
2251 /* Check if any inequality from basic set "i" is or can be relaxed to
2252 * become a bound on the whole set. If so, add the (relaxed) inequality
2253 * to "hull". Relaxation is only allowed if "shift" is set.
2255 static struct isl_basic_set *add_bounds(struct isl_basic_set *bset,
2256 struct sh_data *data, struct isl_set *set, int i, int shift)
2258 int j, k;
2259 unsigned dim = isl_basic_set_total_dim(bset);
2261 for (j = 0; j < set->p[i]->n_eq; ++j) {
2262 for (k = 0; k < 2; ++k) {
2263 isl_seq_neg(set->p[i]->eq[j], set->p[i]->eq[j], 1+dim);
2264 bset = add_bound(bset, data, set, i, set->p[i]->eq[j],
2265 shift);
2268 for (j = 0; j < set->p[i]->n_ineq; ++j)
2269 bset = add_bound(bset, data, set, i, set->p[i]->ineq[j], shift);
2270 return bset;
2273 /* Compute a superset of the convex hull of set that is described
2274 * by only (translates of) the constraints in the constituents of set.
2275 * Translation is only allowed if "shift" is set.
2277 static __isl_give isl_basic_set *uset_simple_hull(__isl_take isl_set *set,
2278 int shift)
2280 struct sh_data *data = NULL;
2281 struct isl_basic_set *hull = NULL;
2282 unsigned n_ineq;
2283 int i;
2285 if (!set)
2286 return NULL;
2288 n_ineq = 0;
2289 for (i = 0; i < set->n; ++i) {
2290 if (!set->p[i])
2291 goto error;
2292 n_ineq += 2 * set->p[i]->n_eq + set->p[i]->n_ineq;
2295 hull = isl_basic_set_alloc_space(isl_space_copy(set->dim), 0, 0, n_ineq);
2296 if (!hull)
2297 goto error;
2299 data = sh_data_alloc(set, n_ineq);
2300 if (!data)
2301 goto error;
2303 for (i = 0; i < set->n; ++i)
2304 hull = add_bounds(hull, data, set, i, shift);
2306 sh_data_free(data);
2307 isl_set_free(set);
2309 return hull;
2310 error:
2311 sh_data_free(data);
2312 isl_basic_set_free(hull);
2313 isl_set_free(set);
2314 return NULL;
2317 /* Compute a superset of the convex hull of map that is described
2318 * by only (translates of) the constraints in the constituents of map.
2319 * Translation is only allowed if "shift" is set.
2321 static __isl_give isl_basic_map *map_simple_hull(__isl_take isl_map *map,
2322 int shift)
2324 struct isl_set *set = NULL;
2325 struct isl_basic_map *model = NULL;
2326 struct isl_basic_map *hull;
2327 struct isl_basic_map *affine_hull;
2328 struct isl_basic_set *bset = NULL;
2330 if (!map)
2331 return NULL;
2332 if (map->n == 0) {
2333 hull = isl_basic_map_empty_like_map(map);
2334 isl_map_free(map);
2335 return hull;
2337 if (map->n == 1) {
2338 hull = isl_basic_map_copy(map->p[0]);
2339 isl_map_free(map);
2340 return hull;
2343 map = isl_map_detect_equalities(map);
2344 affine_hull = isl_map_affine_hull(isl_map_copy(map));
2345 map = isl_map_align_divs(map);
2346 model = map ? isl_basic_map_copy(map->p[0]) : NULL;
2348 set = isl_map_underlying_set(map);
2350 bset = uset_simple_hull(set, shift);
2352 hull = isl_basic_map_overlying_set(bset, model);
2354 hull = isl_basic_map_intersect(hull, affine_hull);
2355 hull = isl_basic_map_remove_redundancies(hull);
2357 if (!hull)
2358 return NULL;
2359 ISL_F_SET(hull, ISL_BASIC_MAP_NO_IMPLICIT);
2360 ISL_F_SET(hull, ISL_BASIC_MAP_ALL_EQUALITIES);
2362 return hull;
2365 /* Compute a superset of the convex hull of map that is described
2366 * by only translates of the constraints in the constituents of map.
2368 __isl_give isl_basic_map *isl_map_simple_hull(__isl_take isl_map *map)
2370 return map_simple_hull(map, 1);
2373 struct isl_basic_set *isl_set_simple_hull(struct isl_set *set)
2375 return (struct isl_basic_set *)
2376 isl_map_simple_hull((struct isl_map *)set);
2379 /* Compute a superset of the convex hull of map that is described
2380 * by only the constraints in the constituents of map.
2382 __isl_give isl_basic_map *isl_map_unshifted_simple_hull(
2383 __isl_take isl_map *map)
2385 return map_simple_hull(map, 0);
2388 __isl_give isl_basic_set *isl_set_unshifted_simple_hull(
2389 __isl_take isl_set *set)
2391 return isl_map_unshifted_simple_hull(set);
2394 /* Given a set "set", return parametric bounds on the dimension "dim".
2396 static struct isl_basic_set *set_bounds(struct isl_set *set, int dim)
2398 unsigned set_dim = isl_set_dim(set, isl_dim_set);
2399 set = isl_set_copy(set);
2400 set = isl_set_eliminate_dims(set, dim + 1, set_dim - (dim + 1));
2401 set = isl_set_eliminate_dims(set, 0, dim);
2402 return isl_set_convex_hull(set);
2405 /* Computes a "simple hull" and then check if each dimension in the
2406 * resulting hull is bounded by a symbolic constant. If not, the
2407 * hull is intersected with the corresponding bounds on the whole set.
2409 struct isl_basic_set *isl_set_bounded_simple_hull(struct isl_set *set)
2411 int i, j;
2412 struct isl_basic_set *hull;
2413 unsigned nparam, left;
2414 int removed_divs = 0;
2416 hull = isl_set_simple_hull(isl_set_copy(set));
2417 if (!hull)
2418 goto error;
2420 nparam = isl_basic_set_dim(hull, isl_dim_param);
2421 for (i = 0; i < isl_basic_set_dim(hull, isl_dim_set); ++i) {
2422 int lower = 0, upper = 0;
2423 struct isl_basic_set *bounds;
2425 left = isl_basic_set_total_dim(hull) - nparam - i - 1;
2426 for (j = 0; j < hull->n_eq; ++j) {
2427 if (isl_int_is_zero(hull->eq[j][1 + nparam + i]))
2428 continue;
2429 if (isl_seq_first_non_zero(hull->eq[j]+1+nparam+i+1,
2430 left) == -1)
2431 break;
2433 if (j < hull->n_eq)
2434 continue;
2436 for (j = 0; j < hull->n_ineq; ++j) {
2437 if (isl_int_is_zero(hull->ineq[j][1 + nparam + i]))
2438 continue;
2439 if (isl_seq_first_non_zero(hull->ineq[j]+1+nparam+i+1,
2440 left) != -1 ||
2441 isl_seq_first_non_zero(hull->ineq[j]+1+nparam,
2442 i) != -1)
2443 continue;
2444 if (isl_int_is_pos(hull->ineq[j][1 + nparam + i]))
2445 lower = 1;
2446 else
2447 upper = 1;
2448 if (lower && upper)
2449 break;
2452 if (lower && upper)
2453 continue;
2455 if (!removed_divs) {
2456 set = isl_set_remove_divs(set);
2457 if (!set)
2458 goto error;
2459 removed_divs = 1;
2461 bounds = set_bounds(set, i);
2462 hull = isl_basic_set_intersect(hull, bounds);
2463 if (!hull)
2464 goto error;
2467 isl_set_free(set);
2468 return hull;
2469 error:
2470 isl_set_free(set);
2471 return NULL;