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