Coverage Report

Created: 2017-04-29 12:21

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