Coverage Report

Created: 2017-03-27 23:01

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