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

source code of polly/lib/External/isl/isl_convex_hull.c