Coverage Report

Created: 2017-03-27 23:01

/Users/buildslave/jenkins/sharedspace/clang-stage2-coverage-R@2/llvm/tools/polly/lib/External/isl/isl_map_simplify.c
Line
Count
Source (jump to first uncovered line)
1
/*
2
 * Copyright 2008-2009 Katholieke Universiteit Leuven
3
 * Copyright 2012-2013 Ecole Normale Superieure
4
 * Copyright 2014-2015 INRIA Rocquencourt
5
 * Copyright 2016      Sven Verdoolaege
6
 *
7
 * Use of this software is governed by the MIT license
8
 *
9
 * Written by Sven Verdoolaege, K.U.Leuven, Departement
10
 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
11
 * and Ecole Normale Superieure, 45 rue d’Ulm, 75230 Paris, France
12
 * and Inria Paris - Rocquencourt, Domaine de Voluceau - Rocquencourt,
13
 * B.P. 105 - 78153 Le Chesnay, France
14
 */
15
16
#include <isl_ctx_private.h>
17
#include <isl_map_private.h>
18
#include "isl_equalities.h"
19
#include <isl/map.h>
20
#include <isl_seq.h>
21
#include "isl_tab.h"
22
#include <isl_space_private.h>
23
#include <isl_mat_private.h>
24
#include <isl_vec_private.h>
25
26
#include <bset_to_bmap.c>
27
#include <bset_from_bmap.c>
28
#include <set_to_map.c>
29
#include <set_from_map.c>
30
31
static void swap_equality(struct isl_basic_map *bmap, int a, int b)
32
362k
{
33
362k
  isl_int *t = bmap->eq[a];
34
362k
  bmap->eq[a] = bmap->eq[b];
35
362k
  bmap->eq[b] = t;
36
362k
}
37
38
static void swap_inequality(struct isl_basic_map *bmap, int a, int b)
39
37.8k
{
40
37.8k
  if (
a != b37.8k
)
{37.8k
41
37.8k
    isl_int *t = bmap->ineq[a];
42
37.8k
    bmap->ineq[a] = bmap->ineq[b];
43
37.8k
    bmap->ineq[b] = t;
44
37.8k
  }
45
37.8k
}
46
47
static void constraint_drop_vars(isl_int *c, unsigned n, unsigned rem)
48
2.79M
{
49
2.79M
  isl_seq_cpy(c, c + n, rem);
50
2.79M
  isl_seq_clr(c + rem, n);
51
2.79M
}
52
53
/* Drop n dimensions starting at first.
54
 *
55
 * In principle, this frees up some extra variables as the number
56
 * of columns remains constant, but we would have to extend
57
 * the div array too as the number of rows in this array is assumed
58
 * to be equal to extra.
59
 */
60
struct isl_basic_set *isl_basic_set_drop_dims(
61
    struct isl_basic_set *bset, unsigned first, unsigned n)
62
104k
{
63
104k
  return isl_basic_map_drop(bset_to_bmap(bset), isl_dim_set, first, n);
64
104k
}
65
66
/* Move "n" divs starting at "first" to the end of the list of divs.
67
 */
68
static struct isl_basic_map *move_divs_last(struct isl_basic_map *bmap,
69
  unsigned first, unsigned n)
70
448
{
71
448
  isl_int **div;
72
448
  int i;
73
448
74
448
  if (first + n == bmap->n_div)
75
376
    return bmap;
76
448
77
72
  
div = 72
isl_alloc_array72
(bmap->ctx, isl_int *, n);
78
72
  if (!div)
79
0
    goto error;
80
144
  
for (i = 0; 72
i < n144
;
++i72
)
81
72
    div[i] = bmap->div[first + i];
82
188
  for (i = 0; 
i < bmap->n_div - first - n188
;
++i116
)
83
116
    bmap->div[first + i] = bmap->div[first + n + i];
84
144
  for (i = 0; 
i < n144
;
++i72
)
85
72
    bmap->div[bmap->n_div - n + i] = div[i];
86
72
  free(div);
87
72
  return bmap;
88
0
error:
89
0
  isl_basic_map_free(bmap);
90
0
  return NULL;
91
72
}
92
93
/* Drop "n" dimensions of type "type" starting at "first".
94
 *
95
 * In principle, this frees up some extra variables as the number
96
 * of columns remains constant, but we would have to extend
97
 * the div array too as the number of rows in this array is assumed
98
 * to be equal to extra.
99
 */
100
struct isl_basic_map *isl_basic_map_drop(struct isl_basic_map *bmap,
101
  enum isl_dim_type type, unsigned first, unsigned n)
102
273k
{
103
273k
  int i;
104
273k
  unsigned dim;
105
273k
  unsigned offset;
106
273k
  unsigned left;
107
273k
108
273k
  if (!bmap)
109
0
    goto error;
110
273k
111
273k
  dim = isl_basic_map_dim(bmap, type);
112
273k
  isl_assert(bmap->ctx, first + n <= dim, goto error);
113
273k
114
273k
  
if (273k
n == 0 && 273k
!isl_space_is_named_or_nested(bmap->dim, type)43.7k
)
115
43.7k
    return bmap;
116
273k
117
229k
  bmap = isl_basic_map_cow(bmap);
118
229k
  if (!bmap)
119
0
    return NULL;
120
229k
121
229k
  offset = isl_basic_map_offset(bmap, type) + first;
122
229k
  left = isl_basic_map_total_dim(bmap) - (offset - 1) - n;
123
231k
  for (i = 0; 
i < bmap->n_eq231k
;
++i1.33k
)
124
1.33k
    constraint_drop_vars(bmap->eq[i]+offset, n, left);
125
229k
126
988k
  for (i = 0; 
i < bmap->n_ineq988k
;
++i758k
)
127
758k
    constraint_drop_vars(bmap->ineq[i]+offset, n, left);
128
229k
129
230k
  for (i = 0; 
i < bmap->n_div230k
;
++i711
)
130
711
    constraint_drop_vars(bmap->div[i]+1+offset, n, left);
131
229k
132
229k
  if (
type == isl_dim_div229k
)
{448
133
448
    bmap = move_divs_last(bmap, first, n);
134
448
    if (!bmap)
135
0
      goto error;
136
448
    
if (448
isl_basic_map_free_div(bmap, n) < 0448
)
137
0
      return isl_basic_map_free(bmap);
138
448
  } else
139
229k
    bmap->dim = isl_space_drop_dims(bmap->dim, type, first, n);
140
229k
  
if (229k
!bmap->dim229k
)
141
0
    goto error;
142
229k
143
229k
  
ISL_F_CLR229k
(bmap, ISL_BASIC_MAP_NORMALIZED);229k
144
229k
  bmap = isl_basic_map_simplify(bmap);
145
229k
  return isl_basic_map_finalize(bmap);
146
0
error:
147
0
  isl_basic_map_free(bmap);
148
0
  return NULL;
149
229k
}
150
151
__isl_give isl_basic_set *isl_basic_set_drop(__isl_take isl_basic_set *bset,
152
  enum isl_dim_type type, unsigned first, unsigned n)
153
168k
{
154
168k
  return bset_from_bmap(isl_basic_map_drop(bset_to_bmap(bset),
155
168k
              type, first, n));
156
168k
}
157
158
struct isl_map *isl_map_drop(struct isl_map *map,
159
  enum isl_dim_type type, unsigned first, unsigned n)
160
472
{
161
472
  int i;
162
472
163
472
  if (!map)
164
0
    goto error;
165
472
166
472
  
isl_assert472
(map->ctx, first + n <= isl_map_dim(map, type), goto error);472
167
472
168
472
  
if (472
n == 0 && 472
!isl_space_is_named_or_nested(map->dim, type)0
)
169
0
    return map;
170
472
  map = isl_map_cow(map);
171
472
  if (!map)
172
0
    goto error;
173
472
  map->dim = isl_space_drop_dims(map->dim, type, first, n);
174
472
  if (!map->dim)
175
0
    goto error;
176
472
177
1.04k
  
for (i = 0; 472
i < map->n1.04k
;
++i568
)
{568
178
568
    map->p[i] = isl_basic_map_drop(map->p[i], type, first, n);
179
568
    if (!map->p[i])
180
0
      goto error;
181
568
  }
182
472
  
ISL_F_CLR472
(map, ISL_MAP_NORMALIZED);472
183
472
184
472
  return map;
185
0
error:
186
0
  isl_map_free(map);
187
0
  return NULL;
188
472
}
189
190
struct isl_set *isl_set_drop(struct isl_set *set,
191
  enum isl_dim_type type, unsigned first, unsigned n)
192
91
{
193
91
  return set_from_map(isl_map_drop(set_to_map(set), type, first, n));
194
91
}
195
196
/*
197
 * We don't cow, as the div is assumed to be redundant.
198
 */
199
__isl_give isl_basic_map *isl_basic_map_drop_div(
200
  __isl_take isl_basic_map *bmap, unsigned div)
201
199k
{
202
199k
  int i;
203
199k
  unsigned pos;
204
199k
205
199k
  if (!bmap)
206
0
    goto error;
207
199k
208
199k
  pos = 1 + isl_space_dim(bmap->dim, isl_dim_all) + div;
209
199k
210
199k
  isl_assert(bmap->ctx, div < bmap->n_div, goto error);
211
199k
212
750k
  
for (i = 0; 199k
i < bmap->n_eq750k
;
++i551k
)
213
551k
    constraint_drop_vars(bmap->eq[i]+pos, 1, bmap->extra-div-1);
214
199k
215
1.30M
  for (i = 0; 
i < bmap->n_ineq1.30M
;
++i1.10M
)
{1.10M
216
1.10M
    if (
!1.10M
isl_int_is_zero1.10M
(bmap->ineq[i][pos]))
{15.6k
217
15.6k
      isl_basic_map_drop_inequality(bmap, i);
218
15.6k
      --i;
219
15.6k
      continue;
220
15.6k
    }
221
1.08M
    constraint_drop_vars(bmap->ineq[i]+pos, 1, bmap->extra-div-1);
222
1.08M
  }
223
199k
224
600k
  for (i = 0; 
i < bmap->n_div600k
;
++i401k
)
225
401k
    constraint_drop_vars(bmap->div[i]+1+pos, 1, bmap->extra-div-1);
226
199k
227
199k
  if (
div != bmap->n_div - 1199k
)
{12.3k
228
12.3k
    int j;
229
12.3k
    isl_int *t = bmap->div[div];
230
12.3k
231
30.6k
    for (j = div; 
j < bmap->n_div - 130.6k
;
++j18.3k
)
232
18.3k
      bmap->div[j] = bmap->div[j+1];
233
12.3k
234
12.3k
    bmap->div[bmap->n_div - 1] = t;
235
12.3k
  }
236
199k
  ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
237
199k
  if (isl_basic_map_free_div(bmap, 1) < 0)
238
0
    return isl_basic_map_free(bmap);
239
199k
240
199k
  return bmap;
241
0
error:
242
0
  isl_basic_map_free(bmap);
243
0
  return NULL;
244
199k
}
245
246
struct isl_basic_map *isl_basic_map_normalize_constraints(
247
  struct isl_basic_map *bmap)
248
1.74M
{
249
1.74M
  int i;
250
1.74M
  isl_int gcd;
251
1.74M
  unsigned total = isl_basic_map_total_dim(bmap);
252
1.74M
253
1.74M
  if (!bmap)
254
0
    return NULL;
255
1.74M
256
1.74M
  
isl_int_init1.74M
(gcd);1.74M
257
3.15M
  for (i = bmap->n_eq - 1; 
i >= 03.15M
;
--i1.41M
)
{1.41M
258
1.41M
    isl_seq_gcd(bmap->eq[i]+1, total, &gcd);
259
1.41M
    if (
isl_int_is_zero1.41M
(gcd))
{128k
260
128k
      if (
!128k
isl_int_is_zero128k
(bmap->eq[i][0]))
{1.93k
261
1.93k
        bmap = isl_basic_map_set_to_empty(bmap);
262
1.93k
        break;
263
1.93k
      }
264
126k
      isl_basic_map_drop_equality(bmap, i);
265
126k
      continue;
266
128k
    }
267
1.28M
    
if (1.28M
ISL_F_ISSET1.28M
(bmap, ISL_BASIC_MAP_RATIONAL))
268
162k
      isl_int_gcd(gcd, gcd, bmap->eq[i][0]);
269
1.28M
    if (isl_int_is_one(gcd))
270
1.25M
      continue;
271
30.9k
    
if (30.9k
!30.9k
isl_int_is_divisible_by30.9k
(bmap->eq[i][0], gcd))
{380
272
380
      bmap = isl_basic_map_set_to_empty(bmap);
273
380
      break;
274
380
    }
275
30.6k
    isl_seq_scale_down(bmap->eq[i], bmap->eq[i], gcd, 1+total);
276
30.6k
  }
277
1.74M
278
8.32M
  for (i = bmap->n_ineq - 1; 
i >= 08.32M
;
--i6.57M
)
{6.58M
279
6.58M
    isl_seq_gcd(bmap->ineq[i]+1, total, &gcd);
280
6.58M
    if (
isl_int_is_zero6.58M
(gcd))
{291k
281
291k
      if (
isl_int_is_neg291k
(bmap->ineq[i][0]))
{12.4k
282
12.4k
        bmap = isl_basic_map_set_to_empty(bmap);
283
12.4k
        break;
284
12.4k
      }
285
279k
      isl_basic_map_drop_inequality(bmap, i);
286
279k
      continue;
287
291k
    }
288
6.29M
    
if (6.29M
ISL_F_ISSET6.29M
(bmap, ISL_BASIC_MAP_RATIONAL))
289
306k
      isl_int_gcd(gcd, gcd, bmap->ineq[i][0]);
290
6.29M
    if (isl_int_is_one(gcd))
291
6.19M
      continue;
292
95.8k
    
isl_int_fdiv_q95.8k
(bmap->ineq[i][0], bmap->ineq[i][0], gcd);95.8k
293
95.8k
    isl_seq_scale_down(bmap->ineq[i]+1, bmap->ineq[i]+1, gcd, total);
294
95.8k
  }
295
1.74M
  isl_int_clear(gcd);
296
1.74M
297
1.74M
  return bmap;
298
1.74M
}
299
300
struct isl_basic_set *isl_basic_set_normalize_constraints(
301
  struct isl_basic_set *bset)
302
124k
{
303
124k
  isl_basic_map *bmap = bset_to_bmap(bset);
304
124k
  return bset_from_bmap(isl_basic_map_normalize_constraints(bmap));
305
124k
}
306
307
/* Reduce the coefficient of the variable at position "pos"
308
 * in integer division "div", such that it lies in the half-open
309
 * interval (1/2,1/2], extracting any excess value from this integer division.
310
 * "pos" is as determined by isl_basic_map_offset, i.e., pos == 0
311
 * corresponds to the constant term.
312
 *
313
 * That is, the integer division is of the form
314
 *
315
 *  floor((... + (c * d + r) * x_pos + ...)/d)
316
 *
317
 * with -d < 2 * r <= d.
318
 * Replace it by
319
 *
320
 *  floor((... + r * x_pos + ...)/d) + c * x_pos
321
 *
322
 * If 2 * ((c * d + r) % d) <= d, then c = floor((c * d + r)/d).
323
 * Otherwise, c = floor((c * d + r)/d) + 1.
324
 *
325
 * This is the same normalization that is performed by isl_aff_floor.
326
 */
327
static __isl_give isl_basic_map *reduce_coefficient_in_div(
328
  __isl_take isl_basic_map *bmap, int div, int pos)
329
6.87k
{
330
6.87k
  isl_int shift;
331
6.87k
  int add_one;
332
6.87k
333
6.87k
  isl_int_init(shift);
334
6.87k
  isl_int_fdiv_r(shift, bmap->div[div][1 + pos], bmap->div[div][0]);
335
6.87k
  isl_int_mul_ui(shift, shift, 2);
336
6.87k
  add_one = isl_int_gt(shift, bmap->div[div][0]);
337
6.87k
  isl_int_fdiv_q(shift, bmap->div[div][1 + pos], bmap->div[div][0]);
338
6.87k
  if (add_one)
339
2.81k
    isl_int_add_ui(shift, shift, 1);
340
6.87k
  isl_int_neg(shift, shift);
341
6.87k
  bmap = isl_basic_map_shift_div(bmap, div, pos, shift);
342
6.87k
  isl_int_clear(shift);
343
6.87k
344
6.87k
  return bmap;
345
6.87k
}
346
347
/* Does the coefficient of the variable at position "pos"
348
 * in integer division "div" need to be reduced?
349
 * That is, does it lie outside the half-open interval (1/2,1/2]?
350
 * The coefficient c/d lies outside this interval if abs(2 * c) >= d and
351
 * 2 * c != d.
352
 */
353
static isl_bool needs_reduction(__isl_keep isl_basic_map *bmap, int div,
354
  int pos)
355
456k
{
356
456k
  isl_bool r;
357
456k
358
456k
  if (isl_int_is_zero(bmap->div[div][1 + pos]))
359
344k
    return isl_bool_false;
360
456k
361
111k
  
isl_int_mul_ui111k
(bmap->div[div][1 + pos], bmap->div[div][1 + pos], 2);111k
362
111k
  r = isl_int_abs_ge(bmap->div[div][1 + pos], bmap->div[div][0]) &&
363
35.2k
      
!35.2k
isl_int_eq35.2k
(bmap->div[div][1 + pos], bmap->div[div][0]);
364
111k
  isl_int_divexact_ui(bmap->div[div][1 + pos],
365
111k
          bmap->div[div][1 + pos], 2);
366
111k
367
111k
  return r;
368
456k
}
369
370
/* Reduce the coefficients (including the constant term) of
371
 * integer division "div", if needed.
372
 * In particular, make sure all coefficients lie in
373
 * the half-open interval (1/2,1/2].
374
 */
375
static __isl_give isl_basic_map *reduce_div_coefficients_of_div(
376
  __isl_take isl_basic_map *bmap, int div)
377
63.1k
{
378
63.1k
  int i;
379
63.1k
  unsigned total = 1 + isl_basic_map_total_dim(bmap);
380
63.1k
381
519k
  for (i = 0; 
i < total519k
;
++i456k
)
{456k
382
456k
    isl_bool reduce;
383
456k
384
456k
    reduce = needs_reduction(bmap, div, i);
385
456k
    if (reduce < 0)
386
0
      return isl_basic_map_free(bmap);
387
456k
    
if (456k
!reduce456k
)
388
449k
      continue;
389
6.87k
    bmap = reduce_coefficient_in_div(bmap, div, i);
390
6.87k
    if (!bmap)
391
0
      break;
392
6.87k
  }
393
63.1k
394
63.1k
  return bmap;
395
63.1k
}
396
397
/* Reduce the coefficients (including the constant term) of
398
 * the known integer divisions, if needed
399
 * In particular, make sure all coefficients lie in
400
 * the half-open interval (1/2,1/2].
401
 */
402
static __isl_give isl_basic_map *reduce_div_coefficients(
403
  __isl_take isl_basic_map *bmap)
404
1.60M
{
405
1.60M
  int i;
406
1.60M
407
1.60M
  if (!bmap)
408
0
    return NULL;
409
1.60M
  
if (1.60M
bmap->n_div == 01.60M
)
410
1.47M
    return bmap;
411
1.60M
412
361k
  
for (i = 0; 130k
i < bmap->n_div361k
;
++i230k
)
{230k
413
230k
    if (isl_int_is_zero(bmap->div[i][0]))
414
166k
      continue;
415
63.1k
    bmap = reduce_div_coefficients_of_div(bmap, i);
416
63.1k
    if (!bmap)
417
0
      break;
418
63.1k
  }
419
130k
420
130k
  return bmap;
421
1.60M
}
422
423
/* Remove any common factor in numerator and denominator of the div expression,
424
 * not taking into account the constant term.
425
 * That is, if the div is of the form
426
 *
427
 *  floor((a + m f(x))/(m d))
428
 *
429
 * then replace it by
430
 *
431
 *  floor((floor(a/m) + f(x))/d)
432
 *
433
 * The difference {a/m}/d in the argument satisfies 0 <= {a/m}/d < 1/d
434
 * and can therefore not influence the result of the floor.
435
 */
436
static void normalize_div_expression(__isl_keep isl_basic_map *bmap, int div)
437
235k
{
438
235k
  unsigned total = isl_basic_map_total_dim(bmap);
439
235k
  isl_ctx *ctx = bmap->ctx;
440
235k
441
235k
  if (isl_int_is_zero(bmap->div[div][0]))
442
166k
    return;
443
68.6k
  isl_seq_gcd(bmap->div[div] + 2, total, &ctx->normalize_gcd);
444
68.6k
  isl_int_gcd(ctx->normalize_gcd, ctx->normalize_gcd, bmap->div[div][0]);
445
68.6k
  if (isl_int_is_one(ctx->normalize_gcd))
446
67.7k
    return;
447
854
  
isl_int_fdiv_q854
(bmap->div[div][1], bmap->div[div][1],854
448
854
      ctx->normalize_gcd);
449
854
  isl_int_divexact(bmap->div[div][0], bmap->div[div][0],
450
854
      ctx->normalize_gcd);
451
854
  isl_seq_scale_down(bmap->div[div] + 2, bmap->div[div] + 2,
452
854
      ctx->normalize_gcd, total);
453
854
}
454
455
/* Remove any common factor in numerator and denominator of a div expression,
456
 * not taking into account the constant term.
457
 * That is, look for any div of the form
458
 *
459
 *  floor((a + m f(x))/(m d))
460
 *
461
 * and replace it by
462
 *
463
 *  floor((floor(a/m) + f(x))/d)
464
 *
465
 * The difference {a/m}/d in the argument satisfies 0 <= {a/m}/d < 1/d
466
 * and can therefore not influence the result of the floor.
467
 */
468
static __isl_give isl_basic_map *normalize_div_expressions(
469
  __isl_take isl_basic_map *bmap)
470
1.60M
{
471
1.60M
  int i;
472
1.60M
473
1.60M
  if (!bmap)
474
0
    return NULL;
475
1.60M
  
if (1.60M
bmap->n_div == 01.60M
)
476
1.47M
    return bmap;
477
1.60M
478
361k
  
for (i = 0; 130k
i < bmap->n_div361k
;
++i230k
)
479
230k
    normalize_div_expression(bmap, i);
480
130k
481
130k
  return bmap;
482
1.60M
}
483
484
/* Assumes divs have been ordered if keep_divs is set.
485
 */
486
static void eliminate_var_using_equality(struct isl_basic_map *bmap,
487
  unsigned pos, isl_int *eq, int keep_divs, int *progress)
488
2.08M
{
489
2.08M
  unsigned total;
490
2.08M
  unsigned space_total;
491
2.08M
  int k;
492
2.08M
  int last_div;
493
2.08M
494
2.08M
  total = isl_basic_map_total_dim(bmap);
495
2.08M
  space_total = isl_space_dim(bmap->dim, isl_dim_all);
496
2.08M
  last_div = isl_seq_last_non_zero(eq + 1 + space_total, bmap->n_div);
497
13.6M
  for (k = 0; 
k < bmap->n_eq13.6M
;
++k11.5M
)
{11.5M
498
11.5M
    if (bmap->eq[k] == eq)
499
2.07M
      continue;
500
9.44M
    
if (9.44M
isl_int_is_zero9.44M
(bmap->eq[k][1+pos]))
501
8.84M
      continue;
502
605k
    
if (605k
progress605k
)
503
95.3k
      *progress = 1;
504
605k
    isl_seq_elim(bmap->eq[k], eq, 1+pos, 1+total, NULL);
505
605k
    isl_seq_normalize(bmap->ctx, bmap->eq[k], 1 + total);
506
605k
  }
507
2.08M
508
9.66M
  for (k = 0; 
k < bmap->n_ineq9.66M
;
++k7.58M
)
{7.58M
509
7.58M
    if (isl_int_is_zero(bmap->ineq[k][1+pos]))
510
7.10M
      continue;
511
480k
    
if (480k
progress480k
)
512
170k
      *progress = 1;
513
480k
    isl_seq_elim(bmap->ineq[k], eq, 1+pos, 1+total, NULL);
514
480k
    isl_seq_normalize(bmap->ctx, bmap->ineq[k], 1 + total);
515
480k
    ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
516
480k
  }
517
2.08M
518
2.52M
  for (k = 0; 
k < bmap->n_div2.52M
;
++k440k
)
{440k
519
440k
    if (isl_int_is_zero(bmap->div[k][0]))
520
298k
      continue;
521
141k
    
if (141k
isl_int_is_zero141k
(bmap->div[k][1+1+pos]))
522
136k
      continue;
523
5.46k
    
if (5.46k
progress5.46k
)
524
2.56k
      *progress = 1;
525
5.46k
    /* We need to be careful about circular definitions,
526
5.46k
     * so for now we just remove the definition of div k
527
5.46k
     * if the equality contains any divs.
528
5.46k
     * If keep_divs is set, then the divs have been ordered
529
5.46k
     * and we can keep the definition as long as the result
530
5.46k
     * is still ordered.
531
5.46k
     */
532
5.46k
    if (
last_div == -1 || 5.46k
(keep_divs && 2.16k
last_div < k2.16k
))
{5.46k
533
5.46k
      isl_seq_elim(bmap->div[k]+1, eq,
534
5.46k
          1+pos, 1+total, &bmap->div[k][0]);
535
5.46k
      normalize_div_expression(bmap, k);
536
5.46k
    } else
537
0
      isl_seq_clr(bmap->div[k], 1 + total);
538
5.46k
    ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
539
5.46k
  }
540
2.08M
}
541
542
/* Assumes divs have been ordered if keep_divs is set.
543
 */
544
static __isl_give isl_basic_map *eliminate_div(__isl_take isl_basic_map *bmap,
545
  isl_int *eq, unsigned div, int keep_divs)
546
118k
{
547
118k
  unsigned pos = isl_space_dim(bmap->dim, isl_dim_all) + div;
548
118k
549
118k
  eliminate_var_using_equality(bmap, pos, eq, keep_divs, NULL);
550
118k
551
118k
  bmap = isl_basic_map_drop_div(bmap, div);
552
118k
553
118k
  return bmap;
554
118k
}
555
556
/* Check if elimination of div "div" using equality "eq" would not
557
 * result in a div depending on a later div.
558
 */
559
static isl_bool ok_to_eliminate_div(struct isl_basic_map *bmap, isl_int *eq,
560
  unsigned div)
561
119k
{
562
119k
  int k;
563
119k
  int last_div;
564
119k
  unsigned space_total = isl_space_dim(bmap->dim, isl_dim_all);
565
119k
  unsigned pos = space_total + div;
566
119k
567
119k
  last_div = isl_seq_last_non_zero(eq + 1 + space_total, bmap->n_div);
568
119k
  if (
last_div < 0 || 119k
last_div <= div119k
)
569
115k
    return isl_bool_true;
570
119k
571
13.7k
  
for (k = 0; 4.02k
k <= last_div13.7k
;
++k9.67k
)
{13.3k
572
13.3k
    if (isl_int_is_zero(bmap->div[k][0]))
573
5.30k
      continue;
574
8.06k
    
if (8.06k
!8.06k
isl_int_is_zero8.06k
(bmap->div[k][1 + 1 + pos]))
575
3.68k
      return isl_bool_false;
576
8.06k
  }
577
4.02k
578
340
  return isl_bool_true;
579
4.02k
}
580
581
/* Eliminate divs based on equalities
582
 */
583
static struct isl_basic_map *eliminate_divs_eq(
584
    struct isl_basic_map *bmap, int *progress)
585
1.66M
{
586
1.66M
  int d;
587
1.66M
  int i;
588
1.66M
  int modified = 0;
589
1.66M
  unsigned off;
590
1.66M
591
1.66M
  bmap = isl_basic_map_order_divs(bmap);
592
1.66M
593
1.66M
  if (!bmap)
594
0
    return NULL;
595
1.66M
596
1.66M
  off = 1 + isl_space_dim(bmap->dim, isl_dim_all);
597
1.66M
598
1.90M
  for (d = bmap->n_div - 1; 
d >= 01.90M
;
--d242k
)
{242k
599
502k
    for (i = 0; 
i < bmap->n_eq502k
;
++i260k
)
{376k
600
376k
      isl_bool ok;
601
376k
602
376k
      if (
!376k
isl_int_is_one376k
(bmap->eq[i][off + d]) &&
603
280k
          
!280k
isl_int_is_negone280k
(bmap->eq[i][off + d]))
604
256k
        continue;
605
119k
      ok = ok_to_eliminate_div(bmap, bmap->eq[i], d);
606
119k
      if (ok < 0)
607
0
        return isl_basic_map_free(bmap);
608
119k
      
if (119k
!ok119k
)
609
3.68k
        continue;
610
115k
      modified = 1;
611
115k
      *progress = 1;
612
115k
      bmap = eliminate_div(bmap, bmap->eq[i], d, 1);
613
115k
      if (isl_basic_map_drop_equality(bmap, i) < 0)
614
0
        return isl_basic_map_free(bmap);
615
115k
      break;
616
115k
    }
617
242k
  }
618
1.66M
  
if (1.66M
modified1.66M
)
619
60.1k
    return eliminate_divs_eq(bmap, progress);
620
1.60M
  return bmap;
621
1.66M
}
622
623
/* Elimininate divs based on inequalities
624
 */
625
static struct isl_basic_map *eliminate_divs_ineq(
626
    struct isl_basic_map *bmap, int *progress)
627
1.60M
{
628
1.60M
  int d;
629
1.60M
  int i;
630
1.60M
  unsigned off;
631
1.60M
  struct isl_ctx *ctx;
632
1.60M
633
1.60M
  if (!bmap)
634
0
    return NULL;
635
1.60M
636
1.60M
  ctx = bmap->ctx;
637
1.60M
  off = 1 + isl_space_dim(bmap->dim, isl_dim_all);
638
1.60M
639
1.71M
  for (d = bmap->n_div - 1; 
d >= 01.71M
;
--d111k
)
{111k
640
167k
    for (i = 0; 
i < bmap->n_eq167k
;
++i55.5k
)
641
89.4k
      
if (89.4k
!89.4k
isl_int_is_zero89.4k
(bmap->eq[i][off + d]))
642
33.9k
        break;
643
111k
    if (i < bmap->n_eq)
644
33.9k
      continue;
645
301k
    
for (i = 0; 77.7k
i < bmap->n_ineq301k
;
++i223k
)
646
256k
      
if (256k
isl_int_abs_gt256k
(bmap->ineq[i][off + d], ctx->one))
647
32.1k
        break;
648
77.7k
    if (i < bmap->n_ineq)
649
32.1k
      continue;
650
45.6k
    *progress = 1;
651
45.6k
    bmap = isl_basic_map_eliminate_vars(bmap, (off-1)+d, 1);
652
45.6k
    if (
!bmap || 45.6k
ISL_F_ISSET45.6k
(bmap, ISL_BASIC_MAP_EMPTY))
653
67
      break;
654
45.5k
    bmap = isl_basic_map_drop_div(bmap, d);
655
45.5k
    if (!bmap)
656
0
      break;
657
45.5k
  }
658
1.60M
  return bmap;
659
1.60M
}
660
661
/* Does the equality constraint at position "eq" in "bmap" involve
662
 * any local variables in the range [first, first + n)
663
 * that are not marked as having an explicit representation?
664
 */
665
static isl_bool bmap_eq_involves_unknown_divs(__isl_keep isl_basic_map *bmap,
666
  int eq, unsigned first, unsigned n)
667
1.67k
{
668
1.67k
  unsigned o_div;
669
1.67k
  int i;
670
1.67k
671
1.67k
  if (!bmap)
672
0
    return isl_bool_error;
673
1.67k
674
1.67k
  o_div = isl_basic_map_offset(bmap, isl_dim_div);
675
2.17k
  for (i = 0; 
i < n2.17k
;
++i500
)
{865
676
865
    isl_bool unknown;
677
865
678
865
    if (isl_int_is_zero(bmap->eq[eq][o_div + first + i]))
679
487
      continue;
680
378
    unknown = isl_basic_map_div_is_marked_unknown(bmap, first + i);
681
378
    if (unknown < 0)
682
0
      return isl_bool_error;
683
378
    
if (378
unknown378
)
684
365
      return isl_bool_true;
685
378
  }
686
1.67k
687
1.30k
  return isl_bool_false;
688
1.67k
}
689
690
/* The last local variable involved in the equality constraint
691
 * at position "eq" in "bmap" is the local variable at position "div".
692
 * It can therefore be used to extract an explicit representation
693
 * for that variable.
694
 * Do so unless the local variable already has an explicit representation or
695
 * the explicit representation would involve any other local variables
696
 * that in turn do not have an explicit representation.
697
 * An equality constraint involving local variables without an explicit
698
 * representation can be used in isl_basic_map_drop_redundant_divs
699
 * to separate out an independent local variable.  Introducing
700
 * an explicit representation here would block this transformation,
701
 * while the partial explicit representation in itself is not very useful.
702
 * Set *progress if anything is changed.
703
 *
704
 * The equality constraint is of the form
705
 *
706
 *  f(x) + n e >= 0
707
 *
708
 * with n a positive number.  The explicit representation derived from
709
 * this constraint is
710
 *
711
 *  floor((-f(x))/n)
712
 */
713
static __isl_give isl_basic_map *set_div_from_eq(__isl_take isl_basic_map *bmap,
714
  int div, int eq, int *progress)
715
38.8k
{
716
38.8k
  unsigned total, o_div;
717
38.8k
  isl_bool involves;
718
38.8k
719
38.8k
  if (!bmap)
720
0
    return NULL;
721
38.8k
722
38.8k
  
if (38.8k
!38.8k
isl_int_is_zero38.8k
(bmap->div[div][0]))
723
37.1k
    return bmap;
724
38.8k
725
1.67k
  involves = bmap_eq_involves_unknown_divs(bmap, eq, 0, div);
726
1.67k
  if (involves < 0)
727
0
    return isl_basic_map_free(bmap);
728
1.67k
  
if (1.67k
involves1.67k
)
729
365
    return bmap;
730
1.67k
731
1.30k
  total = isl_basic_map_dim(bmap, isl_dim_all);
732
1.30k
  o_div = isl_basic_map_offset(bmap, isl_dim_div);
733
1.30k
  isl_seq_neg(bmap->div[div] + 1, bmap->eq[eq], 1 + total);
734
1.30k
  isl_int_set_si(bmap->div[div][1 + o_div + div], 0);
735
1.30k
  isl_int_set(bmap->div[div][0], bmap->eq[eq][o_div + div]);
736
1.30k
  if (progress)
737
1.08k
    *progress = 1;
738
1.30k
  ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
739
1.30k
740
1.30k
  return bmap;
741
1.67k
}
742
743
struct isl_basic_map *isl_basic_map_gauss(
744
  struct isl_basic_map *bmap, int *progress)
745
2.74M
{
746
2.74M
  int k;
747
2.74M
  int done;
748
2.74M
  int last_var;
749
2.74M
  unsigned total_var;
750
2.74M
  unsigned total;
751
2.74M
752
2.74M
  bmap = isl_basic_map_order_divs(bmap);
753
2.74M
754
2.74M
  if (!bmap)
755
0
    return NULL;
756
2.74M
757
2.74M
  total = isl_basic_map_total_dim(bmap);
758
2.74M
  total_var = total - bmap->n_div;
759
2.74M
760
2.74M
  last_var = total - 1;
761
4.70M
  for (done = 0; 
done < bmap->n_eq4.70M
;
++done1.95M
)
{2.10M
762
4.36M
    for (; 
last_var >= 04.36M
;
--last_var2.25M
)
{4.21M
763
13.9M
      for (k = done; 
k < bmap->n_eq13.9M
;
++k9.72M
)
764
11.6M
        
if (11.6M
!11.6M
isl_int_is_zero11.6M
(bmap->eq[k][1+last_var]))
765
1.95M
          break;
766
4.21M
      if (k < bmap->n_eq)
767
1.95M
        break;
768
4.21M
    }
769
2.10M
    if (last_var < 0)
770
148k
      break;
771
1.95M
    
if (1.95M
k != done1.95M
)
772
362k
      swap_equality(bmap, k, done);
773
1.95M
    if (isl_int_is_neg(bmap->eq[done][1+last_var]))
774
233k
      isl_seq_neg(bmap->eq[done], bmap->eq[done], 1+total);
775
1.95M
776
1.95M
    eliminate_var_using_equality(bmap, last_var, bmap->eq[done], 1,
777
1.95M
            progress);
778
1.95M
779
1.95M
    if (last_var >= total_var)
780
38.8k
      bmap = set_div_from_eq(bmap, last_var - total_var,
781
38.8k
            done, progress);
782
1.95M
    if (!bmap)
783
0
      return NULL;
784
1.95M
  }
785
2.74M
  
if (2.74M
done == bmap->n_eq2.74M
)
786
2.59M
    return bmap;
787
452k
  
for (k = done; 148k
k < bmap->n_eq452k
;
++k303k
)
{331k
788
331k
    if (isl_int_is_zero(bmap->eq[k][0]))
789
303k
      continue;
790
27.8k
    return isl_basic_map_set_to_empty(bmap);
791
331k
  }
792
120k
  isl_basic_map_free_equality(bmap, bmap->n_eq-done);
793
120k
  return bmap;
794
148k
}
795
796
struct isl_basic_set *isl_basic_set_gauss(
797
  struct isl_basic_set *bset, int *progress)
798
82.9k
{
799
82.9k
  return bset_from_bmap(isl_basic_map_gauss(bset_to_bmap(bset),
800
82.9k
              progress));
801
82.9k
}
802
803
804
static unsigned int round_up(unsigned int v)
805
1.09M
{
806
1.09M
  int old_v = v;
807
1.09M
808
3.32M
  while (
v3.32M
)
{2.22M
809
2.22M
    old_v = v;
810
2.22M
    v ^= v & -v;
811
2.22M
  }
812
1.09M
  return old_v << 1;
813
1.09M
}
814
815
/* Hash table of inequalities in a basic map.
816
 * "index" is an array of addresses of inequalities in the basic map, some
817
 * of which are NULL.  The inequalities are hashed on the coefficients
818
 * except the constant term.
819
 * "size" is the number of elements in the array and is always a power of two
820
 * "bits" is the number of bits need to represent an index into the array.
821
 * "total" is the total dimension of the basic map.
822
 */
823
struct isl_constraint_index {
824
  unsigned int size;
825
  int bits;
826
  isl_int ***index;
827
  unsigned total;
828
};
829
830
/* Fill in the "ci" data structure for holding the inequalities of "bmap".
831
 */
832
static isl_stat create_constraint_index(struct isl_constraint_index *ci,
833
  __isl_keep isl_basic_map *bmap)
834
1.08M
{
835
1.08M
  isl_ctx *ctx;
836
1.08M
837
1.08M
  ci->index = NULL;
838
1.08M
  if (!bmap)
839
0
    return isl_stat_error;
840
1.08M
  ci->total = isl_basic_set_total_dim(bmap);
841
1.08M
  if (bmap->n_ineq == 0)
842
0
    return isl_stat_ok;
843
1.08M
  ci->size = round_up(4 * (bmap->n_ineq + 1) / 3 - 1);
844
1.08M
  ci->bits = ffs(ci->size) - 1;
845
1.08M
  ctx = isl_basic_map_get_ctx(bmap);
846
1.08M
  ci->index = isl_calloc_array(ctx, isl_int **, ci->size);
847
1.08M
  if (!ci->index)
848
0
    return isl_stat_error;
849
1.08M
850
1.08M
  return isl_stat_ok;
851
1.08M
}
852
853
/* Free the memory allocated by create_constraint_index.
854
 */
855
static void constraint_index_free(struct isl_constraint_index *ci)
856
1.08M
{
857
1.08M
  free(ci->index);
858
1.08M
}
859
860
/* Return the position in ci->index that contains the address of
861
 * an inequality that is equal to *ineq up to the constant term,
862
 * provided this address is not identical to "ineq".
863
 * If there is no such inequality, then return the position where
864
 * such an inequality should be inserted.
865
 */
866
static int hash_index_ineq(struct isl_constraint_index *ci, isl_int **ineq)
867
9.23M
{
868
9.23M
  int h;
869
9.23M
  uint32_t hash = isl_seq_get_hash_bits((*ineq) + 1, ci->total, ci->bits);
870
12.5M
  for (h = hash; 
ci->index[h]12.5M
;
h = (h+1) % ci->size3.27M
)
871
6.68M
    
if (6.68M
ineq != ci->index[h] &&6.68M
872
6.47M
        isl_seq_eq((*ineq) + 1, ci->index[h][0]+1, ci->total))
873
3.40M
      break;
874
9.23M
  return h;
875
9.23M
}
876
877
/* Return the position in ci->index that contains the address of
878
 * an inequality that is equal to the k'th inequality of "bmap"
879
 * up to the constant term, provided it does not point to the very
880
 * same inequality.
881
 * If there is no such inequality, then return the position where
882
 * such an inequality should be inserted.
883
 */
884
static int hash_index(struct isl_constraint_index *ci,
885
  __isl_keep isl_basic_map *bmap, int k)
886
9.21M
{
887
9.21M
  return hash_index_ineq(ci, &bmap->ineq[k]);
888
9.21M
}
889
890
static int set_hash_index(struct isl_constraint_index *ci,
891
        struct isl_basic_set *bset, int k)
892
31.1k
{
893
31.1k
  return hash_index(ci, bset, k);
894
31.1k
}
895
896
/* Fill in the "ci" data structure with the inequalities of "bset".
897
 */
898
static isl_stat setup_constraint_index(struct isl_constraint_index *ci,
899
  __isl_keep isl_basic_set *bset)
900
8.94k
{
901
8.94k
  int k, h;
902
8.94k
903
8.94k
  if (create_constraint_index(ci, bset) < 0)
904
0
    return isl_stat_error;
905
8.94k
906
40.1k
  
for (k = 0; 8.94k
k < bset->n_ineq40.1k
;
++k31.1k
)
{31.1k
907
31.1k
    h = set_hash_index(ci, bset, k);
908
31.1k
    ci->index[h] = &bset->ineq[k];
909
31.1k
  }
910
8.94k
911
8.94k
  return isl_stat_ok;
912
8.94k
}
913
914
/* Is the inequality ineq (obviously) redundant with respect
915
 * to the constraints in "ci"?
916
 *
917
 * Look for an inequality in "ci" with the same coefficients and then
918
 * check if the contant term of "ineq" is greater than or equal
919
 * to the constant term of that inequality.  If so, "ineq" is clearly
920
 * redundant.
921
 *
922
 * Note that hash_index_ineq ignores a stored constraint if it has
923
 * the same address as the passed inequality.  It is ok to pass
924
 * the address of a local variable here since it will never be
925
 * the same as the address of a constraint in "ci".
926
 */
927
static isl_bool constraint_index_is_redundant(struct isl_constraint_index *ci,
928
  isl_int *ineq)
929
23.5k
{
930
23.5k
  int h;
931
23.5k
932
23.5k
  h = hash_index_ineq(ci, &ineq);
933
23.5k
  if (!ci->index[h])
934
10.8k
    return isl_bool_false;
935
12.6k
  
return 12.6k
isl_int_ge12.6k
(ineq[0], (*ci->index[h])[0]);
936
23.5k
}
937
938
/* If we can eliminate more than one div, then we need to make
939
 * sure we do it from last div to first div, in order not to
940
 * change the position of the other divs that still need to
941
 * be removed.
942
 */
943
static struct isl_basic_map *remove_duplicate_divs(
944
  struct isl_basic_map *bmap, int *progress)
945
1.60M
{
946
1.60M
  unsigned int size;
947
1.60M
  int *index;
948
1.60M
  int *elim_for;
949
1.60M
  int k, l, h;
950
1.60M
  int bits;
951
1.60M
  struct isl_blk eq;
952
1.60M
  unsigned total_var;
953
1.60M
  unsigned total;
954
1.60M
  struct isl_ctx *ctx;
955
1.60M
956
1.60M
  bmap = isl_basic_map_order_divs(bmap);
957
1.60M
  if (
!bmap || 1.60M
bmap->n_div <= 11.60M
)
958
1.55M
    return bmap;
959
1.60M
960
48.3k
  total_var = isl_space_dim(bmap->dim, isl_dim_all);
961
48.3k
  total = total_var + bmap->n_div;
962
48.3k
963
48.3k
  ctx = bmap->ctx;
964
157k
  for (k = bmap->n_div - 1; 
k >= 0157k
;
--k109k
)
965
122k
    
if (122k
!122k
isl_int_is_zero122k
(bmap->div[k][0]))
966
13.5k
      break;
967
48.3k
  if (k <= 0)
968
35.2k
    return bmap;
969
48.3k
970
13.0k
  size = round_up(4 * bmap->n_div / 3 - 1);
971
13.0k
  if (size == 0)
972
0
    return bmap;
973
13.0k
  
elim_for = 13.0k
isl_calloc_array13.0k
(ctx, int, bmap->n_div);
974
13.0k
  bits = ffs(size) - 1;
975
13.0k
  index = isl_calloc_array(ctx, int, size);
976
13.0k
  if (
!elim_for || 13.0k
!index13.0k
)
977
0
    goto out;
978
13.0k
  eq = isl_blk_alloc(ctx, 1+total);
979
13.0k
  if (isl_blk_is_error(eq))
980
0
    goto out;
981
13.0k
982
13.0k
  isl_seq_clr(eq.data, 1+total);
983
13.0k
  index[isl_seq_get_hash_bits(bmap->div[k], 2+total, bits)] = k + 1;
984
37.9k
  for (--k; 
k >= 037.9k
;
--k24.8k
)
{24.8k
985
24.8k
    uint32_t hash;
986
24.8k
987
24.8k
    if (isl_int_is_zero(bmap->div[k][0]))
988
5.97k
      continue;
989
24.8k
990
18.8k
    hash = isl_seq_get_hash_bits(bmap->div[k], 2+total, bits);
991
24.7k
    for (h = hash; 
index[h]24.7k
;
h = (h+1) % size5.91k
)
992
8.72k
      
if (8.72k
isl_seq_eq(bmap->div[k],8.72k
993
8.72k
               bmap->div[index[h]-1], 2+total))
994
2.81k
        break;
995
18.8k
    if (
index[h]18.8k
)
{2.81k
996
2.81k
      *progress = 1;
997
2.81k
      l = index[h] - 1;
998
2.81k
      elim_for[l] = k + 1;
999
2.81k
    }
1000
18.8k
    index[h] = k+1;
1001
18.8k
  }
1002
52.0k
  for (l = bmap->n_div - 1; 
l >= 052.0k
;
--l38.9k
)
{38.9k
1003
38.9k
    if (!elim_for[l])
1004
36.1k
      continue;
1005
2.81k
    k = elim_for[l] - 1;
1006
2.81k
    isl_int_set_si(eq.data[1+total_var+k], -1);
1007
2.81k
    isl_int_set_si(eq.data[1+total_var+l], 1);
1008
2.81k
    bmap = eliminate_div(bmap, eq.data, l, 1);
1009
2.81k
    if (!bmap)
1010
0
      break;
1011
2.81k
    
isl_int_set_si2.81k
(eq.data[1+total_var+k], 0);2.81k
1012
2.81k
    isl_int_set_si(eq.data[1+total_var+l], 0);
1013
2.81k
  }
1014
13.0k
1015
13.0k
  isl_blk_free(ctx, eq);
1016
13.0k
out:
1017
13.0k
  free(index);
1018
13.0k
  free(elim_for);
1019
13.0k
  return bmap;
1020
13.0k
}
1021
1022
static int n_pure_div_eq(struct isl_basic_map *bmap)
1023
14.2k
{
1024
14.2k
  int i, j;
1025
14.2k
  unsigned total;
1026
14.2k
1027
14.2k
  total = isl_space_dim(bmap->dim, isl_dim_all);
1028
25.1k
  for (i = 0, j = bmap->n_div-1; 
i < bmap->n_eq25.1k
;
++i10.8k
)
{21.0k
1029
33.7k
    while (
j >= 0 && 33.7k
isl_int_is_zero25.7k
(bmap->eq[i][1 + total + j]))
1030
12.6k
      --j;
1031
21.0k
    if (j < 0)
1032
7.97k
      break;
1033
13.0k
    
if (13.0k
isl_seq_first_non_zero(bmap->eq[i] + 1 + total, j) != -113.0k
)
1034
2.21k
      return 0;
1035
13.0k
  }
1036
12.0k
  return i;
1037
14.2k
}
1038
1039
/* Normalize divs that appear in equalities.
1040
 *
1041
 * In particular, we assume that bmap contains some equalities
1042
 * of the form
1043
 *
1044
 *  a x = m * e_i
1045
 *
1046
 * and we want to replace the set of e_i by a minimal set and
1047
 * such that the new e_i have a canonical representation in terms
1048
 * of the vector x.
1049
 * If any of the equalities involves more than one divs, then
1050
 * we currently simply bail out.
1051
 *
1052
 * Let us first additionally assume that all equalities involve
1053
 * a div.  The equalities then express modulo constraints on the
1054
 * remaining variables and we can use "parameter compression"
1055
 * to find a minimal set of constraints.  The result is a transformation
1056
 *
1057
 *  x = T(x') = x_0 + G x'
1058
 *
1059
 * with G a lower-triangular matrix with all elements below the diagonal
1060
 * non-negative and smaller than the diagonal element on the same row.
1061
 * We first normalize x_0 by making the same property hold in the affine
1062
 * T matrix.
1063
 * The rows i of G with a 1 on the diagonal do not impose any modulo
1064
 * constraint and simply express x_i = x'_i.
1065
 * For each of the remaining rows i, we introduce a div and a corresponding
1066
 * equality.  In particular
1067
 *
1068
 *  g_ii e_j = x_i - g_i(x')
1069
 *
1070
 * where each x'_k is replaced either by x_k (if g_kk = 1) or the
1071
 * corresponding div (if g_kk != 1).
1072
 *
1073
 * If there are any equalities not involving any div, then we
1074
 * first apply a variable compression on the variables x:
1075
 *
1076
 *  x = C x'' x'' = C_2 x
1077
 *
1078
 * and perform the above parameter compression on A C instead of on A.
1079
 * The resulting compression is then of the form
1080
 *
1081
 *  x'' = T(x') = x_0 + G x'
1082
 *
1083
 * and in constructing the new divs and the corresponding equalities,
1084
 * we have to replace each x'', i.e., the x'_k with (g_kk = 1),
1085
 * by the corresponding row from C_2.
1086
 */
1087
static struct isl_basic_map *normalize_divs(
1088
  struct isl_basic_map *bmap, int *progress)
1089
1.60M
{
1090
1.60M
  int i, j, k;
1091
1.60M
  int total;
1092
1.60M
  int div_eq;
1093
1.60M
  struct isl_mat *B;
1094
1.60M
  struct isl_vec *d;
1095
1.60M
  struct isl_mat *T = NULL;
1096
1.60M
  struct isl_mat *C = NULL;
1097
1.60M
  struct isl_mat *C2 = NULL;
1098
1.60M
  isl_int v;
1099
1.60M
  int *pos = NULL;
1100
1.60M
  int dropped, needed;
1101
1.60M
1102
1.60M
  if (!bmap)
1103
0
    return NULL;
1104
1.60M
1105
1.60M
  
if (1.60M
bmap->n_div == 01.60M
)
1106
1.55M
    return bmap;
1107
1.60M
1108
46.7k
  
if (46.7k
bmap->n_eq == 046.7k
)
1109
18.3k
    return bmap;
1110
46.7k
1111
28.3k
  
if (28.3k
ISL_F_ISSET28.3k
(bmap, ISL_BASIC_MAP_NORMALIZED_DIVS))
1112
14.0k
    return bmap;
1113
28.3k
1114
14.2k
  total = isl_space_dim(bmap->dim, isl_dim_all);
1115
14.2k
  div_eq = n_pure_div_eq(bmap);
1116
14.2k
  if (div_eq == 0)
1117
5.87k
    return bmap;
1118
14.2k
1119
8.41k
  
if (8.41k
div_eq < bmap->n_eq8.41k
)
{4.30k
1120
4.30k
    B = isl_mat_sub_alloc6(bmap->ctx, bmap->eq, div_eq,
1121
4.30k
          bmap->n_eq - div_eq, 0, 1 + total);
1122
4.30k
    C = isl_mat_variable_compression(B, &C2);
1123
4.30k
    if (
!C || 4.30k
!C24.30k
)
1124
0
      goto error;
1125
4.30k
    
if (4.30k
C->n_col == 04.30k
)
{4
1126
4
      bmap = isl_basic_map_set_to_empty(bmap);
1127
4
      isl_mat_free(C);
1128
4
      isl_mat_free(C2);
1129
4
      goto done;
1130
4
    }
1131
4.30k
  }
1132
8.41k
1133
8.41k
  d = isl_vec_alloc(bmap->ctx, div_eq);
1134
8.41k
  if (!d)
1135
0
    goto error;
1136
19.2k
  
for (i = 0, j = bmap->n_div-1; 8.41k
i < div_eq19.2k
;
++i10.7k
)
{10.7k
1137
13.8k
    while (
j >= 0 && 13.8k
isl_int_is_zero13.8k
(bmap->eq[i][1 + total + j]))
1138
3.05k
      --j;
1139
10.7k
    isl_int_set(d->block.data[i], bmap->eq[i][1 + total + j]);
1140
10.7k
  }
1141
8.41k
  B = isl_mat_sub_alloc6(bmap->ctx, bmap->eq, 0, div_eq, 0, 1 + total);
1142
8.41k
1143
8.41k
  if (
C8.41k
)
{4.30k
1144
4.30k
    B = isl_mat_product(B, C);
1145
4.30k
    C = NULL;
1146
4.30k
  }
1147
8.41k
1148
8.41k
  T = isl_mat_parameter_compression(B, d);
1149
8.41k
  if (!T)
1150
0
    goto error;
1151
8.41k
  
if (8.41k
T->n_col == 08.41k
)
{587
1152
587
    bmap = isl_basic_map_set_to_empty(bmap);
1153
587
    isl_mat_free(C2);
1154
587
    isl_mat_free(T);
1155
587
    goto done;
1156
587
  }
1157
7.82k
  
isl_int_init7.82k
(v);7.82k
1158
30.5k
  for (i = 0; 
i < T->n_row - 130.5k
;
++i22.7k
)
{22.7k
1159
22.7k
    isl_int_fdiv_q(v, T->row[1 + i][0], T->row[1 + i][1 + i]);
1160
22.7k
    if (isl_int_is_zero(v))
1161
20.7k
      continue;
1162
2.03k
    isl_mat_col_submul(T, 0, v, 1 + i);
1163
2.03k
  }
1164
7.82k
  isl_int_clear(v);
1165
7.82k
  pos = isl_alloc_array(bmap->ctx, int, T->n_row);
1166
7.82k
  if (!pos)
1167
0
    goto error;
1168
7.82k
  /* We have to be careful because dropping equalities may reorder them */
1169
7.82k
  dropped = 0;
1170
18.9k
  for (j = bmap->n_div - 1; 
j >= 018.9k
;
--j11.1k
)
{11.1k
1171
14.6k
    for (i = 0; 
i < bmap->n_eq14.6k
;
++i3.47k
)
1172
13.1k
      
if (13.1k
!13.1k
isl_int_is_zero13.1k
(bmap->eq[i][1 + total + j]))
1173
9.71k
        break;
1174
11.1k
    if (
i < bmap->n_eq11.1k
)
{9.71k
1175
9.71k
      bmap = isl_basic_map_drop_div(bmap, j);
1176
9.71k
      isl_basic_map_drop_equality(bmap, i);
1177
9.71k
      ++dropped;
1178
9.71k
    }
1179
11.1k
  }
1180
7.82k
  pos[0] = 0;
1181
7.82k
  needed = 0;
1182
30.5k
  for (i = 1; 
i < T->n_row30.5k
;
++i22.7k
)
{22.7k
1183
22.7k
    if (isl_int_is_one(T->row[i][i]))
1184
14.1k
      pos[i] = i;
1185
22.7k
    else
1186
8.60k
      needed++;
1187
22.7k
  }
1188
7.82k
  if (
needed > dropped7.82k
)
{13
1189
13
    bmap = isl_basic_map_extend_space(bmap, isl_space_copy(bmap->dim),
1190
13
        needed, needed, 0);
1191
13
    if (!bmap)
1192
0
      goto error;
1193
13
  }
1194
30.5k
  
for (i = 1; 7.82k
i < T->n_row30.5k
;
++i22.7k
)
{22.7k
1195
22.7k
    if (isl_int_is_one(T->row[i][i]))
1196
14.1k
      continue;
1197
8.60k
    k = isl_basic_map_alloc_div(bmap);
1198
8.60k
    pos[i] = 1 + total + k;
1199
8.60k
    isl_seq_clr(bmap->div[k] + 1, 1 + total + bmap->n_div);
1200
8.60k
    isl_int_set(bmap->div[k][0], T->row[i][i]);
1201
8.60k
    if (C2)
1202
4.19k
      isl_seq_cpy(bmap->div[k] + 1, C2->row[i], 1 + total);
1203
8.60k
    else
1204
4.41k
      isl_int_set_si(bmap->div[k][1 + i], 1);
1205
30.9k
    for (j = 0; 
j < i30.9k
;
++j22.3k
)
{22.3k
1206
22.3k
      if (isl_int_is_zero(T->row[i][j]))
1207
16.6k
        continue;
1208
5.71k
      
if (5.71k
pos[j] < T->n_row && 5.71k
C25.63k
)
1209
2.74k
        isl_seq_submul(bmap->div[k] + 1, T->row[i][j],
1210
2.74k
            C2->row[pos[j]], 1 + total);
1211
5.71k
      else
1212
2.96k
        isl_int_neg(bmap->div[k][1 + pos[j]],
1213
5.71k
                T->row[i][j]);
1214
5.71k
    }
1215
8.60k
    j = isl_basic_map_alloc_equality(bmap);
1216
8.60k
    isl_seq_neg(bmap->eq[j], bmap->div[k]+1, 1+total+bmap->n_div);
1217
8.60k
    isl_int_set(bmap->eq[j][pos[i]], bmap->div[k][0]);
1218
8.60k
  }
1219
7.82k
  free(pos);
1220
7.82k
  isl_mat_free(C2);
1221
7.82k
  isl_mat_free(T);
1222
7.82k
1223
7.82k
  if (progress)
1224
7.82k
    *progress = 1;
1225
8.41k
done:
1226
8.41k
  ISL_F_SET(bmap, ISL_BASIC_MAP_NORMALIZED_DIVS);
1227
8.41k
1228
8.41k
  return bmap;
1229
0
error:
1230
0
  free(pos);
1231
0
  isl_mat_free(C);
1232
0
  isl_mat_free(C2);
1233
0
  isl_mat_free(T);
1234
0
  return bmap;
1235
7.82k
}
1236
1237
static struct isl_basic_map *set_div_from_lower_bound(
1238
  struct isl_basic_map *bmap, int div, int ineq)
1239
816
{
1240
816
  unsigned total = 1 + isl_space_dim(bmap->dim, isl_dim_all);
1241
816
1242
816
  isl_seq_neg(bmap->div[div] + 1, bmap->ineq[ineq], total + bmap->n_div);
1243
816
  isl_int_set(bmap->div[div][0], bmap->ineq[ineq][total + div]);
1244
816
  isl_int_add(bmap->div[div][1], bmap->div[div][1], bmap->div[div][0]);
1245
816
  isl_int_sub_ui(bmap->div[div][1], bmap->div[div][1], 1);
1246
816
  isl_int_set_si(bmap->div[div][1 + total + div], 0);
1247
816
1248
816
  return bmap;
1249
816
}
1250
1251
/* Check whether it is ok to define a div based on an inequality.
1252
 * To avoid the introduction of circular definitions of divs, we
1253
 * do not allow such a definition if the resulting expression would refer to
1254
 * any other undefined divs or if any known div is defined in
1255
 * terms of the unknown div.
1256
 */
1257
static isl_bool ok_to_set_div_from_bound(struct isl_basic_map *bmap,
1258
  int div, int ineq)
1259
1.54k
{
1260
1.54k
  int j;
1261
1.54k
  unsigned total = 1 + isl_space_dim(bmap->dim, isl_dim_all);
1262
1.54k
1263
1.54k
  /* Not defined in terms of unknown divs */
1264
5.99k
  for (j = 0; 
j < bmap->n_div5.99k
;
++j4.45k
)
{5.06k
1265
5.06k
    if (div == j)
1266
1.22k
      continue;
1267
3.84k
    
if (3.84k
isl_int_is_zero3.84k
(bmap->ineq[ineq][total + j]))
1268
3.15k
      continue;
1269
686
    
if (686
isl_int_is_zero686
(bmap->div[j][0]))
1270
617
      return isl_bool_false;
1271
686
  }
1272
1.54k
1273
1.54k
  /* No other div defined in terms of this one => avoid loops */
1274
4.11k
  
for (j = 0; 923
j < bmap->n_div4.11k
;
++j3.19k
)
{3.30k
1275
3.30k
    if (div == j)
1276
923
      continue;
1277
2.38k
    
if (2.38k
isl_int_is_zero2.38k
(bmap->div[j][0]))
1278
669
      continue;
1279
1.71k
    
if (1.71k
!1.71k
isl_int_is_zero1.71k
(bmap->div[j][1 + total + div]))
1280
107
      return isl_bool_false;
1281
1.71k
  }
1282
923
1283
816
  return isl_bool_true;
1284
923
}
1285
1286
/* Would an expression for div "div" based on inequality "ineq" of "bmap"
1287
 * be a better expression than the current one?
1288
 *
1289
 * If we do not have any expression yet, then any expression would be better.
1290
 * Otherwise we check if the last variable involved in the inequality
1291
 * (disregarding the div that it would define) is in an earlier position
1292
 * than the last variable involved in the current div expression.
1293
 */
1294
static isl_bool better_div_constraint(__isl_keep isl_basic_map *bmap,
1295
  int div, int ineq)
1296
52.2k
{
1297
52.2k
  unsigned total = 1 + isl_space_dim(bmap->dim, isl_dim_all);
1298
52.2k
  int last_div;
1299
52.2k
  int last_ineq;
1300
52.2k
1301
52.2k
  if (isl_int_is_zero(bmap->div[div][0]))
1302
1.40k
    return isl_bool_true;
1303
52.2k
1304
50.8k
  
if (50.8k
isl_seq_last_non_zero(bmap->ineq[ineq] + total + div + 1,50.8k
1305
50.8k
          bmap->n_div - (div + 1)) >= 0)
1306
2.10k
    return isl_bool_false;
1307
50.8k
1308
48.7k
  last_ineq = isl_seq_last_non_zero(bmap->ineq[ineq], total + div);
1309
48.7k
  last_div = isl_seq_last_non_zero(bmap->div[div] + 1,
1310
48.7k
           total + bmap->n_div);
1311
48.7k
1312
48.7k
  return last_ineq < last_div;
1313
50.8k
}
1314
1315
/* Given two constraints "k" and "l" that are opposite to each other,
1316
 * except for the constant term, check if we can use them
1317
 * to obtain an expression for one of the hitherto unknown divs or
1318
 * a "better" expression for a div for which we already have an expression.
1319
 * "sum" is the sum of the constant terms of the constraints.
1320
 * If this sum is strictly smaller than the coefficient of one
1321
 * of the divs, then this pair can be used define the div.
1322
 * To avoid the introduction of circular definitions of divs, we
1323
 * do not use the pair if the resulting expression would refer to
1324
 * any other undefined divs or if any known div is defined in
1325
 * terms of the unknown div.
1326
 */
1327
static struct isl_basic_map *check_for_div_constraints(
1328
  struct isl_basic_map *bmap, int k, int l, isl_int sum, int *progress)
1329
2.86M
{
1330
2.86M
  int i;
1331
2.86M
  unsigned total = 1 + isl_space_dim(bmap->dim, isl_dim_all);
1332
2.86M
1333
2.99M
  for (i = 0; 
i < bmap->n_div2.99M
;
++i129k
)
{181k
1334
181k
    isl_bool set_div;
1335
181k
1336
181k
    if (isl_int_is_zero(bmap->ineq[k][total + i]))
1337
120k
      continue;
1338
60.3k
    
if (60.3k
isl_int_abs_ge60.3k
(sum, bmap->ineq[k][total + i]))
1339
8.14k
      continue;
1340
52.2k
    set_div = better_div_constraint(bmap, i, k);
1341
52.2k
    if (
set_div >= 0 && 52.2k
set_div52.2k
)
1342
1.53k
      set_div = ok_to_set_div_from_bound(bmap, i, k);
1343
52.2k
    if (set_div < 0)
1344
0
      return isl_basic_map_free(bmap);
1345
52.2k
    
if (52.2k
!set_div52.2k
)
1346
51.4k
      break;
1347
813
    
if (813
isl_int_is_pos813
(bmap->ineq[k][total + i]))
1348
272
      bmap = set_div_from_lower_bound(bmap, i, k);
1349
813
    else
1350
541
      bmap = set_div_from_lower_bound(bmap, i, l);
1351
813
    if (progress)
1352
813
      *progress = 1;
1353
813
    break;
1354
52.2k
  }
1355
2.86M
  return bmap;
1356
2.86M
}
1357
1358
__isl_give isl_basic_map *isl_basic_map_remove_duplicate_constraints(
1359
  __isl_take isl_basic_map *bmap, int *progress, int detect_divs)
1360
1.62M
{
1361
1.62M
  struct isl_constraint_index ci;
1362
1.62M
  int k, l, h;
1363
1.62M
  unsigned total = isl_basic_map_total_dim(bmap);
1364
1.62M
  isl_int sum;
1365
1.62M
1366
1.62M
  if (
!bmap || 1.62M
bmap->n_ineq <= 11.62M
)
1367
547k
    return bmap;
1368
1.62M
1369
1.07M
  
if (1.07M
create_constraint_index(&ci, bmap) < 01.07M
)
1370
0
    return bmap;
1371
1.07M
1372
1.07M
  h = isl_seq_get_hash_bits(bmap->ineq[0] + 1, total, ci.bits);
1373
1.07M
  ci.index[h] = &bmap->ineq[0];
1374
6.08M
  for (k = 1; 
k < bmap->n_ineq6.08M
;
++k5.00M
)
{5.00M
1375
5.00M
    h = hash_index(&ci, bmap, k);
1376
5.00M
    if (
!ci.index[h]5.00M
)
{4.60M
1377
4.60M
      ci.index[h] = &bmap->ineq[k];
1378
4.60M
      continue;
1379
4.60M
    }
1380
407k
    
if (407k
progress407k
)
1381
401k
      *progress = 1;
1382
407k
    l = ci.index[h] - &bmap->ineq[0];
1383
407k
    if (isl_int_lt(bmap->ineq[k][0], bmap->ineq[l][0]))
1384
37.8k
      swap_inequality(bmap, k, l);
1385
407k
    isl_basic_map_drop_inequality(bmap, k);
1386
407k
    --k;
1387
407k
  }
1388
1.07M
  isl_int_init(sum);
1389
5.14M
  for (k = 0; 
k < bmap->n_ineq-15.14M
;
++k4.07M
)
{4.17M
1390
4.17M
    isl_seq_neg(bmap->ineq[k]+1, bmap->ineq[k]+1, total);
1391
4.17M
    h = hash_index(&ci, bmap, k);
1392
4.17M
    isl_seq_neg(bmap->ineq[k]+1, bmap->ineq[k]+1, total);
1393
4.17M
    if (!ci.index[h])
1394
1.18M
      continue;
1395
2.98M
    l = ci.index[h] - &bmap->ineq[0];
1396
2.98M
    isl_int_add(sum, bmap->ineq[k][0], bmap->ineq[l][0]);
1397
2.98M
    if (
isl_int_is_pos2.98M
(sum))
{2.88M
1398
2.88M
      if (detect_divs)
1399
2.86M
        bmap = check_for_div_constraints(bmap, k, l,
1400
2.86M
                 sum, progress);
1401
2.88M
      continue;
1402
2.88M
    }
1403
100k
    
if (100k
isl_int_is_zero100k
(sum))
{9.33k
1404
9.33k
      /* We need to break out of the loop after these
1405
9.33k
       * changes since the contents of the hash
1406
9.33k
       * will no longer be valid.
1407
9.33k
       * Plus, we probably we want to regauss first.
1408
9.33k
       */
1409
9.33k
      if (progress)
1410
9.19k
        *progress = 1;
1411
9.33k
      isl_basic_map_drop_inequality(bmap, l);
1412
9.33k
      isl_basic_map_inequality_to_equality(bmap, k);
1413
9.33k
    } else
1414
91.2k
      bmap = isl_basic_map_set_to_empty(bmap);
1415
100k
    break;
1416
2.98M
  }
1417
1.07M
  isl_int_clear(sum);
1418
1.07M
1419
1.07M
  constraint_index_free(&ci);
1420
1.07M
  return bmap;
1421
1.07M
}
1422
1423
/* Detect all pairs of inequalities that form an equality.
1424
 *
1425
 * isl_basic_map_remove_duplicate_constraints detects at most one such pair.
1426
 * Call it repeatedly while it is making progress.
1427
 */
1428
__isl_give isl_basic_map *isl_basic_map_detect_inequality_pairs(
1429
  __isl_take isl_basic_map *bmap, int *progress)
1430
195
{
1431
195
  int duplicate;
1432
195
1433
310
  do {
1434
310
    duplicate = 0;
1435
310
    bmap = isl_basic_map_remove_duplicate_constraints(bmap,
1436
310
                &duplicate, 0);
1437
310
    if (
progress && 310
duplicate191
)
1438
54
      *progress = 1;
1439
310
  } while (duplicate);
1440
195
1441
195
  return bmap;
1442
195
}
1443
1444
/* Eliminate knowns divs from constraints where they appear with
1445
 * a (positive or negative) unit coefficient.
1446
 *
1447
 * That is, replace
1448
 *
1449
 *  floor(e/m) + f >= 0
1450
 *
1451
 * by
1452
 *
1453
 *  e + m f >= 0
1454
 *
1455
 * and
1456
 *
1457
 *  -floor(e/m) + f >= 0
1458
 *
1459
 * by
1460
 *
1461
 *  -e + m f + m - 1 >= 0
1462
 *
1463
 * The first conversion is valid because floor(e/m) >= -f is equivalent
1464
 * to e/m >= -f because -f is an integral expression.
1465
 * The second conversion follows from the fact that
1466
 *
1467
 *  -floor(e/m) = ceil(-e/m) = floor((-e + m - 1)/m)
1468
 *
1469
 *
1470
 * Note that one of the div constraints may have been eliminated
1471
 * due to being redundant with respect to the constraint that is
1472
 * being modified by this function.  The modified constraint may
1473
 * no longer imply this div constraint, so we add it back to make
1474
 * sure we do not lose any information.
1475
 *
1476
 * We skip integral divs, i.e., those with denominator 1, as we would
1477
 * risk eliminating the div from the div constraints.  We do not need
1478
 * to handle those divs here anyway since the div constraints will turn
1479
 * out to form an equality and this equality can then be used to eliminate
1480
 * the div from all constraints.
1481
 */
1482
static __isl_give isl_basic_map *eliminate_unit_divs(
1483
  __isl_take isl_basic_map *bmap, int *progress)
1484
1.60M
{
1485
1.60M
  int i, j;
1486
1.60M
  isl_ctx *ctx;
1487
1.60M
  unsigned total;
1488
1.60M
1489
1.60M
  if (!bmap)
1490
0
    return NULL;
1491
1.60M
1492
1.60M
  ctx = isl_basic_map_get_ctx(bmap);
1493
1.60M
  total = 1 + isl_space_dim(bmap->dim, isl_dim_all);
1494
1.60M
1495
1.83M
  for (i = 0; 
i < bmap->n_div1.83M
;
++i227k
)
{227k
1496
227k
    if (isl_int_is_zero(bmap->div[i][0]))
1497
166k
      continue;
1498
60.3k
    
if (60.3k
isl_int_is_one60.3k
(bmap->div[i][0]))
1499
203
      continue;
1500
473k
    
for (j = 0; 60.1k
j < bmap->n_ineq473k
;
++j413k
)
{413k
1501
413k
      int s;
1502
413k
1503
413k
      if (
!413k
isl_int_is_one413k
(bmap->ineq[j][total + i]) &&
1504
409k
          
!409k
isl_int_is_negone409k
(bmap->ineq[j][total + i]))
1505
406k
        continue;
1506
413k
1507
6.47k
      *progress = 1;
1508
6.47k
1509
6.47k
      s = isl_int_sgn(bmap->ineq[j][total + i]);
1510
6.47k
      isl_int_set_si(bmap->ineq[j][total + i], 0);
1511
6.47k
      if (s < 0)
1512
3.19k
        isl_seq_combine(bmap->ineq[j],
1513
3.19k
          ctx->negone, bmap->div[i] + 1,
1514
3.19k
          bmap->div[i][0], bmap->ineq[j],
1515
3.19k
          total + bmap->n_div);
1516
6.47k
      else
1517
3.27k
        isl_seq_combine(bmap->ineq[j],
1518
3.27k
          ctx->one, bmap->div[i] + 1,
1519
3.27k
          bmap->div[i][0], bmap->ineq[j],
1520
3.27k
          total + bmap->n_div);
1521
6.47k
      if (
s < 06.47k
)
{3.19k
1522
3.19k
        isl_int_add(bmap->ineq[j][0],
1523
3.19k
          bmap->ineq[j][0], bmap->div[i][0]);
1524
3.19k
        isl_int_sub_ui(bmap->ineq[j][0],
1525
3.19k
          bmap->ineq[j][0], 1);
1526
3.19k
      }
1527
6.47k
1528
6.47k
      bmap = isl_basic_map_extend_constraints(bmap, 0, 1);
1529
6.47k
      if (isl_basic_map_add_div_constraint(bmap, i, s) < 0)
1530
0
        return isl_basic_map_free(bmap);
1531
6.47k
    }
1532
60.1k
  }
1533
1.60M
1534
1.60M
  return bmap;
1535
1.60M
}
1536
1537
struct isl_basic_map *isl_basic_map_simplify(struct isl_basic_map *bmap)
1538
1.36M
{
1539
1.36M
  int progress = 1;
1540
1.36M
  if (!bmap)
1541
0
    return NULL;
1542
2.96M
  
while (1.36M
progress2.96M
)
{1.65M
1543
1.65M
    isl_bool empty;
1544
1.65M
1545
1.65M
    progress = 0;
1546
1.65M
    empty = isl_basic_map_plain_is_empty(bmap);
1547
1.65M
    if (empty < 0)
1548
0
      return isl_basic_map_free(bmap);
1549
1.65M
    
if (1.65M
empty1.65M
)
1550
54.3k
      break;
1551
1.60M
    bmap = isl_basic_map_normalize_constraints(bmap);
1552
1.60M
    bmap = reduce_div_coefficients(bmap);
1553
1.60M
    bmap = normalize_div_expressions(bmap);
1554
1.60M
    bmap = remove_duplicate_divs(bmap, &progress);
1555
1.60M
    bmap = eliminate_unit_divs(bmap, &progress);
1556
1.60M
    bmap = eliminate_divs_eq(bmap, &progress);
1557
1.60M
    bmap = eliminate_divs_ineq(bmap, &progress);
1558
1.60M
    bmap = isl_basic_map_gauss(bmap, &progress);
1559
1.60M
    /* requires equalities in normal form */
1560
1.60M
    bmap = normalize_divs(bmap, &progress);
1561
1.60M
    bmap = isl_basic_map_remove_duplicate_constraints(bmap,
1562
1.60M
                &progress, 1);
1563
1.60M
    if (
bmap && 1.60M
progress1.60M
)
1564
1.60M
      ISL_F_CLR(bmap, ISL_BASIC_MAP_REDUCED_COEFFICIENTS);
1565
1.60M
  }
1566
1.36M
  return bmap;
1567
1.36M
}
1568
1569
struct isl_basic_set *isl_basic_set_simplify(struct isl_basic_set *bset)
1570
498k
{
1571
498k
  return bset_from_bmap(isl_basic_map_simplify(bset_to_bmap(bset)));
1572
498k
}
1573
1574
1575
isl_bool isl_basic_map_is_div_constraint(__isl_keep isl_basic_map *bmap,
1576
  isl_int *constraint, unsigned div)
1577
33.3k
{
1578
33.3k
  unsigned pos;
1579
33.3k
1580
33.3k
  if (!bmap)
1581
0
    return isl_bool_error;
1582
33.3k
1583
33.3k
  pos = 1 + isl_space_dim(bmap->dim, isl_dim_all) + div;
1584
33.3k
1585
33.3k
  if (
isl_int_eq33.3k
(constraint[pos], bmap->div[div][0]))
{14.8k
1586
14.8k
    int neg;
1587
14.8k
    isl_int_sub(bmap->div[div][1],
1588
14.8k
        bmap->div[div][1], bmap->div[div][0]);
1589
14.8k
    isl_int_add_ui(bmap->div[div][1], bmap->div[div][1], 1);
1590
14.8k
    neg = isl_seq_is_neg(constraint, bmap->div[div]+1, pos);
1591
14.8k
    isl_int_sub_ui(bmap->div[div][1], bmap->div[div][1], 1);
1592
14.8k
    isl_int_add(bmap->div[div][1],
1593
14.8k
        bmap->div[div][1], bmap->div[div][0]);
1594
14.8k
    if (!neg)
1595
3.19k
      return isl_bool_false;
1596
11.6k
    
if (11.6k
isl_seq_first_non_zero(constraint+pos+1,11.6k
1597
11.6k
              bmap->n_div-div-1) != -1)
1598
0
      return isl_bool_false;
1599
18.5k
  } else 
if (18.5k
isl_int_abs_eq18.5k
(constraint[pos], bmap->div[div][0]))
{15.3k
1600
15.3k
    if (!isl_seq_eq(constraint, bmap->div[div]+1, pos))
1601
5.16k
      return isl_bool_false;
1602
10.2k
    
if (10.2k
isl_seq_first_non_zero(constraint+pos+1,10.2k
1603
10.2k
              bmap->n_div-div-1) != -1)
1604
1
      return isl_bool_false;
1605
10.2k
  } else
1606
3.11k
    return isl_bool_false;
1607
33.3k
1608
21.8k
  return isl_bool_true;
1609
33.3k
}
1610
1611
isl_bool isl_basic_set_is_div_constraint(__isl_keep isl_basic_set *bset,
1612
  isl_int *constraint, unsigned div)
1613
0
{
1614
0
  return isl_basic_map_is_div_constraint(bset, constraint, div);
1615
0
}
1616
1617
1618
/* If the only constraints a div d=floor(f/m)
1619
 * appears in are its two defining constraints
1620
 *
1621
 *  f - m d >=0
1622
 *  -(f - (m - 1)) + m d >= 0
1623
 *
1624
 * then it can safely be removed.
1625
 */
1626
static isl_bool div_is_redundant(struct isl_basic_map *bmap, int div)
1627
66.1k
{
1628
66.1k
  int i;
1629
66.1k
  unsigned pos = 1 + isl_space_dim(bmap->dim, isl_dim_all) + div;
1630
66.1k
1631
88.2k
  for (i = 0; 
i < bmap->n_eq88.2k
;
++i22.1k
)
1632
53.4k
    
if (53.4k
!53.4k
isl_int_is_zero53.4k
(bmap->eq[i][pos]))
1633
31.2k
      return isl_bool_false;
1634
66.1k
1635
119k
  
for (i = 0; 34.8k
i < bmap->n_ineq119k
;
++i84.2k
)
{95.6k
1636
95.6k
    isl_bool red;
1637
95.6k
1638
95.6k
    if (isl_int_is_zero(bmap->ineq[i][pos]))
1639
62.3k
      continue;
1640
33.2k
    red = isl_basic_map_is_div_constraint(bmap, bmap->ineq[i], div);
1641
33.2k
    if (
red < 0 || 33.2k
!red33.2k
)
1642
11.4k
      return red;
1643
33.2k
  }
1644
34.8k
1645
50.2k
  
for (i = 0; 23.4k
i < bmap->n_div50.2k
;
++i26.8k
)
{26.8k
1646
26.8k
    if (isl_int_is_zero(bmap->div[i][0]))
1647
120
      continue;
1648
26.7k
    
if (26.7k
!26.7k
isl_int_is_zero26.7k
(bmap->div[i][1+pos]))
1649
1
      return isl_bool_false;
1650
26.7k
  }
1651
23.4k
1652
23.4k
  return isl_bool_true;
1653
23.4k
}
1654
1655
/*
1656
 * Remove divs that don't occur in any of the constraints or other divs.
1657
 * These can arise when dropping constraints from a basic map or
1658
 * when the divs of a basic map have been temporarily aligned
1659
 * with the divs of another basic map.
1660
 */
1661
static struct isl_basic_map *remove_redundant_divs(struct isl_basic_map *bmap)
1662
3.09M
{
1663
3.09M
  int i;
1664
3.09M
1665
3.09M
  if (!bmap)
1666
0
    return NULL;
1667
3.09M
1668
3.16M
  
for (i = bmap->n_div-1; 3.09M
i >= 03.16M
;
--i66.1k
)
{66.1k
1669
66.1k
    isl_bool redundant;
1670
66.1k
1671
66.1k
    redundant = div_is_redundant(bmap, i);
1672
66.1k
    if (redundant < 0)
1673
0
      return isl_basic_map_free(bmap);
1674
66.1k
    
if (66.1k
!redundant66.1k
)
1675
42.6k
      continue;
1676
23.4k
    bmap = isl_basic_map_drop_div(bmap, i);
1677
23.4k
  }
1678
3.09M
  return bmap;
1679
3.09M
}
1680
1681
/* Mark "bmap" as final, without checking for obviously redundant
1682
 * integer divisions.  This function should be used when "bmap"
1683
 * is known not to involve any such integer divisions.
1684
 */
1685
__isl_give isl_basic_map *isl_basic_map_mark_final(
1686
  __isl_take isl_basic_map *bmap)
1687
3.10M
{
1688
3.10M
  if (!bmap)
1689
0
    return NULL;
1690
3.10M
  
ISL_F_SET3.10M
(bmap, ISL_BASIC_SET_FINAL);3.10M
1691
3.10M
  return bmap;
1692
3.10M
}
1693
1694
/* Mark "bmap" as final, after removing obviously redundant integer divisions.
1695
 */
1696
struct isl_basic_map *isl_basic_map_finalize(struct isl_basic_map *bmap)
1697
3.09M
{
1698
3.09M
  bmap = remove_redundant_divs(bmap);
1699
3.09M
  bmap = isl_basic_map_mark_final(bmap);
1700
3.09M
  return bmap;
1701
3.09M
}
1702
1703
struct isl_basic_set *isl_basic_set_finalize(struct isl_basic_set *bset)
1704
589k
{
1705
589k
  return bset_from_bmap(isl_basic_map_finalize(bset_to_bmap(bset)));
1706
589k
}
1707
1708
/* Remove definition of any div that is defined in terms of the given variable.
1709
 * The div itself is not removed.  Functions such as
1710
 * eliminate_divs_ineq depend on the other divs remaining in place.
1711
 */
1712
static struct isl_basic_map *remove_dependent_vars(struct isl_basic_map *bmap,
1713
                  int pos)
1714
62.9k
{
1715
62.9k
  int i;
1716
62.9k
1717
62.9k
  if (!bmap)
1718
0
    return NULL;
1719
62.9k
1720
131k
  
for (i = 0; 62.9k
i < bmap->n_div131k
;
++i68.6k
)
{68.6k
1721
68.6k
    if (isl_int_is_zero(bmap->div[i][0]))
1722
67.1k
      continue;
1723
1.44k
    
if (1.44k
isl_int_is_zero1.44k
(bmap->div[i][1+1+pos]))
1724
1.38k
      continue;
1725
65
    bmap = isl_basic_map_mark_div_unknown(bmap, i);
1726
65
    if (!bmap)
1727
0
      return NULL;
1728
65
  }
1729
62.9k
  return bmap;
1730
62.9k
}
1731
1732
/* Eliminate the specified variables from the constraints using
1733
 * Fourier-Motzkin.  The variables themselves are not removed.
1734
 */
1735
struct isl_basic_map *isl_basic_map_eliminate_vars(
1736
  struct isl_basic_map *bmap, unsigned pos, unsigned n)
1737
63.3k
{
1738
63.3k
  int d;
1739
63.3k
  int i, j, k;
1740
63.3k
  unsigned total;
1741
63.3k
  int need_gauss = 0;
1742
63.3k
1743
63.3k
  if (n == 0)
1744
1.57k
    return bmap;
1745
61.7k
  
if (61.7k
!bmap61.7k
)
1746
0
    return NULL;
1747
61.7k
  total = isl_basic_map_total_dim(bmap);
1748
61.7k
1749
61.7k
  bmap = isl_basic_map_cow(bmap);
1750
124k
  for (d = pos + n - 1; 
d >= 0 && 124k
d >= pos112k
;
--d62.9k
)
1751
62.9k
    bmap = remove_dependent_vars(bmap, d);
1752
61.7k
  if (!bmap)
1753
0
    return NULL;
1754
61.7k
1755
61.7k
  for (d = pos + n - 1;
1756
108k
       
d >= 0 && 108k
d >= total - bmap->n_div100k
&&
d >= pos59.1k
;
--d46.9k
)
1757
46.9k
    isl_seq_clr(bmap->div[d-(total-bmap->n_div)], 2+total);
1758
124k
  for (d = pos + n - 1; 
d >= 0 && 124k
d >= pos112k
;
--d62.9k
)
{62.9k
1759
62.9k
    int n_lower, n_upper;
1760
62.9k
    if (!bmap)
1761
0
      return NULL;
1762
78.7k
    
for (i = 0; 62.9k
i < bmap->n_eq78.7k
;
++i15.7k
)
{17.7k
1763
17.7k
      if (isl_int_is_zero(bmap->eq[i][1+d]))
1764
15.7k
        continue;
1765
1.93k
      eliminate_var_using_equality(bmap, d, bmap->eq[i], 0, NULL);
1766
1.93k
      isl_basic_map_drop_equality(bmap, i);
1767
1.93k
      need_gauss = 1;
1768
1.93k
      break;
1769
17.7k
    }
1770
62.9k
    if (i < bmap->n_eq)
1771
726
      continue;
1772
62.2k
    n_lower = 0;
1773
62.2k
    n_upper = 0;
1774
203k
    for (i = 0; 
i < bmap->n_ineq203k
;
++i141k
)
{141k
1775
141k
      if (isl_int_is_pos(bmap->ineq[i][1+d]))
1776
24.3k
        n_lower++;
1777
116k
      else 
if (116k
isl_int_is_neg116k
(bmap->ineq[i][1+d]))
1778
23.0k
        n_upper++;
1779
141k
    }
1780
62.2k
    bmap = isl_basic_map_extend_constraints(bmap,
1781
62.2k
        0, n_lower * n_upper);
1782
62.2k
    if (!bmap)
1783
0
      goto error;
1784
161k
    
for (i = bmap->n_ineq - 1; 62.2k
i >= 0161k
;
--i98.8k
)
{98.8k
1785
98.8k
      int last;
1786
98.8k
      if (isl_int_is_zero(bmap->ineq[i][1+d]))
1787
51.5k
        continue;
1788
47.3k
      last = -1;
1789
186k
      for (j = 0; 
j < i186k
;
++j139k
)
{139k
1790
139k
        if (isl_int_is_zero(bmap->ineq[j][1+d]))
1791
99.4k
          continue;
1792
40.1k
        last = j;
1793
40.1k
        if (
isl_int_sgn40.1k
(bmap->ineq[i][1+d]) ==40.1k
1794
40.1k
            isl_int_sgn(bmap->ineq[j][1+d]))
1795
9.37k
          continue;
1796
30.7k
        k = isl_basic_map_alloc_inequality(bmap);
1797
30.7k
        if (k < 0)
1798
0
          goto error;
1799
30.7k
        isl_seq_cpy(bmap->ineq[k], bmap->ineq[i],
1800
30.7k
            1+total);
1801
30.7k
        isl_seq_elim(bmap->ineq[k], bmap->ineq[j],
1802
30.7k
            1+d, 1+total, NULL);
1803
30.7k
      }
1804
47.3k
      isl_basic_map_drop_inequality(bmap, i);
1805
47.3k
      i = last + 1;
1806
47.3k
    }
1807
62.2k
    
if (62.2k
n_lower > 0 && 62.2k
n_upper > 020.5k
)
{19.5k
1808
19.5k
      bmap = isl_basic_map_normalize_constraints(bmap);
1809
19.5k
      bmap = isl_basic_map_remove_duplicate_constraints(bmap,
1810
19.5k
                    NULL, 0);
1811
19.5k
      bmap = isl_basic_map_gauss(bmap, NULL);
1812
19.5k
      bmap = isl_basic_map_remove_redundancies(bmap);
1813
19.5k
      need_gauss = 0;
1814
19.5k
      if (!bmap)
1815
0
        goto error;
1816
19.5k
      
if (19.5k
ISL_F_ISSET19.5k
(bmap, ISL_BASIC_MAP_EMPTY))
1817
67
        break;
1818
19.5k
    }
1819
62.2k
  }
1820
61.7k
  
ISL_F_CLR61.7k
(bmap, ISL_BASIC_MAP_NORMALIZED);61.7k
1821
61.7k
  if (need_gauss)
1822
927
    bmap = isl_basic_map_gauss(bmap, NULL);
1823
61.7k
  return bmap;
1824
0
error:
1825
0
  isl_basic_map_free(bmap);
1826
0
  return NULL;
1827
61.7k
}
1828
1829
struct isl_basic_set *isl_basic_set_eliminate_vars(
1830
  struct isl_basic_set *bset, unsigned pos, unsigned n)
1831
14.4k
{
1832
14.4k
  return bset_from_bmap(isl_basic_map_eliminate_vars(bset_to_bmap(bset),
1833
14.4k
                pos, n));
1834
14.4k
}
1835
1836
/* Eliminate the specified n dimensions starting at first from the
1837
 * constraints, without removing the dimensions from the space.
1838
 * If the set is rational, the dimensions are eliminated using Fourier-Motzkin.
1839
 * Otherwise, they are projected out and the original space is restored.
1840
 */
1841
__isl_give isl_basic_map *isl_basic_map_eliminate(
1842
  __isl_take isl_basic_map *bmap,
1843
  enum isl_dim_type type, unsigned first, unsigned n)
1844
8.21k
{
1845
8.21k
  isl_space *space;
1846
8.21k
1847
8.21k
  if (!bmap)
1848
0
    return NULL;
1849
8.21k
  
if (8.21k
n == 08.21k
)
1850
0
    return bmap;
1851
8.21k
1852
8.21k
  
if (8.21k
first + n > isl_basic_map_dim(bmap, type) || 8.21k
first + n < first8.21k
)
1853
0
    isl_die(bmap->ctx, isl_error_invalid,
1854
8.21k
      "index out of bounds", goto error);
1855
8.21k
1856
8.21k
  
if (8.21k
ISL_F_ISSET8.21k
(bmap, ISL_BASIC_MAP_RATIONAL))
{0
1857
0
    first += isl_basic_map_offset(bmap, type) - 1;
1858
0
    bmap = isl_basic_map_eliminate_vars(bmap, first, n);
1859
0
    return isl_basic_map_finalize(bmap);
1860
0
  }
1861
8.21k
1862
8.21k
  space = isl_basic_map_get_space(bmap);
1863
8.21k
  bmap = isl_basic_map_project_out(bmap, type, first, n);
1864
8.21k
  bmap = isl_basic_map_insert_dims(bmap, type, first, n);
1865
8.21k
  bmap = isl_basic_map_reset_space(bmap, space);
1866
8.21k
  return bmap;
1867
0
error:
1868
0
  isl_basic_map_free(bmap);
1869
0
  return NULL;
1870
8.21k
}
1871
1872
__isl_give isl_basic_set *isl_basic_set_eliminate(
1873
  __isl_take isl_basic_set *bset,
1874
  enum isl_dim_type type, unsigned first, unsigned n)
1875
2.56k
{
1876
2.56k
  return isl_basic_map_eliminate(bset, type, first, n);
1877
2.56k
}
1878
1879
/* Remove all constraints from "bmap" that reference any unknown local
1880
 * variables (directly or indirectly).
1881
 *
1882
 * Dropping all constraints on a local variable will make it redundant,
1883
 * so it will get removed implicitly by
1884
 * isl_basic_map_drop_constraints_involving_dims.  Some other local
1885
 * variables may also end up becoming redundant if they only appear
1886
 * in constraints together with the unknown local variable.
1887
 * Therefore, start over after calling
1888
 * isl_basic_map_drop_constraints_involving_dims.
1889
 */
1890
__isl_give isl_basic_map *isl_basic_map_drop_constraint_involving_unknown_divs(
1891
  __isl_take isl_basic_map *bmap)
1892
5.06k
{
1893
5.06k
  isl_bool known;
1894
5.06k
  int i, n_div, o_div;
1895
5.06k
1896
5.06k
  known = isl_basic_map_divs_known(bmap);
1897
5.06k
  if (known < 0)
1898
0
    return isl_basic_map_free(bmap);
1899
5.06k
  
if (5.06k
known5.06k
)
1900
5.06k
    return bmap;
1901
5.06k
1902
0
  n_div = isl_basic_map_dim(bmap, isl_dim_div);
1903
0
  o_div = isl_basic_map_offset(bmap, isl_dim_div) - 1;
1904
0
1905
0
  for (i = 0; 
i < n_div0
;
++i0
)
{0
1906
0
    known = isl_basic_map_div_is_known(bmap, i);
1907
0
    if (known < 0)
1908
0
      return isl_basic_map_free(bmap);
1909
0
    
if (0
known0
)
1910
0
      continue;
1911
0
    bmap = remove_dependent_vars(bmap, o_div + i);
1912
0
    bmap = isl_basic_map_drop_constraints_involving_dims(bmap,
1913
0
                  isl_dim_div, i, 1);
1914
0
    if (!bmap)
1915
0
      return NULL;
1916
0
    n_div = isl_basic_map_dim(bmap, isl_dim_div);
1917
0
    i = -1;
1918
0
  }
1919
0
1920
0
  return bmap;
1921
0
}
1922
1923
/* Remove all constraints from "map" that reference any unknown local
1924
 * variables (directly or indirectly).
1925
 *
1926
 * Since constraints may get dropped from the basic maps,
1927
 * they may no longer be disjoint from each other.
1928
 */
1929
__isl_give isl_map *isl_map_drop_constraint_involving_unknown_divs(
1930
  __isl_take isl_map *map)
1931
933
{
1932
933
  int i;
1933
933
  isl_bool known;
1934
933
1935
933
  known = isl_map_divs_known(map);
1936
933
  if (known < 0)
1937
0
    return isl_map_free(map);
1938
933
  
if (933
known933
)
1939
933
    return map;
1940
933
1941
0
  map = isl_map_cow(map);
1942
0
  if (!map)
1943
0
    return NULL;
1944
0
1945
0
  
for (i = 0; 0
i < map->n0
;
++i0
)
{0
1946
0
    map->p[i] =
1947
0
        isl_basic_map_drop_constraint_involving_unknown_divs(
1948
0
                    map->p[i]);
1949
0
    if (!map->p[i])
1950
0
      return isl_map_free(map);
1951
0
  }
1952
0
1953
0
  
if (0
map->n > 10
)
1954
0
    ISL_F_CLR(map, ISL_MAP_DISJOINT);
1955
0
1956
0
  return map;
1957
0
}
1958
1959
/* Don't assume equalities are in order, because align_divs
1960
 * may have changed the order of the divs.
1961
 */
1962
static void compute_elimination_index(struct isl_basic_map *bmap, int *elim)
1963
13.4k
{
1964
13.4k
  int d, i;
1965
13.4k
  unsigned total;
1966
13.4k
1967
13.4k
  total = isl_space_dim(bmap->dim, isl_dim_all);
1968
113k
  for (d = 0; 
d < total113k
;
++d100k
)
1969
100k
    elim[d] = -1;
1970
28.6k
  for (i = 0; 
i < bmap->n_eq28.6k
;
++i15.1k
)
{15.1k
1971
52.4k
    for (d = total - 1; 
d >= 052.4k
;
--d37.2k
)
{52.4k
1972
52.4k
      if (isl_int_is_zero(bmap->eq[i][1+d]))
1973
37.2k
        continue;
1974
15.1k
      elim[d] = i;
1975
15.1k
      break;
1976
52.4k
    }
1977
15.1k
  }
1978
13.4k
}
1979
1980
static void set_compute_elimination_index(struct isl_basic_set *bset, int *elim)
1981
242
{
1982
242
  compute_elimination_index(bset_to_bmap(bset), elim);
1983
242
}
1984
1985
static int reduced_using_equalities(isl_int *dst, isl_int *src,
1986
  struct isl_basic_map *bmap, int *elim)
1987
116k
{
1988
116k
  int d;
1989
116k
  int copied = 0;
1990
116k
  unsigned total;
1991
116k
1992
116k
  total = isl_space_dim(bmap->dim, isl_dim_all);
1993
1.11M
  for (d = total - 1; 
d >= 01.11M
;
--d1.00M
)
{1.00M
1994
1.00M
    if (isl_int_is_zero(src[1+d]))
1995
872k
      continue;
1996
129k
    
if (129k
elim[d] == -1129k
)
1997
119k
      continue;
1998
9.96k
    
if (9.96k
!copied9.96k
)
{9.70k
1999
9.70k
      isl_seq_cpy(dst, src, 1 + total);
2000
9.70k
      copied = 1;
2001
9.70k
    }
2002
9.96k
    isl_seq_elim(dst, bmap->eq[elim[d]], 1 + d, 1 + total, NULL);
2003
9.96k
  }
2004
116k
  return copied;
2005
116k
}
2006
2007
static int set_reduced_using_equalities(isl_int *dst, isl_int *src,
2008
  struct isl_basic_set *bset, int *elim)
2009
525
{
2010
525
  return reduced_using_equalities(dst, src,
2011
525
          bset_to_bmap(bset), elim);
2012
525
}
2013
2014
static struct isl_basic_set *isl_basic_set_reduce_using_equalities(
2015
  struct isl_basic_set *bset, struct isl_basic_set *context)
2016
6.42k
{
2017
6.42k
  int i;
2018
6.42k
  int *elim;
2019
6.42k
2020
6.42k
  if (
!bset || 6.42k
!context6.42k
)
2021
0
    goto error;
2022
6.42k
2023
6.42k
  
if (6.42k
context->n_eq == 06.42k
)
{6.18k
2024
6.18k
    isl_basic_set_free(context);
2025
6.18k
    return bset;
2026
6.18k
  }
2027
6.42k
2028
242
  bset = isl_basic_set_cow(bset);
2029
242
  if (!bset)
2030
0
    goto error;
2031
242
2032
242
  
elim = 242
isl_alloc_array242
(bset->ctx, int, isl_basic_set_n_dim(bset));
2033
242
  if (!elim)
2034
0
    goto error;
2035
242
  set_compute_elimination_index(context, elim);
2036
698
  for (i = 0; 
i < bset->n_eq698
;
++i456
)
2037
456
    set_reduced_using_equalities(bset->eq[i], bset->eq[i],
2038
456
              context, elim);
2039
311
  for (i = 0; 
i < bset->n_ineq311
;
++i69
)
2040
69
    set_reduced_using_equalities(bset->ineq[i], bset->ineq[i],
2041
69
              context, elim);
2042
242
  isl_basic_set_free(context);
2043
242
  free(elim);
2044
242
  bset = isl_basic_set_simplify(bset);
2045
242
  bset = isl_basic_set_finalize(bset);
2046
242
  return bset;
2047
0
error:
2048
0
  isl_basic_set_free(bset);
2049
0
  isl_basic_set_free(context);
2050
0
  return NULL;
2051
242
}
2052
2053
/* For each inequality in "ineq" that is a shifted (more relaxed)
2054
 * copy of an inequality in "context", mark the corresponding entry
2055
 * in "row" with -1.
2056
 * If an inequality only has a non-negative constant term, then
2057
 * mark it as well.
2058
 */
2059
static isl_stat mark_shifted_constraints(__isl_keep isl_mat *ineq,
2060
  __isl_keep isl_basic_set *context, int *row)
2061
8.91k
{
2062
8.91k
  struct isl_constraint_index ci;
2063
8.91k
  int n_ineq;
2064
8.91k
  unsigned total;
2065
8.91k
  int k;
2066
8.91k
2067
8.91k
  if (
!ineq || 8.91k
!context8.91k
)
2068
0
    return isl_stat_error;
2069
8.91k
  
if (8.91k
context->n_ineq == 08.91k
)
2070
0
    return isl_stat_ok;
2071
8.91k
  
if (8.91k
setup_constraint_index(&ci, context) < 08.91k
)
2072
0
    return isl_stat_error;
2073
8.91k
2074
8.91k
  n_ineq = isl_mat_rows(ineq);
2075
8.91k
  total = isl_mat_cols(ineq) - 1;
2076
32.1k
  for (k = 0; 
k < n_ineq32.1k
;
++k23.2k
)
{23.2k
2077
23.2k
    int l;
2078
23.2k
    isl_bool redundant;
2079
23.2k
2080
23.2k
    l = isl_seq_first_non_zero(ineq->row[k] + 1, total);
2081
23.2k
    if (
l < 0 && 23.2k
isl_int_is_nonneg4
(ineq->row[k][0]))
{0
2082
0
      row[k] = -1;
2083
0
      continue;
2084
0
    }
2085
23.2k
    redundant = constraint_index_is_redundant(&ci, ineq->row[k]);
2086
23.2k
    if (redundant < 0)
2087
0
      goto error;
2088
23.2k
    
if (23.2k
!redundant23.2k
)
2089
16.3k
      continue;
2090
6.95k
    row[k] = -1;
2091
6.95k
  }
2092
8.91k
  constraint_index_free(&ci);
2093
8.91k
  return isl_stat_ok;
2094
0
error:
2095
0
  constraint_index_free(&ci);
2096
0
  return isl_stat_error;
2097
8.91k
}
2098
2099
static struct isl_basic_set *remove_shifted_constraints(
2100
  struct isl_basic_set *bset, struct isl_basic_set *context)
2101
33
{
2102
33
  struct isl_constraint_index ci;
2103
33
  int k;
2104
33
2105
33
  if (
!bset || 33
!context33
)
2106
0
    return bset;
2107
33
2108
33
  
if (33
context->n_ineq == 033
)
2109
0
    return bset;
2110
33
  
if (33
setup_constraint_index(&ci, context) < 033
)
2111
0
    return bset;
2112
33
2113
271
  
for (k = 0; 33
k < bset->n_ineq271
;
++k238
)
{238
2114
238
    isl_bool redundant;
2115
238
2116
238
    redundant = constraint_index_is_redundant(&ci, bset->ineq[k]);
2117
238
    if (redundant < 0)
2118
0
      goto error;
2119
238
    
if (238
!redundant238
)
2120
81
      continue;
2121
157
    bset = isl_basic_set_cow(bset);
2122
157
    if (!bset)
2123
0
      goto error;
2124
157
    isl_basic_set_drop_inequality(bset, k);
2125
157
    --k;
2126
157
  }
2127
33
  constraint_index_free(&ci);
2128
33
  return bset;
2129
0
error:
2130
0
  constraint_index_free(&ci);
2131
0
  return bset;
2132
33
}
2133
2134
/* Remove constraints from "bmap" that are identical to constraints
2135
 * in "context" or that are more relaxed (greater constant term).
2136
 *
2137
 * We perform the test for shifted copies on the pure constraints
2138
 * in remove_shifted_constraints.
2139
 */
2140
static __isl_give isl_basic_map *isl_basic_map_remove_shifted_constraints(
2141
  __isl_take isl_basic_map *bmap, __isl_take isl_basic_map *context)
2142
403
{
2143
403
  isl_basic_set *bset, *bset_context;
2144
403
2145
403
  if (
!bmap || 403
!context403
)
2146
0
    goto error;
2147
403
2148
403
  
if (403
bmap->n_ineq == 0 || 403
context->n_ineq == 033
)
{370
2149
370
    isl_basic_map_free(context);
2150
370
    return bmap;
2151
370
  }
2152
403
2153
33
  context = isl_basic_map_align_divs(context, bmap);
2154
33
  bmap = isl_basic_map_align_divs(bmap, context);
2155
33
2156
33
  bset = isl_basic_map_underlying_set(isl_basic_map_copy(bmap));
2157
33
  bset_context = isl_basic_map_underlying_set(context);
2158
33
  bset = remove_shifted_constraints(bset, bset_context);
2159
33
  isl_basic_set_free(bset_context);
2160
33
2161
33
  bmap = isl_basic_map_overlying_set(bset, bmap);
2162
33
2163
33
  return bmap;
2164
0
error:
2165
0
  isl_basic_map_free(bmap);
2166
0
  isl_basic_map_free(context);
2167
0
  return NULL;
2168
403
}
2169
2170
/* Does the (linear part of a) constraint "c" involve any of the "len"
2171
 * "relevant" dimensions?
2172
 */
2173
static int is_related(isl_int *c, int len, int *relevant)
2174
134k
{
2175
134k
  int i;
2176
134k
2177
1.85M
  for (i = 0; 
i < len1.85M
;
++i1.72M
)
{1.77M
2178
1.77M
    if (!relevant[i])
2179
1.38M
      continue;
2180
395k
    
if (395k
!395k
isl_int_is_zero395k
(c[i]))
2181
53.1k
      return 1;
2182
395k
  }
2183
134k
2184
81.7k
  return 0;
2185
134k
}
2186
2187
/* Drop constraints from "bmap" that do not involve any of
2188
 * the dimensions marked "relevant".
2189
 */
2190
static __isl_give isl_basic_map *drop_unrelated_constraints(
2191
  __isl_take isl_basic_map *bmap, int *relevant)
2192
40.8k
{
2193
40.8k
  int i, dim;
2194
40.8k
2195
40.8k
  dim = isl_basic_map_dim(bmap, isl_dim_all);
2196
175k
  for (i = 0; 
i < dim175k
;
++i134k
)
2197
158k
    
if (158k
!relevant[i]158k
)
2198
23.8k
      break;
2199
40.8k
  if (i >= dim)
2200
16.9k
    return bmap;
2201
40.8k
2202
35.6k
  
for (i = bmap->n_eq - 1; 23.8k
i >= 035.6k
;
--i11.7k
)
2203
11.7k
    
if (11.7k
!is_related(bmap->eq[i] + 1, dim, relevant)11.7k
)
{6.29k
2204
6.29k
      bmap = isl_basic_map_cow(bmap);
2205
6.29k
      if (isl_basic_map_drop_equality(bmap, i) < 0)
2206
0
        return isl_basic_map_free(bmap);
2207
6.29k
    }
2208
23.8k
2209
147k
  
for (i = bmap->n_ineq - 1; 23.8k
i >= 0147k
;
--i123k
)
2210
123k
    
if (123k
!is_related(bmap->ineq[i] + 1, dim, relevant)123k
)
{75.4k
2211
75.4k
      bmap = isl_basic_map_cow(bmap);
2212
75.4k
      if (isl_basic_map_drop_inequality(bmap, i) < 0)
2213
0
        return isl_basic_map_free(bmap);
2214
75.4k
    }
2215
23.8k
2216
23.8k
  return bmap;
2217
23.8k
}
2218
2219
/* Update the groups in "group" based on the (linear part of a) constraint "c".
2220
 *
2221
 * In particular, for any variable involved in the constraint,
2222
 * find the actual group id from before and replace the group
2223
 * of the corresponding variable by the minimal group of all
2224
 * the variables involved in the constraint considered so far
2225
 * (if this minimum is smaller) or replace the minimum by this group
2226
 * (if the minimum is larger).
2227
 *
2228
 * At the end, all the variables in "c" will (indirectly) point
2229
 * to the minimal of the groups that they referred to originally.
2230
 */
2231
static void update_groups(int dim, int *group, isl_int *c)
2232
276k
{
2233
276k
  int j;
2234
276k
  int min = dim;
2235
276k
2236
3.75M
  for (j = 0; 
j < dim3.75M
;
++j3.47M
)
{3.47M
2237
3.47M
    if (isl_int_is_zero(c[j]))
2238
3.11M
      continue;
2239
364k
    
while (364k
group[j] >= 0 && 364k
group[group[j]] != group[j]131k
)
2240
48
      group[j] = group[group[j]];
2241
364k
    if (group[j] == min)
2242
49.8k
      continue;
2243
314k
    
if (314k
group[j] < min314k
)
{279k
2244
279k
      if (
min >= 0 && 279k
min < dim279k
)
2245
3.37k
        group[min] = group[j];
2246
279k
      min = group[j];
2247
279k
    } else
2248
35.1k
      group[group[j]] = min;
2249
314k
  }
2250
276k
}
2251
2252
/* Allocate an array of groups of variables, one for each variable
2253
 * in "context", initialized to zero.
2254
 */
2255
static int *alloc_groups(__isl_keep isl_basic_set *context)
2256
21.3k
{
2257
21.3k
  isl_ctx *ctx;
2258
21.3k
  int dim;
2259
21.3k
2260
21.3k
  dim = isl_basic_set_dim(context, isl_dim_set);
2261
21.3k
  ctx = isl_basic_set_get_ctx(context);
2262
21.3k
  return isl_calloc_array(ctx, int, dim);
2263
21.3k
}
2264
2265
/* Drop constraints from "bmap" that only involve variables that are
2266
 * not related to any of the variables marked with a "-1" in "group".
2267
 *
2268
 * We construct groups of variables that collect variables that
2269
 * (indirectly) appear in some common constraint of "bmap".
2270
 * Each group is identified by the first variable in the group,
2271
 * except for the special group of variables that was already identified
2272
 * in the input as -1 (or are related to those variables).
2273
 * If group[i] is equal to i (or -1), then the group of i is i (or -1),
2274
 * otherwise the group of i is the group of group[i].
2275
 *
2276
 * We first initialize groups for the remaining variables.
2277
 * Then we iterate over the constraints of "bmap" and update the
2278
 * group of the variables in the constraint by the smallest group.
2279
 * Finally, we resolve indirect references to groups by running over
2280
 * the variables.
2281
 *
2282
 * After computing the groups, we drop constraints that do not involve
2283
 * any variables in the -1 group.
2284
 */
2285
__isl_give isl_basic_map *isl_basic_map_drop_unrelated_constraints(
2286
  __isl_take isl_basic_map *bmap, __isl_take int *group)
2287
51.1k
{
2288
51.1k
  int dim;
2289
51.1k
  int i;
2290
51.1k
  int last;
2291
51.1k
2292
51.1k
  if (!bmap)
2293
0
    return NULL;
2294
51.1k
2295
51.1k
  dim = isl_basic_map_dim(bmap, isl_dim_all);
2296
51.1k
2297
51.1k
  last = -1;
2298
304k
  for (i = 0; 
i < dim304k
;
++i253k
)
2299
253k
    
if (253k
group[i] >= 0253k
)
2300
109k
      last = group[i] = i;
2301
51.1k
  if (
last < 051.1k
)
{10.2k
2302
10.2k
    free(group);
2303
10.2k
    return bmap;
2304
10.2k
  }
2305
51.1k
2306
80.9k
  
for (i = 0; 40.8k
i < bmap->n_eq80.9k
;
++i40.1k
)
2307
40.1k
    update_groups(dim, group, bmap->eq[i] + 1);
2308
276k
  for (i = 0; 
i < bmap->n_ineq276k
;
++i236k
)
2309
236k
    update_groups(dim, group, bmap->ineq[i] + 1);
2310
40.8k
2311
269k
  for (i = 0; 
i < dim269k
;
++i228k
)
2312
228k
    
if (228k
group[i] >= 0228k
)
2313
73.2k
      group[i] = group[group[i]];
2314
40.8k
2315
269k
  for (i = 0; 
i < dim269k
;
++i228k
)
2316
228k
    group[i] = group[i] == -1;
2317
40.8k
2318
40.8k
  bmap = drop_unrelated_constraints(bmap, group);
2319
40.8k
2320
40.8k
  free(group);
2321
40.8k
  return bmap;
2322
51.1k
}
2323
2324
/* Drop constraints from "context" that are irrelevant for computing
2325
 * the gist of "bset".
2326
 *
2327
 * In particular, drop constraints in variables that are not related
2328
 * to any of the variables involved in the constraints of "bset"
2329
 * in the sense that there is no sequence of constraints that connects them.
2330
 *
2331
 * We first mark all variables that appear in "bset" as belonging
2332
 * to a "-1" group and then continue with group_and_drop_irrelevant_constraints.
2333
 */
2334
static __isl_give isl_basic_set *drop_irrelevant_constraints(
2335
  __isl_take isl_basic_set *context, __isl_keep isl_basic_set *bset)
2336
13.4k
{
2337
13.4k
  int *group;
2338
13.4k
  int dim;
2339
13.4k
  int i, j;
2340
13.4k
2341
13.4k
  if (
!context || 13.4k
!bset13.4k
)
2342
0
    return isl_basic_set_free(context);
2343
13.4k
2344
13.4k
  group = alloc_groups(context);
2345
13.4k
2346
13.4k
  if (!group)
2347
0
    return isl_basic_set_free(context);
2348
13.4k
2349
13.4k
  dim = isl_basic_set_dim(bset, isl_dim_set);
2350
82.7k
  for (i = 0; 
i < dim82.7k
;
++i69.3k
)
{69.3k
2351
115k
    for (j = 0; 
j < bset->n_eq115k
;
++j46.1k
)
2352
63.6k
      
if (63.6k
!63.6k
isl_int_is_zero63.6k
(bset->eq[j][1 + i]))
2353
17.4k
        break;
2354
69.3k
    if (
j < bset->n_eq69.3k
)
{17.4k
2355
17.4k
      group[i] = -1;
2356
17.4k
      continue;
2357
17.4k
    }
2358
172k
    
for (j = 0; 51.8k
j < bset->n_ineq172k
;
++j120k
)
2359
138k
      
if (138k
!138k
isl_int_is_zero138k
(bset->ineq[j][1 + i]))
2360
17.9k
        break;
2361
51.8k
    if (j < bset->n_ineq)
2362
17.9k
      group[i] = -1;
2363
51.8k
  }
2364
13.4k
2365
13.4k
  return isl_basic_map_drop_unrelated_constraints(context, group);
2366
13.4k
}
2367
2368
/* Drop constraints from "context" that are irrelevant for computing
2369
 * the gist of the inequalities "ineq".
2370
 * Inequalities in "ineq" for which the corresponding element of row
2371
 * is set to -1 have already been marked for removal and should be ignored.
2372
 *
2373
 * In particular, drop constraints in variables that are not related
2374
 * to any of the variables involved in "ineq"
2375
 * in the sense that there is no sequence of constraints that connects them.
2376
 *
2377
 * We first mark all variables that appear in "bset" as belonging
2378
 * to a "-1" group and then continue with group_and_drop_irrelevant_constraints.
2379
 */
2380
static __isl_give isl_basic_set *drop_irrelevant_constraints_marked(
2381
  __isl_take isl_basic_set *context, __isl_keep isl_mat *ineq, int *row)
2382
7.86k
{
2383
7.86k
  int *group;
2384
7.86k
  int dim;
2385
7.86k
  int i, j, n;
2386
7.86k
2387
7.86k
  if (
!context || 7.86k
!ineq7.86k
)
2388
0
    return isl_basic_set_free(context);
2389
7.86k
2390
7.86k
  group = alloc_groups(context);
2391
7.86k
2392
7.86k
  if (!group)
2393
0
    return isl_basic_set_free(context);
2394
7.86k
2395
7.86k
  dim = isl_basic_set_dim(context, isl_dim_set);
2396
7.86k
  n = isl_mat_rows(ineq);
2397
48.5k
  for (i = 0; 
i < dim48.5k
;
++i40.6k
)
{40.6k
2398
163k
    for (j = 0; 
j < n163k
;
++j123k
)
{139k
2399
139k
      if (row[j] < 0)
2400
23.1k
        continue;
2401
116k
      
if (116k
!116k
isl_int_is_zero116k
(ineq->row[j][1 + i]))
2402
16.2k
        break;
2403
116k
    }
2404
40.6k
    if (j < n)
2405
16.2k
      group[i] = -1;
2406
40.6k
  }
2407
7.86k
2408
7.86k
  return isl_basic_map_drop_unrelated_constraints(context, group);
2409
7.86k
}
2410
2411
/* Do all "n" entries of "row" contain a negative value?
2412
 */
2413
static int all_neg(int *row, int n)
2414
8.91k
{
2415
8.91k
  int i;
2416
8.91k
2417
14.8k
  for (i = 0; 
i < n14.8k
;
++i5.90k
)
2418
13.7k
    
if (13.7k
row[i] >= 013.7k
)
2419
7.86k
      return 0;
2420
8.91k
2421
1.05k
  return 1;
2422
8.91k
}
2423
2424
/* Update the inequalities in "bset" based on the information in "row"
2425
 * and "tab".
2426
 *
2427
 * In particular, the array "row" contains either -1, meaning that
2428
 * the corresponding inequality of "bset" is redundant, or the index
2429
 * of an inequality in "tab".
2430
 *
2431
 * If the row entry is -1, then drop the inequality.
2432
 * Otherwise, if the constraint is marked redundant in the tableau,
2433
 * then drop the inequality.  Similarly, if it is marked as an equality
2434
 * in the tableau, then turn the inequality into an equality and
2435
 * perform Gaussian elimination.
2436
 */
2437
static __isl_give isl_basic_set *update_ineq(__isl_take isl_basic_set *bset,
2438
  __isl_keep int *row, struct isl_tab *tab)
2439
8.91k
{
2440
8.91k
  int i;
2441
8.91k
  unsigned n_ineq;
2442
8.91k
  unsigned n_eq;
2443
8.91k
  int found_equality = 0;
2444
8.91k
2445
8.91k
  if (!bset)
2446
0
    return NULL;
2447
8.91k
  
if (8.91k
tab && 8.91k
tab->empty7.63k
)
2448
2.18k
    return isl_basic_set_set_to_empty(bset);
2449
8.91k
2450
6.73k
  n_ineq = bset->n_ineq;
2451
24.4k
  for (i = n_ineq - 1; 
i >= 024.4k
;
--i17.7k
)
{17.7k
2452
17.7k
    if (
row[i] < 017.7k
)
{6.37k
2453
6.37k
      if (isl_basic_set_drop_inequality(bset, i) < 0)
2454
0
        return isl_basic_set_free(bset);
2455
6.37k
      continue;
2456
6.37k
    }
2457
11.3k
    
if (11.3k
!tab11.3k
)
2458
408
      continue;
2459
10.9k
    n_eq = tab->n_eq;
2460
10.9k
    if (
isl_tab_is_equality(tab, n_eq + row[i])10.9k
)
{151
2461
151
      isl_basic_map_inequality_to_equality(bset, i);
2462
151
      found_equality = 1;
2463
10.8k
    } else 
if (10.8k
isl_tab_is_redundant(tab, n_eq + row[i])10.8k
)
{244
2464
244
      if (isl_basic_set_drop_inequality(bset, i) < 0)
2465
0
        return isl_basic_set_free(bset);
2466
244
    }
2467
10.9k
  }
2468
6.73k
2469
6.73k
  
if (6.73k
found_equality6.73k
)
2470
144
    bset = isl_basic_set_gauss(bset, NULL);
2471
6.73k
  bset = isl_basic_set_finalize(bset);
2472
6.73k
  return bset;
2473
6.73k
}
2474
2475
/* Update the inequalities in "bset" based on the information in "row"
2476
 * and "tab" and free all arguments (other than "bset").
2477
 */
2478
static __isl_give isl_basic_set *update_ineq_free(
2479
  __isl_take isl_basic_set *bset, __isl_take isl_mat *ineq,
2480
  __isl_take isl_basic_set *context, __isl_take int *row,
2481
  struct isl_tab *tab)
2482
8.91k
{
2483
8.91k
  isl_mat_free(ineq);
2484
8.91k
  isl_basic_set_free(context);
2485
8.91k
2486
8.91k
  bset = update_ineq(bset, row, tab);
2487
8.91k
2488
8.91k
  free(row);
2489
8.91k
  isl_tab_free(tab);
2490
8.91k
  return bset;
2491
8.91k
}
2492
2493
/* Remove all information from bset that is redundant in the context
2494
 * of context.
2495
 * "ineq" contains the (possibly transformed) inequalities of "bset",
2496
 * in the same order.
2497
 * The (explicit) equalities of "bset" are assumed to have been taken
2498
 * into account by the transformation such that only the inequalities
2499
 * are relevant.
2500
 * "context" is assumed not to be empty.
2501
 *
2502
 * "row" keeps track of the constraint index of a "bset" inequality in "tab".
2503
 * A value of -1 means that the inequality is obviously redundant and may
2504
 * not even appear in  "tab".
2505
 *
2506
 * We first mark the inequalities of "bset"
2507
 * that are obviously redundant with respect to some inequality in "context".
2508
 * Then we remove those constraints from "context" that have become
2509
 * irrelevant for computing the gist of "bset".
2510
 * Note that this removal of constraints cannot be replaced by
2511
 * a factorization because factors in "bset" may still be connected
2512
 * to each other through constraints in "context".
2513
 *
2514
 * If there are any inequalities left, we construct a tableau for
2515
 * the context and then add the inequalities of "bset".
2516
 * Before adding these inequalities, we freeze all constraints such that
2517
 * they won't be considered redundant in terms of the constraints of "bset".
2518
 * Then we detect all redundant constraints (among the
2519
 * constraints that weren't frozen), first by checking for redundancy in the
2520
 * the tableau and then by checking if replacing a constraint by its negation
2521
 * would lead to an empty set.  This last step is fairly expensive
2522
 * and could be optimized by more reuse of the tableau.
2523
 * Finally, we update bset according to the results.
2524
 */
2525
static __isl_give isl_basic_set *uset_gist_full(__isl_take isl_basic_set *bset,
2526
  __isl_take isl_mat *ineq, __isl_take isl_basic_set *context)
2527
13.4k
{
2528
13.4k
  int i, r;
2529
13.4k
  int *row = NULL;
2530
13.4k
  isl_ctx *ctx;
2531
13.4k
  isl_basic_set *combined = NULL;
2532
13.4k
  struct isl_tab *tab = NULL;
2533
13.4k
  unsigned n_eq, context_ineq;
2534
13.4k
2535
13.4k
  if (
!bset || 13.4k
!ineq13.4k
||
!context13.4k
)
2536
0
    goto error;
2537
13.4k
2538
13.4k
  
if (13.4k
bset->n_ineq == 0 || 13.4k
isl_basic_set_plain_is_universe(context)9.42k
)
{4.53k
2539
4.53k
    isl_basic_set_free(context);
2540
4.53k
    isl_mat_free(ineq);
2541
4.53k
    return bset;
2542
4.53k
  }
2543
13.4k
2544
8.91k
  ctx = isl_basic_set_get_ctx(context);
2545
8.91k
  row = isl_calloc_array(ctx, int, bset->n_ineq);
2546
8.91k
  if (!row)
2547
0
    goto error;
2548
8.91k
2549
8.91k
  
if (8.91k
mark_shifted_constraints(ineq, context, row) < 08.91k
)
2550
0
    goto error;
2551
8.91k
  
if (8.91k
all_neg(row, bset->n_ineq)8.91k
)
2552
1.05k
    return update_ineq_free(bset, ineq, context, row, NULL);
2553
8.91k
2554
7.86k
  context = drop_irrelevant_constraints_marked(context, ineq, row);
2555
7.86k
  if (!context)
2556
0
    goto error;
2557
7.86k
  
if (7.86k
isl_basic_set_plain_is_universe(context)7.86k
)
2558
228
    return update_ineq_free(bset, ineq, context, row, NULL);
2559
7.86k
2560
7.63k
  n_eq = context->n_eq;
2561
7.63k
  context_ineq = context->n_ineq;
2562
7.63k
  combined = isl_basic_set_cow(isl_basic_set_copy(context));
2563
7.63k
  combined = isl_basic_set_extend_constraints(combined, 0, bset->n_ineq);
2564
7.63k
  tab = isl_tab_from_basic_set(combined, 0);
2565
31.3k
  for (i = 0; 
i < context_ineq31.3k
;
++i23.7k
)
2566
23.7k
    
if (23.7k
isl_tab_freeze_constraint(tab, n_eq + i) < 023.7k
)
2567
0
      goto error;
2568
7.63k
  
if (7.63k
isl_tab_extend_cons(tab, bset->n_ineq) < 07.63k
)
2569
0
    goto error;
2570
7.63k
  r = context_ineq;
2571
26.9k
  for (i = 0; 
i < bset->n_ineq26.9k
;
++i19.2k
)
{19.2k
2572
19.2k
    if (row[i] < 0)
2573
3.35k
      continue;
2574
15.9k
    combined = isl_basic_set_add_ineq(combined, ineq->row[i]);
2575
15.9k
    if (isl_tab_add_ineq(tab, ineq->row[i]) < 0)
2576
0
      goto error;
2577
15.9k
    row[i] = r++;
2578
15.9k
  }
2579
7.63k
  
if (7.63k
isl_tab_detect_implicit_equalities(tab) < 07.63k
)
2580
0
    goto error;
2581
7.63k
  
if (7.63k
isl_tab_detect_redundant(tab) < 07.63k
)
2582
0
    goto error;
2583
26.9k
  
for (i = bset->n_ineq - 1; 7.63k
i >= 026.9k
;
--i19.2k
)
{19.2k
2584
19.2k
    isl_basic_set *test;
2585
19.2k
    int is_empty;
2586
19.2k
2587
19.2k
    if (row[i] < 0)
2588
3.35k
      continue;
2589
15.9k
    r = row[i];
2590
15.9k
    if (tab->con[n_eq + r].is_redundant)
2591
385
      continue;
2592
15.5k
    test = isl_basic_set_dup(combined);
2593
15.5k
    if (isl_inequality_negate(test, r) < 0)
2594
0
      test = isl_basic_set_free(test);
2595
15.5k
    test = isl_basic_set_update_from_tab(test, tab);
2596
15.5k
    is_empty = isl_basic_set_is_empty(test);
2597
15.5k
    isl_basic_set_free(test);
2598
15.5k
    if (is_empty < 0)
2599
0
      goto error;
2600
15.5k
    
if (15.5k
is_empty15.5k
)
2601
4.95k
      tab->con[n_eq + r].is_redundant = 1;
2602
15.5k
  }
2603
7.63k
  bset = update_ineq_free(bset, ineq, context, row, tab);
2604
7.63k
  if (
bset7.63k
)
{7.63k
2605
7.63k
    ISL_F_SET(bset, ISL_BASIC_SET_NO_IMPLICIT);
2606
7.63k
    ISL_F_SET(bset, ISL_BASIC_SET_NO_REDUNDANT);
2607
7.63k
  }
2608
7.63k
2609
7.63k
  isl_basic_set_free(combined);
2610
7.63k
  return bset;
2611
0
error:
2612
0
  free(row);
2613
0
  isl_mat_free(ineq);
2614
0
  isl_tab_free(tab);
2615
0
  isl_basic_set_free(combined);
2616
0
  isl_basic_set_free(context);
2617
0
  isl_basic_set_free(bset);
2618
0
  return NULL;
2619
7.63k
}
2620
2621
/* Extract the inequalities of "bset" as an isl_mat.
2622
 */
2623
static __isl_give isl_mat *extract_ineq(__isl_keep isl_basic_set *bset)
2624
13.4k
{
2625
13.4k
  unsigned total;
2626
13.4k
  isl_ctx *ctx;
2627
13.4k
  isl_mat *ineq;
2628
13.4k
2629
13.4k
  if (!bset)
2630
0
    return NULL;
2631
13.4k
2632
13.4k
  ctx = isl_basic_set_get_ctx(bset);
2633
13.4k
  total = isl_basic_set_total_dim(bset);
2634
13.4k
  ineq = isl_mat_sub_alloc6(ctx, bset->ineq, 0, bset->n_ineq,
2635
13.4k
            0, 1 + total);
2636
13.4k
2637
13.4k
  return ineq;
2638
13.4k
}
2639
2640
/* Remove all information from "bset" that is redundant in the context
2641
 * of "context", for the case where both "bset" and "context" are
2642
 * full-dimensional.
2643
 */
2644
static __isl_give isl_basic_set *uset_gist_uncompressed(
2645
  __isl_take isl_basic_set *bset, __isl_take isl_basic_set *context)
2646
7.02k
{
2647
7.02k
  isl_mat *ineq;
2648
7.02k
2649
7.02k
  ineq = extract_ineq(bset);
2650
7.02k
  return uset_gist_full(bset, ineq, context);
2651
7.02k
}
2652
2653
/* Remove all information from "bset" that is redundant in the context
2654
 * of "context", for the case where the combined equalities of
2655
 * "bset" and "context" allow for a compression that can be obtained
2656
 * by preapplication of "T".
2657
 *
2658
 * "bset" itself is not transformed by "T".  Instead, the inequalities
2659
 * are extracted from "bset" and those are transformed by "T".
2660
 * uset_gist_full then determines which of the transformed inequalities
2661
 * are redundant with respect to the transformed "context" and removes
2662
 * the corresponding inequalities from "bset".
2663
 *
2664
 * After preapplying "T" to the inequalities, any common factor is
2665
 * removed from the coefficients.  If this results in a tightening
2666
 * of the constant term, then the same tightening is applied to
2667
 * the corresponding untransformed inequality in "bset".
2668
 * That is, if after plugging in T, a constraint f(x) >= 0 is of the form
2669
 *
2670
 *  g f'(x) + r >= 0
2671
 *
2672
 * with 0 <= r < g, then it is equivalent to
2673
 *
2674
 *  f'(x) >= 0
2675
 *
2676
 * This means that f(x) >= 0 is equivalent to f(x) - r >= 0 in the affine
2677
 * subspace compressed by T since the latter would be transformed to
2678
 *
2679
 *  g f'(x) >= 0
2680
 */
2681
static __isl_give isl_basic_set *uset_gist_compressed(
2682
  __isl_take isl_basic_set *bset, __isl_take isl_basic_set *context,
2683
  __isl_take isl_mat *T)
2684
6.42k
{
2685
6.42k
  isl_ctx *ctx;
2686
6.42k
  isl_mat *ineq;
2687
6.42k
  int i, n_row, n_col;
2688
6.42k
  isl_int rem;
2689
6.42k
2690
6.42k
  ineq = extract_ineq(bset);
2691
6.42k
  ineq = isl_mat_product(ineq, isl_mat_copy(T));
2692
6.42k
  context = isl_basic_set_preimage(context, T);
2693
6.42k
2694
6.42k
  if (
!ineq || 6.42k
!context6.42k
)
2695
0
    goto error;
2696
6.42k
  
if (6.42k
isl_basic_set_plain_is_empty(context)6.42k
)
{7
2697
7
    isl_mat_free(ineq);
2698
7
    isl_basic_set_free(context);
2699
7
    return isl_basic_set_set_to_empty(bset);
2700
7
  }
2701
6.42k
2702
6.42k
  ctx = isl_mat_get_ctx(ineq);
2703
6.42k
  n_row = isl_mat_rows(ineq);
2704
6.42k
  n_col = isl_mat_cols(ineq);
2705
6.42k
  isl_int_init(rem);
2706
12.7k
  for (i = 0; 
i < n_row12.7k
;
++i6.30k
)
{6.30k
2707
6.30k
    isl_seq_gcd(ineq->row[i] + 1, n_col - 1, &ctx->normalize_gcd);
2708
6.30k
    if (isl_int_is_zero(ctx->normalize_gcd))
2709
17
      continue;
2710
6.29k
    
if (6.29k
isl_int_is_one6.29k
(ctx->normalize_gcd))
2711
5.90k
      continue;
2712
389
    isl_seq_scale_down(ineq->row[i] + 1, ineq->row[i] + 1,
2713
389
            ctx->normalize_gcd, n_col - 1);
2714
389
    isl_int_fdiv_r(rem, ineq->row[i][0], ctx->normalize_gcd);
2715
389
    isl_int_fdiv_q(ineq->row[i][0],
2716
389
        ineq->row[i][0], ctx->normalize_gcd);
2717
389
    if (isl_int_is_zero(rem))
2718
281
      continue;
2719
108
    bset = isl_basic_set_cow(bset);
2720
108
    if (!bset)
2721
0
      break;
2722
108
    
isl_int_sub108
(bset->ineq[i][0], bset->ineq[i][0], rem);108
2723
108
  }
2724
6.42k
  isl_int_clear(rem);
2725
6.42k
2726
6.42k
  return uset_gist_full(bset, ineq, context);
2727
0
error:
2728
0
  isl_mat_free(ineq);
2729
0
  isl_basic_set_free(context);
2730
0
  isl_basic_set_free(bset);
2731
0
  return NULL;
2732
6.42k
}
2733
2734
/* Project "bset" onto the variables that are involved in "template".
2735
 */
2736
static __isl_give isl_basic_set *project_onto_involved(
2737
  __isl_take isl_basic_set *bset, __isl_keep isl_basic_set *template)
2738
6.42k
{
2739
6.42k
  int i, n;
2740
6.42k
2741
6.42k
  if (
!bset || 6.42k
!template6.42k
)
2742
0
    return isl_basic_set_free(bset);
2743
6.42k
2744
6.42k
  n = isl_basic_set_dim(template, isl_dim_set);
2745
6.42k
2746
40.5k
  for (i = 0; 
i < n40.5k
;
++i34.0k
)
{34.0k
2747
34.0k
    isl_bool involved;
2748
34.0k
2749
34.0k
    involved = isl_basic_set_involves_dims(template,
2750
34.0k
              isl_dim_set, i, 1);
2751
34.0k
    if (involved < 0)
2752
0
      return isl_basic_set_free(bset);
2753
34.0k
    
if (34.0k
involved34.0k
)
2754
19.6k
      continue;
2755
14.4k
    bset = isl_basic_set_eliminate_vars(bset, i, 1);
2756
14.4k
  }
2757
6.42k
2758
6.42k
  return bset;
2759
6.42k
}
2760
2761
/* Remove all information from bset that is redundant in the context
2762
 * of context.  In particular, equalities that are linear combinations
2763
 * of those in context are removed.  Then the inequalities that are
2764
 * redundant in the context of the equalities and inequalities of
2765
 * context are removed.
2766
 *
2767
 * First of all, we drop those constraints from "context"
2768
 * that are irrelevant for computing the gist of "bset".
2769
 * Alternatively, we could factorize the intersection of "context" and "bset".
2770
 *
2771
 * We first compute the intersection of the integer affine hulls
2772
 * of "bset" and "context",
2773
 * compute the gist inside this intersection and then reduce
2774
 * the constraints with respect to the equalities of the context
2775
 * that only involve variables already involved in the input.
2776
 *
2777
 * If two constraints are mutually redundant, then uset_gist_full
2778
 * will remove the second of those constraints.  We therefore first
2779
 * sort the constraints so that constraints not involving existentially
2780
 * quantified variables are given precedence over those that do.
2781
 * We have to perform this sorting before the variable compression,
2782
 * because that may effect the order of the variables.
2783
 */
2784
static __isl_give isl_basic_set *uset_gist(__isl_take isl_basic_set *bset,
2785
  __isl_take isl_basic_set *context)
2786
13.4k
{
2787
13.4k
  isl_mat *eq;
2788
13.4k
  isl_mat *T;
2789
13.4k
  isl_basic_set *aff;
2790
13.4k
  isl_basic_set *aff_context;
2791
13.4k
  unsigned total;
2792
13.4k
2793
13.4k
  if (
!bset || 13.4k
!context13.4k
)
2794
0
    goto error;
2795
13.4k
2796
13.4k
  context = drop_irrelevant_constraints(context, bset);
2797
13.4k
2798
13.4k
  bset = isl_basic_set_detect_equalities(bset);
2799
13.4k
  aff = isl_basic_set_copy(bset);
2800
13.4k
  aff = isl_basic_set_plain_affine_hull(aff);
2801
13.4k
  context = isl_basic_set_detect_equalities(context);
2802
13.4k
  aff_context = isl_basic_set_copy(context);
2803
13.4k
  aff_context = isl_basic_set_plain_affine_hull(aff_context);
2804
13.4k
  aff = isl_basic_set_intersect(aff, aff_context);
2805
13.4k
  if (!aff)
2806
0
    goto error;
2807
13.4k
  
if (13.4k
isl_basic_set_plain_is_empty(aff)13.4k
)
{1
2808
1
    isl_basic_set_free(bset);
2809
1
    isl_basic_set_free(context);
2810
1
    return aff;
2811
1
  }
2812
13.4k
  bset = isl_basic_set_sort_constraints(bset);
2813
13.4k
  if (
aff->n_eq == 013.4k
)
{7.02k
2814
7.02k
    isl_basic_set_free(aff);
2815
7.02k
    return uset_gist_uncompressed(bset, context);
2816
7.02k
  }
2817
6.42k
  total = isl_basic_set_total_dim(bset);
2818
6.42k
  eq = isl_mat_sub_alloc6(bset->ctx, aff->eq, 0, aff->n_eq, 0, 1 + total);
2819
6.42k
  eq = isl_mat_cow(eq);
2820
6.42k
  T = isl_mat_variable_compression(eq, NULL);
2821
6.42k
  isl_basic_set_free(aff);
2822
6.42k
  if (
T && 6.42k
T->n_col == 06.42k
)
{0
2823
0
    isl_mat_free(T);
2824
0
    isl_basic_set_free(context);
2825
0
    return isl_basic_set_set_to_empty(bset);
2826
0
  }
2827
6.42k
2828
6.42k
  aff_context = isl_basic_set_affine_hull(isl_basic_set_copy(context));
2829
6.42k
  aff_context = project_onto_involved(aff_context, bset);
2830
6.42k
2831
6.42k
  bset = uset_gist_compressed(bset, context, T);
2832
6.42k
  bset = isl_basic_set_reduce_using_equalities(bset, aff_context);
2833
6.42k
2834
6.42k
  if (
bset6.42k
)
{6.42k
2835
6.42k
    ISL_F_SET(bset, ISL_BASIC_SET_NO_IMPLICIT);
2836
6.42k
    ISL_F_SET(bset, ISL_BASIC_SET_NO_REDUNDANT);
2837
6.42k
  }
2838
6.42k
2839
6.42k
  return bset;
2840
0
error:
2841
0
  isl_basic_set_free(bset);
2842
0
  isl_basic_set_free(context);
2843
0
  return NULL;
2844
6.42k
}
2845
2846
/* Return the number of equality constraints in "bmap" that involve
2847
 * local variables.  This function assumes that Gaussian elimination
2848
 * has been applied to the equality constraints.
2849
 */
2850
static int n_div_eq(__isl_keep isl_basic_map *bmap)
2851
806
{
2852
806
  int i;
2853
806
  int total, n_div;
2854
806
2855
806
  if (!bmap)
2856
0
    return -1;
2857
806
2858
806
  
if (806
bmap->n_eq == 0806
)
2859
287
    return 0;
2860
806
2861
519
  total = isl_basic_map_dim(bmap, isl_dim_all);
2862
519
  n_div = isl_basic_map_dim(bmap, isl_dim_div);
2863
519
  total -= n_div;
2864
519
2865
907
  for (i = 0; 
i < bmap->n_eq907
;
++i388
)
2866
594
    
if (594
isl_seq_first_non_zero(bmap->eq[i] + 1 + total,594
2867
594
              n_div) == -1)
2868
206
      return i;
2869
519
2870
313
  return bmap->n_eq;
2871
519
}
2872
2873
/* Construct a basic map in "space" defined by the equality constraints in "eq".
2874
 * The constraints are assumed not to involve any local variables.
2875
 */
2876
static __isl_give isl_basic_map *basic_map_from_equalities(
2877
  __isl_take isl_space *space, __isl_take isl_mat *eq)
2878
2
{
2879
2
  int i, k;
2880
2
  isl_basic_map *bmap = NULL;
2881
2
2882
2
  if (
!space || 2
!eq2
)
2883
0
    goto error;
2884
2
2885
2
  
if (2
1 + isl_space_dim(space, isl_dim_all) != eq->n_col2
)
2886
0
    isl_die(isl_space_get_ctx(space), isl_error_internal,
2887
2
      "unexpected number of columns", goto error);
2888
2
2889
2
  bmap = isl_basic_map_alloc_space(isl_space_copy(space),
2890
2
              0, eq->n_row, 0);
2891
6
  for (i = 0; 
i < eq->n_row6
;
++i4
)
{4
2892
4
    k = isl_basic_map_alloc_equality(bmap);
2893
4
    if (k < 0)
2894
0
      goto error;
2895
4
    isl_seq_cpy(bmap->eq[k], eq->row[i], eq->n_col);
2896
4
  }
2897
2
2898
2
  isl_space_free(space);
2899
2
  isl_mat_free(eq);
2900
2
  return bmap;
2901
0
error:
2902
0
  isl_space_free(space);
2903
0
  isl_mat_free(eq);
2904
0
  isl_basic_map_free(bmap);
2905
0
  return NULL;
2906
2
}
2907
2908
/* Construct and return a variable compression based on the equality
2909
 * constraints in "bmap1" and "bmap2" that do not involve the local variables.
2910
 * "n1" is the number of (initial) equality constraints in "bmap1"
2911
 * that do involve local variables.
2912
 * "n2" is the number of (initial) equality constraints in "bmap2"
2913
 * that do involve local variables.
2914
 * "total" is the total number of other variables.
2915
 * This function assumes that Gaussian elimination
2916
 * has been applied to the equality constraints in both "bmap1" and "bmap2"
2917
 * such that the equality constraints not involving local variables
2918
 * are those that start at "n1" or "n2".
2919
 *
2920
 * If either of "bmap1" and "bmap2" does not have such equality constraints,
2921
 * then simply compute the compression based on the equality constraints
2922
 * in the other basic map.
2923
 * Otherwise, combine the equality constraints from both into a new
2924
 * basic map such that Gaussian elimination can be applied to this combination
2925
 * and then construct a variable compression from the resulting
2926
 * equality constraints.
2927
 */
2928
static __isl_give isl_mat *combined_variable_compression(
2929
  __isl_keep isl_basic_map *bmap1, int n1,
2930
  __isl_keep isl_basic_map *bmap2, int n2, int total)
2931
7
{
2932
7
  isl_ctx *ctx;
2933
7
  isl_mat *E1, *E2, *V;
2934
7
  isl_basic_map *bmap;
2935
7
2936
7
  ctx = isl_basic_map_get_ctx(bmap1);
2937
7
  if (
bmap1->n_eq == n17
)
{3
2938
3
    E2 = isl_mat_sub_alloc6(ctx, bmap2->eq,
2939
3
          n2, bmap2->n_eq - n2, 0, 1 + total);
2940
3
    return isl_mat_variable_compression(E2, NULL);
2941
3
  }
2942
4
  
if (4
bmap2->n_eq == n24
)
{2
2943
2
    E1 = isl_mat_sub_alloc6(ctx, bmap1->eq,
2944
2
          n1, bmap1->n_eq - n1, 0, 1 + total);
2945
2
    return isl_mat_variable_compression(E1, NULL);
2946
2
  }
2947
2
  E1 = isl_mat_sub_alloc6(ctx, bmap1->eq,
2948
2
        n1, bmap1->n_eq - n1, 0, 1 + total);
2949
2
  E2 = isl_mat_sub_alloc6(ctx, bmap2->eq,
2950
2
        n2, bmap2->n_eq - n2, 0, 1 + total);
2951
2
  E1 = isl_mat_concat(E1, E2);
2952
2
  bmap = basic_map_from_equalities(isl_basic_map_get_space(bmap1), E1);
2953
2
  bmap = isl_basic_map_gauss(bmap, NULL);
2954
2
  if (!bmap)
2955
0
    return NULL;
2956
2
  E1 = isl_mat_sub_alloc6(ctx, bmap->eq, 0, bmap->n_eq, 0, 1 + total);
2957
2
  V = isl_mat_variable_compression(E1, NULL);
2958
2
  isl_basic_map_free(bmap);
2959
2
2960
2
  return V;
2961
2
}
2962
2963
/* Extract the stride constraints from "bmap", compressed
2964
 * with respect to both the stride constraints in "context" and
2965
 * the remaining equality constraints in both "bmap" and "context".
2966
 * "bmap_n_eq" is the number of (initial) stride constraints in "bmap".
2967
 * "context_n_eq" is the number of (initial) stride constraints in "context".
2968
 *
2969
 * Let x be all variables in "bmap" (and "context") other than the local
2970
 * variables.  First compute a variable compression
2971
 *
2972
 *  x = V x'
2973
 *
2974
 * based on the non-stride equality constraints in "bmap" and "context".
2975
 * Consider the stride constraints of "context",
2976
 *
2977
 *  A(x) + B(y) = 0
2978
 *
2979
 * with y the local variables and plug in the variable compression,
2980
 * resulting in
2981
 *
2982
 *  A(V x') + B(y) = 0
2983
 *
2984
 * Use these constraints to compute a parameter compression on x'
2985
 *
2986
 *  x' = T x''
2987
 *
2988
 * Now consider the stride constraints of "bmap"
2989
 *
2990
 *  C(x) + D(y) = 0
2991
 *
2992
 * and plug in x = V*T x''.
2993
 * That is, return A = [C*V*T D].
2994
 */
2995
static __isl_give isl_mat *extract_compressed_stride_constraints(
2996
  __isl_keep isl_basic_map *bmap, int bmap_n_eq,
2997
  __isl_keep isl_basic_map *context, int context_n_eq)
2998
7
{
2999
7
  int total, n_div;
3000
7
  isl_ctx *ctx;
3001
7
  isl_mat *A, *B, *T, *V;
3002
7
3003
7
  total = isl_basic_map_dim(context, isl_dim_all);
3004
7
  n_div = isl_basic_map_dim(context, isl_dim_div);
3005
7
  total -= n_div;
3006
7
3007
7
  ctx = isl_basic_map_get_ctx(bmap);
3008
7
3009
7
  V = combined_variable_compression(bmap, bmap_n_eq,
3010
7
            context, context_n_eq, total);
3011
7
3012
7
  A = isl_mat_sub_alloc6(ctx, context->eq, 0, context_n_eq, 0, 1 + total);
3013
7
  B = isl_mat_sub_alloc6(ctx, context->eq,
3014
7
        0, context_n_eq, 1 + total, n_div);
3015
7
  A = isl_mat_product(A, isl_mat_copy(V));
3016
7
  T = isl_mat_parameter_compression_ext(A, B);
3017
7
  T = isl_mat_product(V, T);
3018
7
3019
7
  n_div = isl_basic_map_dim(bmap, isl_dim_div);
3020
7
  T = isl_mat_diagonal(T, isl_mat_identity(ctx, n_div));
3021
7
3022
7
  A = isl_mat_sub_alloc6(ctx, bmap->eq,
3023
7
        0, bmap_n_eq, 0, 1 + total + n_div);
3024
7
  A = isl_mat_product(A, T);
3025
7
3026
7
  return A;
3027
7
}
3028
3029
/* Remove the prime factors from *g that have an exponent that
3030
 * is strictly smaller than the exponent in "c".
3031
 * All exponents in *g are known to be smaller than or equal
3032
 * to those in "c".
3033
 *
3034
 * That is, if *g is equal to
3035
 *
3036
 *  p_1^{e_1} p_2^{e_2} ... p_n^{e_n}
3037
 *
3038
 * and "c" is equal to
3039
 *
3040
 *  p_1^{f_1} p_2^{f_2} ... p_n^{f_n}
3041
 *
3042
 * then update *g to
3043
 *
3044
 *  p_1^{e_1 * (e_1 = f_1)} p_2^{e_2 * (e_2 = f_2)} ...
3045
 *    p_n^{e_n * (e_n = f_n)}
3046
 *
3047
 * If e_i = f_i, then c / *g does not have any p_i factors and therefore
3048
 * neither does the gcd of *g and c / *g.
3049
 * If e_i < f_i, then the gcd of *g and c / *g has a positive
3050
 * power min(e_i, s_i) of p_i with s_i = f_i - e_i among its factors.
3051
 * Dividing *g by this gcd therefore strictly reduces the exponent
3052
 * of the prime factors that need to be removed, while leaving the
3053
 * other prime factors untouched.
3054
 * Repeating this process until gcd(*g, c / *g) = 1 therefore
3055
 * removes all undesired factors, without removing any others.
3056
 */
3057
static void remove_incomplete_powers(isl_int *g, isl_int c)
3058
6
{
3059
6
  isl_int t;
3060
6
3061
6
  isl_int_init(t);
3062
9
  for (;;) {
3063
9
    isl_int_divexact(t, c, *g);
3064
9
    isl_int_gcd(t, t, *g);
3065
9
    if (isl_int_is_one(t))
3066
6
      break;
3067
3
    
isl_int_divexact3
(*g, *g, t);3
3068
3
  }
3069
6
  isl_int_clear(t);
3070
6
}
3071
3072
/* Reduce the "n" stride constraints in "bmap" based on a copy "A"
3073
 * of the same stride constraints in a compressed space that exploits
3074
 * all equalities in the context and the other equalities in "bmap".
3075
 *
3076
 * If the stride constraints of "bmap" are of the form
3077
 *
3078
 *  C(x) + D(y) = 0
3079
 *
3080
 * then A is of the form
3081
 *
3082
 *  B(x') + D(y) = 0
3083
 *
3084
 * If any of these constraints involves only a single local variable y,
3085
 * then the constraint appears as
3086
 *
3087
 *  f(x) + m y_i = 0
3088
 *
3089
 * in "bmap" and as
3090
 *
3091
 *  h(x') + m y_i = 0
3092
 *
3093
 * in "A".
3094
 *
3095
 * Let g be the gcd of m and the coefficients of h.
3096
 * Then, in particular, g is a divisor of the coefficients of h and
3097
 *
3098
 *  f(x) = h(x')
3099
 *
3100
 * is known to be a multiple of g.
3101
 * If some prime factor in m appears with the same exponent in g,
3102
 * then it can be removed from m because f(x) is already known
3103
 * to be a multiple of g and therefore in particular of this power
3104
 * of the prime factors.
3105
 * Prime factors that appear with a smaller exponent in g cannot
3106
 * be removed from m.
3107
 * Let g' be the divisor of g containing all prime factors that
3108
 * appear with the same exponent in m and g, then
3109
 *
3110
 *  f(x) + m y_i = 0
3111
 *
3112
 * can be replaced by
3113
 *
3114
 *  f(x) + m/g' y_i' = 0
3115
 *
3116
 * Note that (if g' != 1) this changes the explicit representation
3117
 * of y_i to that of y_i', so the integer division at position i
3118
 * is marked unknown and later recomputed by a call to
3119
 * isl_basic_map_gauss.
3120
 */
3121
static __isl_give isl_basic_map *reduce_stride_constraints(
3122
  __isl_take isl_basic_map *bmap, int n, __isl_keep isl_mat *A)
3123
7
{
3124
7
  int i;
3125
7
  int total, n_div;
3126
7
  int any = 0;
3127
7
  isl_int gcd;
3128
7
3129
7
  if (
!bmap || 7
!A7
)
3130
0
    return isl_basic_map_free(bmap);
3131
7
3132
7
  total = isl_basic_map_dim(bmap, isl_dim_all);
3133
7
  n_div = isl_basic_map_dim(bmap, isl_dim_div);
3134
7
  total -= n_div;
3135
7
3136
7
  isl_int_init(gcd);
3137
14
  for (i = 0; 
i < n14
;
++i7
)
{7
3138
7
    int div;
3139
7
3140
7
    div = isl_seq_first_non_zero(bmap->eq[i] + 1 + total, n_div);
3141
7
    if (div < 0)
3142
0
      isl_die(isl_basic_map_get_ctx(bmap), isl_error_internal,
3143
7
        "equality constraints modified unexpectedly",
3144
7
        goto error);
3145
7
    
if (7
isl_seq_first_non_zero(bmap->eq[i] + 1 + total + div + 1,7
3146
7
            n_div - div - 1) != -1)
3147
0
      continue;
3148
7
    
if (7
isl_mat_row_gcd(A, i, &gcd) < 07
)
3149
0
      goto error;
3150
7
    
if (7
isl_int_is_one7
(gcd))
3151
1
      continue;
3152
6
    remove_incomplete_powers(&gcd, bmap->eq[i][1 + total + div]);
3153
6
    if (isl_int_is_one(gcd))
3154
1
      continue;
3155
5
    
isl_int_divexact5
(bmap->eq[i][1 + total + div],5
3156
5
        bmap->eq[i][1 + total + div], gcd);
3157
5
    bmap = isl_basic_map_mark_div_unknown(bmap, div);
3158
5
    if (!bmap)
3159
0
      goto error;
3160
5
    any = 1;
3161
5
  }
3162
7
  
isl_int_clear7
(gcd);7
3163
7
3164
7
  if (any)
3165
5
    bmap = isl_basic_map_gauss(bmap, NULL);
3166
7
3167
7
  return bmap;
3168
0
error:
3169
0
  isl_int_clear(gcd);
3170
0
  isl_basic_map_free(bmap);
3171
0
  return NULL;
3172
7
}
3173
3174
/* Simplify the stride constraints in "bmap" based on
3175
 * the remaining equality constraints in "bmap" and all equality
3176
 * constraints in "context".
3177
 * Only do this if both "bmap" and "context" have stride constraints.
3178
 *
3179
 * First extract a copy of the stride constraints in "bmap" in a compressed
3180
 * space exploiting all the other equality constraints and then
3181
 * use this compressed copy to simplify the original stride constraints.
3182
 */
3183
static __isl_give isl_basic_map *gist_strides(__isl_take isl_basic_map *bmap,
3184
  __isl_keep isl_basic_map *context)
3185
403
{
3186
403
  int bmap_n_eq, context_n_eq;
3187
403
  isl_mat *A;
3188
403
3189
403
  if (
!bmap || 403
!context403
)
3190
0
    return isl_basic_map_free(bmap);
3191
403
3192
403
  bmap_n_eq = n_div_eq(bmap);
3193
403
  context_n_eq = n_div_eq(context);
3194
403
3195
403
  if (
bmap_n_eq < 0 || 403
context_n_eq < 0403
)
3196
0
    return isl_basic_map_free(bmap);
3197
403
  
if (403
bmap_n_eq == 0 || 403
context_n_eq == 0245
)
3198
396
    return bmap;
3199
403
3200
7
  A = extract_compressed_stride_constraints(bmap, bmap_n_eq,
3201
7
                context, context_n_eq);
3202
7
  bmap = reduce_stride_constraints(bmap, bmap_n_eq, A);
3203
7
3204
7
  isl_mat_free(A);
3205
7
3206
7
  return bmap;
3207
403
}
3208
3209
/* Return a basic map that has the same intersection with "context" as "bmap"
3210
 * and that is as "simple" as possible.
3211
 *
3212
 * The core computation is performed on the pure constraints.
3213
 * When we add back the meaning of the integer divisions, we need
3214
 * to (re)introduce the div constraints.  If we happen to have
3215
 * discovered that some of these integer divisions are equal to
3216
 * some affine combination of other variables, then these div
3217
 * constraints may end up getting simplified in terms of the equalities,
3218
 * resulting in extra inequalities on the other variables that
3219
 * may have been removed already or that may not even have been
3220
 * part of the input.  We try and remove those constraints of
3221
 * this form that are most obviously redundant with respect to
3222
 * the context.  We also remove those div constraints that are
3223
 * redundant with respect to the other constraints in the result.
3224
 *
3225
 * The stride constraints among the equality constraints in "bmap" are
3226
 * also simplified with respecting to the other equality constraints
3227
 * in "bmap" and with respect to all equality constraints in "context".
3228
 */
3229
struct isl_basic_map *isl_basic_map_gist(struct isl_basic_map *bmap,
3230
  struct isl_basic_map *context)
3231
14.7k
{
3232
14.7k
  isl_basic_set *bset, *eq;
3233
14.7k
  isl_basic_map *eq_bmap;
3234
14.7k
  unsigned total, n_div, extra, n_eq, n_ineq;
3235
14.7k
3236
14.7k
  if (
!bmap || 14.7k
!context14.7k
)
3237
0
    goto error;
3238
14.7k
3239
14.7k
  
if (14.7k
isl_basic_map_plain_is_universe(bmap)14.7k
)
{1.24k
3240
1.24k
    isl_basic_map_free(context);
3241
1.24k
    return bmap;
3242
1.24k
  }
3243
13.4k
  
if (13.4k
isl_basic_map_plain_is_empty(context)13.4k
)
{0
3244
0
    isl_space *space = isl_basic_map_get_space(bmap);
3245
0
    isl_basic_map_free(bmap);
3246
0
    isl_basic_map_free(context);
3247
0
    return isl_basic_map_universe(space);
3248
0
  }
3249
13.4k
  
if (13.4k
isl_basic_map_plain_is_empty(bmap)13.4k
)
{0
3250
0
    isl_basic_map_free(context);
3251
0
    return bmap;
3252
0
  }
3253
13.4k
3254
13.4k
  bmap = isl_basic_map_remove_redundancies(bmap);
3255
13.4k
  context = isl_basic_map_remove_redundancies(context);
3256
13.4k
  if (!context)
3257
0
    goto error;
3258
13.4k
3259
13.4k
  context = isl_basic_map_align_divs(context, bmap);
3260
13.4k
  n_div = isl_basic_map_dim(context, isl_dim_div);
3261
13.4k
  total = isl_basic_map_dim(bmap, isl_dim_all);
3262
13.4k
  extra = n_div - isl_basic_map_dim(bmap, isl_dim_div);
3263
13.4k
3264
13.4k
  bset = isl_basic_map_underlying_set(isl_basic_map_copy(bmap));
3265
13.4k
  bset = isl_basic_set_add_dims(bset, isl_dim_set, extra);
3266
13.4k
  bset = uset_gist(bset,
3267
13.4k
        isl_basic_map_underlying_set(isl_basic_map_copy(context)));
3268
13.4k
  bset = isl_basic_set_project_out(bset, isl_dim_set, total, extra);
3269
13.4k
3270
13.4k
  if (
!bset || 13.4k
bset->n_eq == 013.4k
||
n_div == 08.15k
||
3271
13.0k
      
isl_basic_set_plain_is_empty(bset)504
)
{13.0k
3272
13.0k
    isl_basic_map_free(context);
3273
13.0k
    return isl_basic_map_overlying_set(bset, bmap);
3274
13.0k
  }
3275
13.4k
3276
403
  n_eq = bset->n_eq;
3277
403
  n_ineq = bset->n_ineq;
3278
403
  eq = isl_basic_set_copy(bset);
3279
403
  eq = isl_basic_set_cow(eq);
3280
403
  if (isl_basic_set_free_inequality(eq, n_ineq) < 0)
3281
0
    eq = isl_basic_set_free(eq);
3282
403
  if (isl_basic_set_free_equality(bset, n_eq) < 0)
3283
0
    bset = isl_basic_set_free(bset);
3284
403
3285
403
  eq_bmap = isl_basic_map_overlying_set(eq, isl_basic_map_copy(bmap));
3286
403
  eq_bmap = gist_strides(eq_bmap, context);
3287
403
  eq_bmap = isl_basic_map_remove_shifted_constraints(eq_bmap, context);
3288
403
  bmap = isl_basic_map_overlying_set(bset, bmap);
3289
403
  bmap = isl_basic_map_intersect(bmap, eq_bmap);
3290
403
  bmap = isl_basic_map_remove_redundancies(bmap);
3291
403
3292
403
  return bmap;
3293
0
error:
3294
0
  isl_basic_map_free(bmap);
3295
0
  isl_basic_map_free(context);
3296
0
  return NULL;
3297
13.4k
}
3298
3299
/*
3300
 * Assumes context has no implicit divs.
3301
 */
3302
__isl_give isl_map *isl_map_gist_basic_map(__isl_take isl_map *map,
3303
  __isl_take isl_basic_map *context)
3304
13.5k
{
3305
13.5k
  int i;
3306
13.5k
3307
13.5k
  if (
!map || 13.5k
!context13.5k
)
3308
0
    goto error;
3309
13.5k
3310
13.5k
  
if (13.5k
isl_basic_map_plain_is_empty(context)13.5k
)
{0
3311
0
    isl_space *space = isl_map_get_space(map);
3312
0
    isl_map_free(map);
3313
0
    isl_basic_map_free(context);
3314
0
    return isl_map_universe(space);
3315
0
  }
3316
13.5k
3317
13.5k
  context = isl_basic_map_remove_redundancies(context);
3318
13.5k
  map = isl_map_cow(map);
3319
13.5k
  if (
!map || 13.5k
!context13.5k
)
3320
0
    goto error;
3321
13.5k
  
isl_assert13.5k
(map->ctx, isl_space_is_equal(map->dim, context->dim), goto error);13.5k
3322
13.5k
  map = isl_map_compute_divs(map);
3323
13.5k
  if (!map)
3324
0
    goto error;
3325
27.5k
  
for (i = map->n - 1; 13.5k
i >= 027.5k
;
--i13.9k
)
{13.9k
3326
13.9k
    map->p[i] = isl_basic_map_gist(map->p[i],
3327
13.9k
            isl_basic_map_copy(context));
3328
13.9k
    if (!map->p[i])
3329
0
      goto error;
3330
13.9k
    
if (13.9k
isl_basic_map_plain_is_empty(map->p[i])13.9k
)
{2.20k
3331
2.20k
      isl_basic_map_free(map->p[i]);
3332
2.20k
      if (i != map->n - 1)
3333
934
        map->p[i] = map->p[map->n - 1];
3334
2.20k
      map->n--;
3335
2.20k
    }
3336
13.9k
  }
3337
13.5k
  isl_basic_map_free(context);
3338
13.5k
  ISL_F_CLR(map, ISL_MAP_NORMALIZED);
3339
13.5k
  return map;
3340
0
error:
3341
0
  isl_map_free(map);
3342
0
  isl_basic_map_free(context);
3343
0
  return NULL;
3344
13.5k
}
3345
3346
/* Drop all inequalities from "bmap" that also appear in "context".
3347
 * "context" is assumed to have only known local variables and
3348
 * the initial local variables of "bmap" are assumed to be the same
3349
 * as those of "context".
3350
 * The constraints of both "bmap" and "context" are assumed
3351
 * to have been sorted using isl_basic_map_sort_constraints.
3352
 *
3353
 * Run through the inequality constraints of "bmap" and "context"
3354
 * in sorted order.
3355
 * If a constraint of "bmap" involves variables not in "context",
3356
 * then it cannot appear in "context".
3357
 * If a matching constraint is found, it is removed from "bmap".
3358
 */
3359
static __isl_give isl_basic_map *drop_inequalities(
3360
  __isl_take isl_basic_map *bmap, __isl_keep isl_basic_map *context)
3361
329
{
3362
329
  int i1, i2;
3363
329
  unsigned total, extra;
3364
329
3365
329
  if (
!bmap || 329
!context329
)
3366
0
    return isl_basic_map_free(bmap);
3367
329
3368
329
  total = isl_basic_map_total_dim(context);
3369
329
  extra = isl_basic_map_total_dim(bmap) - total;
3370
329
3371
329
  i1 = bmap->n_ineq - 1;
3372
329
  i2 = context->n_ineq - 1;
3373
1.38k
  while (
bmap && 1.38k
i1 >= 01.38k
&&
i2 >= 01.17k
)
{1.05k
3374
1.05k
    int cmp;
3375
1.05k
3376
1.05k
    if (isl_seq_first_non_zero(bmap->ineq[i1] + 1 + total,
3377
18
              extra) != -1) {
3378
18
      --i1;
3379
18
      continue;
3380
18
    }
3381
1.03k
    cmp = isl_basic_map_constraint_cmp(context, bmap->ineq[i1],
3382
1.03k
              context->ineq[i2]);
3383
1.03k
    if (
cmp < 01.03k
)
{0
3384
0
      --i2;
3385
0
      continue;
3386
0
    }
3387
1.03k
    
if (1.03k
cmp > 01.03k
)
{289
3388
289
      --i1;
3389
289
      continue;
3390
289
    }
3391
745
    
if (745
isl_int_eq745
(bmap->ineq[i1][0], context->ineq[i2][0]))
{673
3392
673
      bmap = isl_basic_map_cow(bmap);
3393
673
      if (isl_basic_map_drop_inequality(bmap, i1) < 0)
3394
0
        bmap = isl_basic_map_free(bmap);
3395
673
    }
3396
745
    --i1;
3397
745
    --i2;
3398
745
  }
3399
329
3400
329
  return bmap;
3401
329
}
3402
3403
/* Drop all equalities from "bmap" that also appear in "context".
3404
 * "context" is assumed to have only known local variables and
3405
 * the initial local variables of "bmap" are assumed to be the same
3406
 * as those of "context".
3407
 *
3408
 * Run through the equality constraints of "bmap" and "context"
3409
 * in sorted order.
3410
 * If a constraint of "bmap" involves variables not in "context",
3411
 * then it cannot appear in "context".
3412
 * If a matching constraint is found, it is removed from "bmap".
3413
 */
3414
static __isl_give isl_basic_map *drop_equalities(
3415
  __isl_take isl_basic_map *bmap, __isl_keep isl_basic_map *context)
3416
329
{
3417
329
  int i1, i2;
3418
329
  unsigned total, extra;
3419
329
3420
329
  if (
!bmap || 329
!context329
)
3421
0
    return isl_basic_map_free(bmap);
3422
329
3423
329
  total = isl_basic_map_total_dim(context);
3424
329
  extra = isl_basic_map_total_dim(bmap) - total;
3425
329
3426
329
  i1 = bmap->n_eq - 1;
3427
329
  i2 = context->n_eq - 1;
3428
329
3429
355
  while (
bmap && 355
i1 >= 0355
&&
i2 >= 040
)
{26
3430
26
    int last1, last2;
3431
26
3432
26
    if (isl_seq_first_non_zero(bmap->eq[i1] + 1 + total,
3433
26
              extra) != -1)
3434
0
      break;
3435
26
    last1 = isl_seq_last_non_zero(bmap->eq[i1] + 1, total);
3436
26
    last2 = isl_seq_last_non_zero(context->eq[i2] + 1, total);
3437
26
    if (
last1 > last226
)
{0
3438
0
      --i2;
3439
0
      continue;
3440
0
    }
3441
26
    
if (26
last1 < last226
)
{3
3442
3
      --i1;
3443
3
      continue;
3444
3
    }
3445
23
    
if (23
isl_seq_eq(bmap->eq[i1], context->eq[i2], 1 + total)23
)
{23
3446
23
      bmap = isl_basic_map_cow(bmap);
3447
23
      if (isl_basic_map_drop_equality(bmap, i1) < 0)
3448
0
        bmap = isl_basic_map_free(bmap);
3449
23
    }
3450
23
    --i1;
3451
23
    --i2;
3452
23
  }
3453
329
3454
329
  return bmap;
3455
329
}
3456
3457
/* Remove the constraints in "context" from "bmap".
3458
 * "context" is assumed to have explicit representations
3459
 * for all local variables.
3460
 *
3461
 * First align the divs of "bmap" to those of "context" and
3462
 * sort the constraints.  Then drop all constraints from "bmap"
3463
 * that appear in "context".
3464
 */
3465
__isl_give isl_basic_map *isl_basic_map_plain_gist(
3466
  __isl_take isl_basic_map *bmap, __isl_take isl_basic_map *context)
3467
329
{
3468
329
  isl_bool done, known;
3469
329
3470
329
  done = isl_basic_map_plain_is_universe(context);
3471
329
  if (done == isl_bool_false)
3472
329
    done = isl_basic_map_plain_is_universe(bmap);
3473
329
  if (done == isl_bool_false)
3474
329
    done = isl_basic_map_plain_is_empty(context);
3475
329
  if (done == isl_bool_false)
3476
329
    done = isl_basic_map_plain_is_empty(bmap);
3477
329
  if (done < 0)
3478
0
    goto error;
3479
329
  
if (329
done329
)
{0
3480
0
    isl_basic_map_free(context);
3481
0
    return bmap;
3482
0
  }
3483
329
  known = isl_basic_map_divs_known(context);
3484
329
  if (known < 0)
3485
0
    goto error;
3486
329
  
if (329
!known329
)
3487
0
    isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
3488
329
      "context has unknown divs", goto error);
3489
329
3490
329
  bmap = isl_basic_map_align_divs(bmap, context);
3491
329
  bmap = isl_basic_map_gauss(bmap, NULL);
3492
329
  bmap = isl_basic_map_sort_constraints(bmap);
3493
329
  context = isl_basic_map_sort_constraints(context);
3494
329
3495
329
  bmap = drop_inequalities(bmap, context);
3496
329
  bmap = drop_equalities(bmap, context);
3497
329
3498
329
  isl_basic_map_free(context);
3499
329
  bmap = isl_basic_map_finalize(bmap);
3500
329
  return bmap;
3501
0
error:
3502
0
  isl_basic_map_free(bmap);
3503
0
  isl_basic_map_free(context);
3504
0
  return NULL;
3505
329
}
3506
3507
/* Replace "map" by the disjunct at position "pos" and free "context".
3508
 */
3509
static __isl_give isl_map *replace_by_disjunct(__isl_take isl_map *map,
3510
  int pos, __isl_take isl_basic_map *context)
3511
7
{
3512
7
  isl_basic_map *bmap;
3513
7
3514
7
  bmap = isl_basic_map_copy(map->p[pos]);
3515
7
  isl_map_free(map);
3516
7
  isl_basic_map_free(context);
3517
7
  return isl_map_from_basic_map(bmap);
3518
7
}
3519
3520
/* Remove the constraints in "context" from "map".
3521
 * If any of the disjuncts in the result turns out to be the universe,
3522
 * then return this universe.
3523
 * "context" is assumed to have explicit representations
3524
 * for all local variables.
3525
 */
3526
__isl_give isl_map *isl_map_plain_gist_basic_map(__isl_take isl_map *map,
3527
  __isl_take isl_basic_map *context)
3528
152
{
3529
152
  int i;
3530
152
  isl_bool univ, known;
3531
152
3532
152
  univ = isl_basic_map_plain_is_universe(context);
3533
152
  if (univ < 0)
3534
0
    goto error;
3535
152
  
if (152
univ152
)
{0
3536
0
    isl_basic_map_free(context);
3537
0
    return map;
3538
0
  }
3539
152
  known = isl_basic_map_divs_known(context);
3540
152
  if (known < 0)
3541
0
    goto error;
3542
152
  
if (152
!known152
)
3543
0
    isl_die(isl_map_get_ctx(map), isl_error_invalid,
3544
152
      "context has unknown divs", goto error);
3545
152
3546
152
  map = isl_map_cow(map);
3547
152
  if (!map)
3548
0
    goto error;
3549
474
  
for (i = 0; 152
i < map->n474
;
++i322
)
{329
3550
329
    map->p[i] = isl_basic_map_plain_gist(map->p[i],
3551
329
            isl_basic_map_copy(context));
3552
329
    univ = isl_basic_map_plain_is_universe(map->p[i]);
3553
329
    if (univ < 0)
3554
0
      goto error;
3555
329
    
if (329
univ && 329
map->n > 17
)
3556
7
      return replace_by_disjunct(map, i, context);
3557
329
  }
3558
152
3559
145
  isl_basic_map_free(context);
3560
145
  ISL_F_CLR(map, ISL_MAP_NORMALIZED);
3561
145
  if (map->n > 1)
3562
145
    ISL_F_CLR(map, ISL_MAP_DISJOINT);
3563
145
  return map;
3564
0
error:
3565
0
  isl_map_free(map);
3566
0
  isl_basic_map_free(context);
3567
0
  return NULL;
3568
152
}
3569
3570
/* Replace "map" by a universe map in the same space and free "drop".
3571
 */
3572
static __isl_give isl_map *replace_by_universe(__isl_take isl_map *map,
3573
  __isl_take isl_map *drop)
3574
4.77k
{
3575
4.77k
  isl_map *res;
3576
4.77k
3577
4.77k
  res = isl_map_universe(isl_map_get_space(map));
3578
4.77k
  isl_map_free(map);
3579
4.77k
  isl_map_free(drop);
3580
4.77k
  return res;
3581
4.77k
}
3582
3583
/* Return a map that has the same intersection with "context" as "map"
3584
 * and that is as "simple" as possible.
3585
 *
3586
 * If "map" is already the universe, then we cannot make it any simpler.
3587
 * Similarly, if "context" is the universe, then we cannot exploit it
3588
 * to simplify "map"
3589
 * If "map" and "context" are identical to each other, then we can
3590
 * return the corresponding universe.
3591
 *
3592
 * If either "map" or "context" consists of multiple disjuncts,
3593
 * then check if "context" happens to be a subset of "map",
3594
 * in which case all constraints can be removed.
3595
 * In case of multiple disjuncts, the standard procedure
3596
 * may not be able to detect that all constraints can be removed.
3597
 *
3598
 * If none of these cases apply, we have to work a bit harder.
3599
 * During this computation, we make use of a single disjunct context,
3600
 * so if the original context consists of more than one disjunct
3601
 * then we need to approximate the context by a single disjunct set.
3602
 * Simply taking the simple hull may drop constraints that are
3603
 * only implicitly available in each disjunct.  We therefore also
3604
 * look for constraints among those defining "map" that are valid
3605
 * for the context.  These can then be used to simplify away
3606
 * the corresponding constraints in "map".
3607
 */
3608
static __isl_give isl_map *map_gist(__isl_take isl_map *map,
3609
  __isl_take isl_map *context)
3610
38.2k
{
3611
38.2k
  int equal;
3612
38.2k
  int is_universe;
3613
38.2k
  int single_disjunct_map, single_disjunct_context;
3614
38.2k
  isl_bool subset;
3615
38.2k
  isl_basic_map *hull;
3616
38.2k
3617
38.2k
  is_universe = isl_map_plain_is_universe(map);
3618
38.2k
  if (
is_universe >= 0 && 38.2k
!is_universe38.2k
)
3619
23.3k
    is_universe = isl_map_plain_is_universe(context);
3620
38.2k
  if (is_universe < 0)
3621
0
    goto error;
3622
38.2k
  
if (38.2k
is_universe38.2k
)
{21.6k
3623
21.6k
    isl_map_free(context);
3624
21.6k
    return map;
3625
21.6k
  }
3626
38.2k
3627
16.5k
  equal = isl_map_plain_is_equal(map, context);
3628
16.5k
  if (equal < 0)
3629
0
    goto error;
3630
16.5k
  
if (16.5k
equal16.5k
)
3631
4.74k
    return replace_by_universe(map, context);
3632
16.5k
3633
11.8k
  single_disjunct_map = isl_map_n_basic_map(map) == 1;
3634
11.8k
  single_disjunct_context = isl_map_n_basic_map(context) == 1;
3635
11.8k
  if (
!single_disjunct_map || 11.8k
!single_disjunct_context7.87k
)
{4.19k
3636
4.19k
    subset = isl_map_is_subset(context, map);
3637
4.19k
    if (subset < 0)
3638
0
      goto error;
3639
4.19k
    
if (4.19k
subset4.19k
)
3640
24
      return replace_by_universe(map, context);
3641
4.19k
  }
3642
11.8k
3643
11.8k
  context = isl_map_compute_divs(context);
3644
11.8k
  if (!context)
3645
0
    goto error;
3646
11.8k
  
if (11.8k
single_disjunct_context11.8k
)
{11.4k
3647
11.4k
    hull = isl_map_simple_hull(context);
3648
397
  } else {
3649
397
    isl_ctx *ctx;
3650
397
    isl_map_list *list;
3651
397
3652
397
    ctx = isl_map_get_ctx(map);
3653
397
    list = isl_map_list_alloc(ctx, 2);
3654
397
    list = isl_map_list_add(list, isl_map_copy(context));
3655
397
    list = isl_map_list_add(list, isl_map_copy(map));
3656
397
    hull = isl_map_unshifted_simple_hull_from_map_list(context,
3657
397
                    list);
3658
397
  }
3659
11.8k
  return isl_map_gist_basic_map(map, hull);
3660
0
error:
3661
0
  isl_map_free(map);
3662
0
  isl_map_free(context);
3663
0
  return NULL;
3664
11.8k
}
3665
3666
__isl_give isl_map *isl_map_gist(__isl_take isl_map *map,
3667
  __isl_take isl_map *context)
3668
38.2k
{
3669
38.2k
  return isl_map_align_params_map_map_and(map, context, &map_gist);
3670
38.2k
}
3671
3672
struct isl_basic_set *isl_basic_set_gist(struct isl_basic_set *bset,
3673
            struct isl_basic_set *context)
3674
719
{
3675
719
  return bset_from_bmap(isl_basic_map_gist(bset_to_bmap(bset),
3676
719
            bset_to_bmap(context)));
3677
719
}
3678
3679
__isl_give isl_set *isl_set_gist_basic_set(__isl_take isl_set *set,
3680
  __isl_take isl_basic_set *context)
3681
1.76k
{
3682
1.76k
  return set_from_map(isl_map_gist_basic_map(set_to_map(set),
3683
1.76k
          bset_to_bmap(context)));
3684
1.76k
}
3685
3686
__isl_give isl_set *isl_set_gist_params_basic_set(__isl_take isl_set *set,
3687
  __isl_take isl_basic_set *context)
3688
0
{
3689
0
  isl_space *space = isl_set_get_space(set);
3690
0
  isl_basic_set *dom_context = isl_basic_set_universe(space);
3691
0
  dom_context = isl_basic_set_intersect_params(dom_context, context);
3692
0
  return isl_set_gist_basic_set(set, dom_context);
3693
0
}
3694
3695
__isl_give isl_set *isl_set_gist(__isl_take isl_set *set,
3696
  __isl_take isl_set *context)
3697
10.3k
{
3698
10.3k
  return set_from_map(isl_map_gist(set_to_map(set), set_to_map(context)));
3699
10.3k
}
3700
3701
/* Compute the gist of "bmap" with respect to the constraints "context"
3702
 * on the domain.
3703
 */
3704
__isl_give isl_basic_map *isl_basic_map_gist_domain(
3705
  __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *context)
3706
0
{
3707
0
  isl_space *space = isl_basic_map_get_space(bmap);
3708
0
  isl_basic_map *bmap_context = isl_basic_map_universe(space);
3709
0
3710
0
  bmap_context = isl_basic_map_intersect_domain(bmap_context, context);
3711
0
  return isl_basic_map_gist(bmap, bmap_context);
3712
0
}
3713
3714
__isl_give isl_map *isl_map_gist_domain(__isl_take isl_map *map,
3715
  __isl_take isl_set *context)
3716
4.55k
{
3717
4.55k
  isl_map *map_context = isl_map_universe(isl_map_get_space(map));
3718
4.55k
  map_context = isl_map_intersect_domain(map_context, context);
3719
4.55k
  return isl_map_gist(map, map_context);
3720
4.55k
}
3721
3722
__isl_give isl_map *isl_map_gist_range(__isl_take isl_map *map,
3723
  __isl_take isl_set *context)
3724
74
{
3725
74
  isl_map *map_context = isl_map_universe(isl_map_get_space(map));
3726
74
  map_context = isl_map_intersect_range(map_context, context);
3727
74
  return isl_map_gist(map, map_context);
3728
74
}
3729
3730
__isl_give isl_map *isl_map_gist_params(__isl_take isl_map *map,
3731
  __isl_take isl_set *context)
3732
20.1k
{
3733
20.1k
  isl_map *map_context = isl_map_universe(isl_map_get_space(map));
3734
20.1k
  map_context = isl_map_intersect_params(map_context, context);
3735
20.1k
  return isl_map_gist(map, map_context);
3736
20.1k
}
3737
3738
__isl_give isl_set *isl_set_gist_params(__isl_take isl_set *set,
3739
  __isl_take isl_set *context)
3740
15.1k
{
3741
15.1k
  return isl_map_gist_params(set, context);
3742
15.1k
}
3743
3744
/* Quick check to see if two basic maps are disjoint.
3745
 * In particular, we reduce the equalities and inequalities of
3746
 * one basic map in the context of the equalities of the other
3747
 * basic map and check if we get a contradiction.
3748
 */
3749
isl_bool isl_basic_map_plain_is_disjoint(__isl_keep isl_basic_map *bmap1,
3750
  __isl_keep isl_basic_map *bmap2)
3751
15.1k
{
3752
15.1k
  struct isl_vec *v = NULL;
3753
15.1k
  int *elim = NULL;
3754
15.1k
  unsigned total;
3755
15.1k
  int i;
3756
15.1k
3757
15.1k
  if (
!bmap1 || 15.1k
!bmap215.1k
)
3758
0
    return isl_bool_error;
3759
15.1k
  
isl_assert15.1k
(bmap1->ctx, isl_space_is_equal(bmap1->dim, bmap2->dim),15.1k
3760
15.1k
      return isl_bool_error);
3761
15.1k
  
if (15.1k
bmap1->n_div || 15.1k
bmap2->n_div14.3k
)
3762
944
    return isl_bool_false;
3763
14.2k
  
if (14.2k
!bmap1->n_eq && 14.2k
!bmap2->n_eq6.84k
)
3764
6.08k
    return isl_bool_false;
3765
14.2k
3766
8.13k
  total = isl_space_dim(bmap1->dim, isl_dim_all);
3767
8.13k
  if (total == 0)
3768
0
    return isl_bool_false;
3769
8.13k
  v = isl_vec_alloc(bmap1->ctx, 1 + total);
3770
8.13k
  if (!v)
3771
0
    goto error;
3772
8.13k
  
elim = 8.13k
isl_alloc_array8.13k
(bmap1->ctx, int, total);
3773
8.13k
  if (!elim)
3774
0
    goto error;
3775
8.13k
  compute_elimination_index(bmap1, elim);
3776
16.9k
  for (i = 0; 
i < bmap2->n_eq16.9k
;
++i8.84k
)
{9.02k
3777
9.02k
    int reduced;
3778
9.02k
    reduced = reduced_using_equalities(v->block.data, bmap2->eq[i],
3779
9.02k
              bmap1, elim);
3780
9.02k
    if (
reduced && 9.02k
!3.59k
isl_int_is_zero3.59k
(v->block.data[0]) &&
3781
457
        isl_seq_first_non_zero(v->block.data + 1, total) == -1)
3782
172
      goto disjoint;
3783
9.02k
  }
3784
74.1k
  
for (i = 0; 7.96k
i < bmap2->n_ineq74.1k
;
++i66.2k
)
{69.0k
3785
69.0k
    int reduced;
3786
69.0k
    reduced = reduced_using_equalities(v->block.data,
3787
69.0k
            bmap2->ineq[i], bmap1, elim);
3788
69.0k
    if (
reduced && 69.0k
isl_int_is_neg3.76k
(v->block.data[0]) &&
3789
3.04k
        isl_seq_first_non_zero(v->block.data + 1, total) == -1)
3790
2.85k
      goto disjoint;
3791
69.0k
  }
3792
5.10k
  compute_elimination_index(bmap2, elim);
3793
41.3k
  for (i = 0; 
i < bmap1->n_ineq41.3k
;
++i36.2k
)
{37.6k
3794
37.6k
    int reduced;
3795
37.6k
    reduced = reduced_using_equalities(v->block.data,
3796
37.6k
            bmap1->ineq[i], bmap2, elim);
3797
37.6k
    if (
reduced && 37.6k
isl_int_is_neg2.09k
(v->block.data[0]) &&
3798
1.53k
        isl_seq_first_non_zero(v->block.data + 1, total) == -1)
3799
1.42k
      goto disjoint;
3800
37.6k
  }
3801
3.68k
  isl_vec_free(v);
3802
3.68k
  free(elim);
3803
3.68k
  return isl_bool_false;
3804
4.45k
disjoint:
3805
4.45k
  isl_vec_free(v);
3806
4.45k
  free(elim);
3807
4.45k
  return isl_bool_true;
3808
0
error:
3809
0
  isl_vec_free(v);
3810
0
  free(elim);
3811
0
  return isl_bool_error;
3812
5.10k
}
3813
3814
int isl_basic_set_plain_is_disjoint(__isl_keep isl_basic_set *bset1,
3815
  __isl_keep isl_basic_set *bset2)
3816
0
{
3817
0
  return isl_basic_map_plain_is_disjoint(bset_to_bmap(bset1),
3818
0
                bset_to_bmap(bset2));
3819
0
}
3820
3821
/* Does "test" hold for all pairs of basic maps in "map1" and "map2"?
3822
 */
3823
static isl_bool all_pairs(__isl_keep isl_map *map1, __isl_keep isl_map *map2,
3824
  isl_bool (*test)(__isl_keep isl_basic_map *bmap1,
3825
    __isl_keep isl_basic_map *bmap2))
3826
11.8k
{
3827
11.8k
  int i, j;
3828
11.8k
3829
11.8k
  if (
!map1 || 11.8k
!map211.8k
)
3830
0
    return isl_bool_error;
3831
11.8k
3832
17.8k
  
for (i = 0; 11.8k
i < map1->n17.8k
;
++i5.95k
)
{13.9k
3833
21.1k
    for (j = 0; 
j < map2->n21.1k
;
++j7.18k
)
{15.1k
3834
15.1k
      isl_bool d = test(map1->p[i], map2->p[j]);
3835
15.1k
      if (d != isl_bool_true)
3836
7.96k
        return d;
3837
15.1k
    }
3838
13.9k
  }
3839
11.8k
3840
3.93k
  return isl_bool_true;
3841
11.8k
}
3842
3843
/* Are "map1" and "map2" obviously disjoint, based on information
3844
 * that can be derived without looking at the individual basic maps?
3845
 *
3846
 * In particular, if one of them is empty or if they live in different spaces
3847
 * (ignoring parameters), then they are clearly disjoint.
3848
 */
3849
static isl_bool isl_map_plain_is_disjoint_global(__isl_keep isl_map *map1,
3850
  __isl_keep isl_map *map2)
3851
24.2k
{
3852
24.2k
  isl_bool disjoint;
3853
24.2k
  isl_bool match;
3854
24.2k
3855
24.2k
  if (
!map1 || 24.2k
!map224.2k
)
3856
0
    return isl_bool_error;
3857
24.2k
3858
24.2k
  disjoint = isl_map_plain_is_empty(map1);
3859
24.2k
  if (
disjoint < 0 || 24.2k
disjoint24.2k
)
3860
3.38k
    return disjoint;
3861
24.2k
3862
20.8k
  disjoint = isl_map_plain_is_empty(map2);
3863
20.8k
  if (
disjoint < 0 || 20.8k
disjoint20.8k
)
3864
4.77k
    return disjoint;
3865
20.8k
3866
16.0k
  match = isl_space_tuple_is_equal(map1->dim, isl_dim_in,
3867
16.0k
        map2->dim, isl_dim_in);
3868
16.0k
  if (
match < 0 || 16.0k
!match16.0k
)
3869
0
    
return match < 0 ? 0
isl_bool_error0
:
isl_bool_true0
;
3870
16.0k
3871
16.0k
  match = isl_space_tuple_is_equal(map1->dim, isl_dim_out,
3872
16.0k
        map2->dim, isl_dim_out);
3873
16.0k
  if (
match < 0 || 16.0k
!match16.0k
)
3874
0
    
return match < 0 ? 0
isl_bool_error0
:
isl_bool_true0
;
3875
16.0k
3876
16.0k
  return isl_bool_false;
3877
16.0k
}
3878
3879
/* Are "map1" and "map2" obviously disjoint?
3880
 *
3881
 * If one of them is empty or if they live in different spaces (ignoring
3882
 * parameters), then they are clearly disjoint.
3883
 * This is checked by isl_map_plain_is_disjoint_global.
3884
 *
3885
 * If they have different parameters, then we skip any further tests.
3886
 *
3887
 * If they are obviously equal, but not obviously empty, then we will
3888
 * not be able to detect if they are disjoint.
3889
 *
3890
 * Otherwise we check if each basic map in "map1" is obviously disjoint
3891
 * from each basic map in "map2".
3892
 */
3893
isl_bool isl_map_plain_is_disjoint(__isl_keep isl_map *map1,
3894
  __isl_keep isl_map *map2)
3895
0
{
3896
0
  isl_bool disjoint;
3897
0
  isl_bool intersect;
3898
0
  isl_bool match;
3899
0
3900
0
  disjoint = isl_map_plain_is_disjoint_global(map1, map2);
3901
0
  if (
disjoint < 0 || 0
disjoint0
)
3902
0
    return disjoint;
3903
0
3904
0
  match = isl_map_has_equal_params(map1, map2);
3905
0
  if (
match < 0 || 0
!match0
)
3906
0
    
return match < 0 ? 0
isl_bool_error0
:
isl_bool_false0
;
3907
0
3908
0
  intersect = isl_map_plain_is_equal(map1, map2);
3909
0
  if (
intersect < 0 || 0
intersect0
)
3910
0
    
return intersect < 0 ? 0
isl_bool_error0
:
isl_bool_false0
;
3911
0
3912
0
  return all_pairs(map1, map2, &isl_basic_map_plain_is_disjoint);
3913
0
}
3914
3915
/* Are "map1" and "map2" disjoint?
3916
 *
3917
 * They are disjoint if they are "obviously disjoint" or if one of them
3918
 * is empty.  Otherwise, they are not disjoint if one of them is universal.
3919
 * If the two inputs are (obviously) equal and not empty, then they are
3920
 * not disjoint.
3921
 * If none of these cases apply, then check if all pairs of basic maps
3922
 * are disjoint.
3923
 */
3924
isl_bool isl_map_is_disjoint(__isl_keep isl_map *map1, __isl_keep isl_map *map2)
3925
24.2k
{
3926
24.2k
  isl_bool disjoint;
3927
24.2k
  isl_bool intersect;
3928
24.2k
3929
24.2k
  disjoint = isl_map_plain_is_disjoint_global(map1, map2);
3930
24.2k
  if (
disjoint < 0 || 24.2k
disjoint24.2k
)
3931
8.15k
    return disjoint;
3932
24.2k
3933
16.0k
  disjoint = isl_map_is_empty(map1);
3934
16.0k
  if (
disjoint < 0 || 16.0k
disjoint16.0k
)
3935
0
    return disjoint;
3936
16.0k
3937
16.0k
  disjoint = isl_map_is_empty(map2);
3938
16.0k
  if (
disjoint < 0 || 16.0k
disjoint16.0k
)
3939
0
    return disjoint;
3940
16.0k
3941
16.0k
  intersect = isl_map_plain_is_universe(map1);
3942
16.0k
  if (
intersect < 0 || 16.0k
intersect16.0k
)
3943
2.66k
    
return intersect < 0 ? 2.66k
isl_bool_error0
:
isl_bool_false2.66k
;
3944
16.0k
3945
13.4k
  intersect = isl_map_plain_is_universe(map2);
3946
13.4k
  if (
intersect < 0 || 13.4k
intersect13.4k
)
3947
1.12k
    
return intersect < 0 ? 1.12k
isl_bool_error0
:
isl_bool_false1.12k
;
3948
13.4k
3949
12.3k
  intersect = isl_map_plain_is_equal(map1, map2);
3950
12.3k
  if (
intersect < 0 || 12.3k
intersect12.3k
)
3951
416
    return isl_bool_not(intersect);
3952
12.3k
3953
11.8k
  return all_pairs(map1, map2, &isl_basic_map_is_disjoint);
3954
12.3k
}
3955
3956
/* Are "bmap1" and "bmap2" disjoint?
3957
 *
3958
 * They are disjoint if they are "obviously disjoint" or if one of them
3959
 * is empty.  Otherwise, they are not disjoint if one of them is universal.
3960
 * If none of these cases apply, we compute the intersection and see if
3961
 * the result is empty.
3962
 */
3963
isl_bool isl_basic_map_is_disjoint(__isl_keep isl_basic_map *bmap1,
3964
  __isl_keep isl_basic_map *bmap2)
3965
15.1k
{
3966
15.1k
  isl_bool disjoint;
3967
15.1k
  isl_bool intersect;
3968
15.1k
  isl_basic_map *test;
3969
15.1k
3970
15.1k
  disjoint = isl_basic_map_plain_is_disjoint(bmap1, bmap2);
3971
15.1k
  if (
disjoint < 0 || 15.1k
disjoint15.1k
)
3972
4.45k
    return disjoint;
3973
15.1k
3974
10.7k
  disjoint = isl_basic_map_is_empty(bmap1);
3975
10.7k
  if (
disjoint < 0 || 10.7k
disjoint10.7k
)
3976
0
    return disjoint;
3977
10.7k
3978
10.7k
  disjoint = isl_basic_map_is_empty(bmap2);
3979
10.7k
  if (
disjoint < 0 || 10.7k
disjoint10.7k
)
3980
0
    return disjoint;
3981
10.7k
3982
10.7k
  intersect = isl_basic_map_plain_is_universe(bmap1);
3983
10.7k
  if (
intersect < 0 || 10.7k
intersect10.7k
)
3984
0
    
return intersect < 0 ? 0
isl_bool_error0
:
isl_bool_false0
;
3985
10.7k
3986
10.7k
  intersect = isl_basic_map_plain_is_universe(bmap2);
3987
10.7k
  if (
intersect < 0 || 10.7k
intersect10.7k
)
3988
0
    
return intersect < 0 ? 0
isl_bool_error0
:
isl_bool_false0
;
3989
10.7k
3990
10.7k
  test = isl_basic_map_intersect(isl_basic_map_copy(bmap1),
3991
10.7k
    isl_basic_map_copy(bmap2));
3992
10.7k
  disjoint = isl_basic_map_is_empty(test);
3993
10.7k
  isl_basic_map_free(test);
3994
10.7k
3995
10.7k
  return disjoint;
3996
10.7k
}
3997
3998
/* Are "bset1" and "bset2" disjoint?
3999
 */
4000
isl_bool isl_basic_set_is_disjoint(__isl_keep isl_basic_set *bset1,
4001
  __isl_keep isl_basic_set *bset2)
4002
21
{
4003
21
  return isl_basic_map_is_disjoint(bset1, bset2);
4004
21
}
4005
4006
isl_bool isl_set_plain_is_disjoint(__isl_keep isl_set *set1,
4007
  __isl_keep isl_set *set2)
4008
0
{
4009
0
  return isl_map_plain_is_disjoint(set_to_map(set1), set_to_map(set2));
4010
0
}
4011
4012
/* Are "set1" and "set2" disjoint?
4013
 */
4014
isl_bool isl_set_is_disjoint(__isl_keep isl_set *set1, __isl_keep isl_set *set2)
4015
8.47k
{
4016
8.47k
  return isl_map_is_disjoint(set1, set2);
4017
8.47k
}
4018
4019
/* Is "v" equal to 0, 1 or -1?
4020
 */
4021
static int is_zero_or_one(isl_int v)
4022
118
{
4023
118
  return 
isl_int_is_zero118
(v) ||
isl_int_is_one51
(v) ||
isl_int_is_negone29
(v);
4024
118
}
4025
4026
/* Check if we can combine a given div with lower bound l and upper
4027
 * bound u with some other div and if so return that other div.
4028
 * Otherwise return -1.
4029
 *
4030
 * We first check that
4031
 *  - the bounds are opposites of each other (except for the constant
4032
 *    term)
4033
 *  - the bounds do not reference any other div
4034
 *  - no div is defined in terms of this div
4035
 *
4036
 * Let m be the size of the range allowed on the div by the bounds.
4037
 * That is, the bounds are of the form
4038
 *
4039
 *  e <= a <= e + m - 1
4040
 *
4041
 * with e some expression in the other variables.
4042
 * We look for another div b such that no third div is defined in terms
4043
 * of this second div b and such that in any constraint that contains
4044
 * a (except for the given lower and upper bound), also contains b
4045
 * with a coefficient that is m times that of b.
4046
 * That is, all constraints (execpt for the lower and upper bound)
4047
 * are of the form
4048
 *
4049
 *  e + f (a + m b) >= 0
4050
 *
4051
 * Furthermore, in the constraints that only contain b, the coefficient
4052
 * of b should be equal to 1 or -1.
4053
 * If so, we return b so that "a + m b" can be replaced by
4054
 * a single div "c = a + m b".
4055
 */
4056
static int div_find_coalesce(struct isl_basic_map *bmap, int *pairs,
4057
  unsigned div, unsigned l, unsigned u)
4058
474
{
4059
474
  int i, j;
4060
474
  unsigned dim;
4061
474
  int coalesce = -1;
4062
474
4063
474
  if (bmap->n_div <= 1)
4064
262
    return -1;
4065
212
  dim = isl_space_dim(bmap->dim, isl_dim_all);
4066
212
  if (isl_seq_first_non_zero(bmap->ineq[l] + 1 + dim, div) != -1)
4067
6
    return -1;
4068
206
  
if (206
isl_seq_first_non_zero(bmap->ineq[l] + 1 + dim + div + 1,206
4069
206
           bmap->n_div - div - 1) != -1)
4070
19
    return -1;
4071
187
  
if (187
!isl_seq_is_neg(bmap->ineq[l] + 1, bmap->ineq[u] + 1,187
4072
187
          dim + bmap->n_div))
4073
119
    return -1;
4074
187
4075
224
  
for (i = 0; 68
i < bmap->n_div224
;
++i156
)
{156
4076
156
    if (isl_int_is_zero(bmap->div[i][0]))
4077
102
      continue;
4078
54
    
if (54
!54
isl_int_is_zero54
(bmap->div[i][1 + 1 + dim + div]))
4079
0
      return -1;
4080
54
  }
4081
68
4082
68
  
isl_int_add68
(bmap->ineq[l][0], bmap->ineq[l][0], bmap->ineq[u][0]);68
4083
68
  if (
isl_int_is_neg68
(bmap->ineq[l][0]))
{0
4084
0
    isl_int_sub(bmap->ineq[l][0],
4085
0
          bmap->ineq[l][0], bmap->ineq[u][0]);
4086
0
    bmap = isl_basic_map_copy(bmap);
4087
0
    bmap = isl_basic_map_set_to_empty(bmap);
4088
0
    isl_basic_map_free(bmap);
4089
0
    return -1;
4090
0
  }
4091
68
  
isl_int_add_ui68
(bmap->ineq[l][0], bmap->ineq[l][0], 1);68
4092
222
  for (i = 0; 
i < bmap->n_div222
;
++i154
)
{155
4093
155
    if (i == div)
4094
67
      continue;
4095
88
    
if (88
!pairs[i]88
)
4096
48
      continue;
4097
150
    
for (j = 0; 40
j < bmap->n_div150
;
++j110
)
{110
4098
110
      if (isl_int_is_zero(bmap->div[j][0]))
4099
79
        continue;
4100
31
      
if (31
!31
isl_int_is_zero31
(bmap->div[j][1 + 1 + dim + i]))
4101
0
        break;
4102
31
    }
4103
40
    if (j < bmap->n_div)
4104
0
      continue;
4105
164
    
for (j = 0; 40
j < bmap->n_ineq164
;
++j124
)
{163
4106
163
      int valid;
4107
163
      if (
j == l || 163
j == u152
)
4108
22
        continue;
4109
141
      
if (141
isl_int_is_zero141
(bmap->ineq[j][1 + dim + div]))
{118
4110
118
        if (is_zero_or_one(bmap->ineq[j][1 + dim + i]))
4111
98
          continue;
4112
20
        break;
4113
118
      }
4114
23
      
if (23
isl_int_is_zero23
(bmap->ineq[j][1 + dim + i]))
4115
11
        break;
4116
12
      
isl_int_mul12
(bmap->ineq[j][1 + dim + div],12
4117
12
            bmap->ineq[j][1 + dim + div],
4118
12
            bmap->ineq[l][0]);
4119
12
      valid = isl_int_eq(bmap->ineq[j][1 + dim + div],
4120
12
             bmap->ineq[j][1 + dim + i]);
4121
12
      isl_int_divexact(bmap->ineq[j][1 + dim + div],
4122
12
           bmap->ineq[j][1 + dim + div],
4123
12
           bmap->ineq[l][0]);
4124
12
      if (!valid)
4125
8
        break;
4126
12
    }
4127
40
    if (j < bmap->n_ineq)
4128
39
      continue;
4129
1
    coalesce = i;
4130
1
    break;
4131
40
  }
4132
68
  isl_int_sub_ui(bmap->ineq[l][0], bmap->ineq[l][0], 1);
4133
68
  isl_int_sub(bmap->ineq[l][0], bmap->ineq[l][0], bmap->ineq[u][0]);
4134
68
  return coalesce;
4135
68
}
4136
4137
/* Internal data structure used during the construction and/or evaluation of
4138
 * an inequality that ensures that a pair of bounds always allows
4139
 * for an integer value.
4140
 *
4141
 * "tab" is the tableau in which the inequality is evaluated.  It may
4142
 * be NULL until it is actually needed.
4143
 * "v" contains the inequality coefficients.
4144
 * "g", "fl" and "fu" are temporary scalars used during the construction and
4145
 * evaluation.
4146
 */
4147
struct test_ineq_data {
4148
  struct isl_tab *tab;
4149
  isl_vec *v;
4150
  isl_int g;
4151
  isl_int fl;
4152
  isl_int fu;
4153
};
4154
4155
/* Free all the memory allocated by the fields of "data".
4156
 */
4157
static void test_ineq_data_clear(struct test_ineq_data *data)
4158
708
{
4159
708
  isl_tab_free(data->tab);
4160
708
  isl_vec_free(data->v);
4161
708
  isl_int_clear(data->g);
4162
708
  isl_int_clear(data->fl);
4163
708
  isl_int_clear(data->fu);
4164
708
}
4165
4166
/* Is the inequality stored in data->v satisfied by "bmap"?
4167
 * That is, does it only attain non-negative values?
4168
 * data->tab is a tableau corresponding to "bmap".
4169
 */
4170
static isl_bool test_ineq_is_satisfied(__isl_keep isl_basic_map *bmap,
4171
  struct test_ineq_data *data)
4172
515
{
4173
515
  isl_ctx *ctx;
4174
515
  enum isl_lp_result res;
4175
515
4176
515
  ctx = isl_basic_map_get_ctx(bmap);
4177
515
  if (!data->tab)
4178
411
    data->tab = isl_tab_from_basic_map(bmap, 0);
4179
515
  res = isl_tab_min(data->tab, data->v->el, ctx->one, &data->g, NULL, 0);
4180
515
  if (res == isl_lp_error)
4181
0
    return isl_bool_error;
4182
515
  
return res == isl_lp_ok && 515
isl_int_is_nonneg515
(data->g);
4183
515
}
4184
4185
/* Given a lower and an upper bound on div i, do they always allow
4186
 * for an integer value of the given div?
4187
 * Determine this property by constructing an inequality
4188
 * such that the property is guaranteed when the inequality is nonnegative.
4189
 * The lower bound is inequality l, while the upper bound is inequality u.
4190
 * The constructed inequality is stored in data->v.
4191
 *
4192
 * Let the upper bound be
4193
 *
4194
 *  -n_u a + e_u >= 0
4195
 *
4196
 * and the lower bound
4197
 *
4198
 *  n_l a + e_l >= 0
4199
 *
4200
 * Let n_u = f_u g and n_l = f_l g, with g = gcd(n_u, n_l).
4201
 * We have
4202
 *
4203
 *  - f_u e_l <= f_u f_l g a <= f_l e_u
4204
 *
4205
 * Since all variables are integer valued, this is equivalent to
4206
 *
4207
 *  - f_u e_l - (f_u - 1) <= f_u f_l g a <= f_l e_u + (f_l - 1)
4208
 *
4209
 * If this interval is at least f_u f_l g, then it contains at least
4210
 * one integer value for a.
4211
 * That is, the test constraint is
4212
 *
4213
 *  f_l e_u + f_u e_l + f_l - 1 + f_u - 1 + 1 >= f_u f_l g
4214
 *
4215
 * or
4216
 *
4217
 *  f_l e_u + f_u e_l + f_l - 1 + f_u - 1 + 1 - f_u f_l g >= 0
4218
 *
4219
 * If the coefficients of f_l e_u + f_u e_l have a common divisor g',
4220
 * then the constraint can be scaled down by a factor g',
4221
 * with the constant term replaced by
4222
 * floor((f_l e_{u,0} + f_u e_{l,0} + f_l - 1 + f_u - 1 + 1 - f_u f_l g)/g').
4223
 * Note that the result of applying Fourier-Motzkin to this pair
4224
 * of constraints is
4225
 *
4226
 *  f_l e_u + f_u e_l >= 0
4227
 *
4228
 * If the constant term of the scaled down version of this constraint,
4229
 * i.e., floor((f_l e_{u,0} + f_u e_{l,0})/g') is equal to the constant
4230
 * term of the scaled down test constraint, then the test constraint
4231
 * is known to hold and no explicit evaluation is required.
4232
 * This is essentially the Omega test.
4233
 *
4234
 * If the test constraint consists of only a constant term, then
4235
 * it is sufficient to look at the sign of this constant term.
4236
 */
4237
static isl_bool int_between_bounds(__isl_keep isl_basic_map *bmap, int i,
4238
  int l, int u, struct test_ineq_data *data)
4239
1.44k
{
4240
1.44k
  unsigned offset, n_div;
4241
1.44k
  offset = isl_basic_map_offset(bmap, isl_dim_div);
4242
1.44k
  n_div = isl_basic_map_dim(bmap, isl_dim_div);
4243
1.44k
4244
1.44k
  isl_int_gcd(data->g,
4245
1.44k
        bmap->ineq[l][offset + i], bmap->ineq[u][offset + i]);
4246
1.44k
  isl_int_divexact(data->fl, bmap->ineq[l][offset + i], data->g);
4247
1.44k
  isl_int_divexact(data->fu, bmap->ineq[u][offset + i], data->g);
4248
1.44k
  isl_int_neg(data->fu, data->fu);
4249
1.44k
  isl_seq_combine(data->v->el, data->fl, bmap->ineq[u],
4250
1.44k
      data->fu, bmap->ineq[l], offset + n_div);
4251
1.44k
  isl_int_mul(data->g, data->g, data->fl);
4252
1.44k
  isl_int_mul(data->g, data->g, data->fu);
4253
1.44k
  isl_int_sub(data->g, data->g, data->fl);
4254
1.44k
  isl_int_sub(data->g, data->g, data->fu);
4255
1.44k
  isl_int_add_ui(data->g, data->g, 1);
4256
1.44k
  isl_int_sub(data->fl, data->v->el[0], data->g);
4257
1.44k
4258
1.44k
  isl_seq_gcd(data->v->el + 1, offset - 1 + n_div, &data->g);
4259
1.44k
  if (isl_int_is_zero(data->g))
4260
423
    
return 423
isl_int_is_nonneg423
(data->fl);
4261
1.01k
  
if (1.01k
isl_int_is_one1.01k
(data->g))
{349
4262
349
    isl_int_set(data->v->el[0], data->fl);
4263
349
    return test_ineq_is_satisfied(bmap, data);
4264
349
  }
4265
670
  
isl_int_fdiv_q670
(data->fl, data->fl, data->g);670
4266
670
  isl_int_fdiv_q(data->v->el[0], data->v->el[0], data->g);
4267
670
  if (isl_int_eq(data->fl, data->v->el[0]))
4268
504
    return isl_bool_true;
4269
166
  
isl_int_set166
(data->v->el[0], data->fl);166
4270
166
  isl_seq_scale_down(data->v->el + 1, data->v->el + 1, data->g,
4271
166
          offset - 1 + n_div);
4272
166
4273
166
  return test_ineq_is_satisfied(bmap, data);
4274
670
}
4275
4276
/* Remove more kinds of divs that are not strictly needed.
4277
 * In particular, if all pairs of lower and upper bounds on a div
4278
 * are such that they allow at least one integer value of the div,
4279
 * then we can eliminate the div using Fourier-Motzkin without
4280
 * introducing any spurious solutions.
4281
 *
4282
 * If at least one of the two constraints has a unit coefficient for the div,
4283
 * then the presence of such a value is guaranteed so there is no need to check.
4284
 * In particular, the value attained by the bound with unit coefficient
4285
 * can serve as this intermediate value.
4286
 */
4287
static struct isl_basic_map *drop_more_redundant_divs(
4288
  struct isl_basic_map *bmap, int *pairs, int n)
4289
708
{
4290
708
  isl_ctx *ctx;
4291
708
  struct test_ineq_data data = { NULL, NULL };
4292
708
  unsigned off, n_div;
4293
708
  int remove = -1;
4294
708
4295
708
  isl_int_init(data.g);
4296
708
  isl_int_init(data.fl);
4297
708
  isl_int_init(data.fu);
4298
708
4299
708
  if (!bmap)
4300
0
    goto error;
4301
708
4302
708
  ctx = isl_basic_map_get_ctx(bmap);
4303
708
  off = isl_basic_map_offset(bmap, isl_dim_div);
4304
708
  n_div = isl_basic_map_dim(bmap, isl_dim_div);
4305
708
  data.v = isl_vec_alloc(ctx, off + n_div);
4306
708
  if (!data.v)
4307
0
    goto error;
4308
708
4309
1.18k
  
while (708
n > 01.18k
)
{838
4310
838
    int i, l, u;
4311
838
    int best = -1;
4312
838
    isl_bool has_int;
4313
838
4314
2.47k
    for (i = 0; 
i < n_div2.47k
;
++i1.63k
)
{1.63k
4315
1.63k
      if (!pairs[i])
4316
512
        continue;
4317
1.12k
      
if (1.12k
best >= 0 && 1.12k
pairs[best] <= pairs[i]285
)
4318
215
        continue;
4319
908
      best = i;
4320
908
    }
4321
838
4322
838
    i = best;
4323
6.04k
    for (l = 0; 
l < bmap->n_ineq6.04k
;
++l5.20k
)
{5.68k
4324
5.68k
      if (
!5.68k
isl_int_is_pos5.68k
(bmap->ineq[l][off + i]))
4325
4.29k
        continue;
4326
1.39k
      
if (1.39k
isl_int_is_one1.39k
(bmap->ineq[l][off + i]))
4327
402
        continue;
4328
9.40k
      
for (u = 0; 990
u < bmap->n_ineq9.40k
;
++u8.41k
)
{8.88k
4329
8.88k
        if (
!8.88k
isl_int_is_neg8.88k
(bmap->ineq[u][off + i]))
4330
7.06k
          continue;
4331
1.82k
        
if (1.82k
isl_int_is_negone1.82k
(bmap->ineq[u][off + i]))
4332
380
          continue;
4333
1.44k
        has_int = int_between_bounds(bmap, i, l, u,
4334
1.44k
                &data);
4335
1.44k
        if (has_int < 0)
4336
0
          goto error;
4337
1.44k
        
if (1.44k
data.tab && 1.44k
data.tab->empty1.08k
)
4338
0
          break;
4339
1.44k
        
if (1.44k
!has_int1.44k
)
4340
474
          break;
4341
1.44k
      }
4342
990
      
if (990
u < bmap->n_ineq990
)
4343
474
        break;
4344
990
    }
4345
838
    
if (838
data.tab && 838
data.tab->empty530
)
{0
4346
0
      bmap = isl_basic_map_set_to_empty(bmap);
4347
0
      break;
4348
0
    }
4349
838
    
if (838
l == bmap->n_ineq838
)
{364
4350
364
      remove = i;
4351
364
      break;
4352
364
    }
4353
474
    pairs[i] = 0;
4354
474
    --n;
4355
474
  }
4356
708
4357
708
  test_ineq_data_clear(&data);
4358
708
4359
708
  free(pairs);
4360
708
4361
708
  if (remove < 0)
4362
344
    return bmap;
4363
708
4364
364
  bmap = isl_basic_map_remove_dims(bmap, isl_dim_div, remove, 1);
4365
364
  return isl_basic_map_drop_redundant_divs(bmap);
4366
0
error:
4367
0
  free(pairs);
4368
0
  isl_basic_map_free(bmap);
4369
0
  test_ineq_data_clear(&data);
4370
0
  return NULL;
4371
708
}
4372
4373
/* Given a pair of divs div1 and div2 such that, except for the lower bound l
4374
 * and the upper bound u, div1 always occurs together with div2 in the form
4375
 * (div1 + m div2), where m is the constant range on the variable div1
4376
 * allowed by l and u, replace the pair div1 and div2 by a single
4377
 * div that is equal to div1 + m div2.
4378
 *
4379
 * The new div will appear in the location that contains div2.
4380
 * We need to modify all constraints that contain
4381
 * div2 = (div - div1) / m
4382
 * The coefficient of div2 is known to be equal to 1 or -1.
4383
 * (If a constraint does not contain div2, it will also not contain div1.)
4384
 * If the constraint also contains div1, then we know they appear
4385
 * as f (div1 + m div2) and we can simply replace (div1 + m div2) by div,
4386
 * i.e., the coefficient of div is f.
4387
 *
4388
 * Otherwise, we first need to introduce div1 into the constraint.
4389
 * Let the l be
4390
 *
4391
 *  div1 + f >=0
4392
 *
4393
 * and u
4394
 *
4395
 *  -div1 + f' >= 0
4396
 *
4397
 * A lower bound on div2
4398
 *
4399
 *  div2 + t >= 0
4400
 *
4401
 * can be replaced by
4402
 *
4403
 *  m div2 + div1 + m t + f >= 0
4404
 *
4405
 * An upper bound
4406
 *
4407
 *  -div2 + t >= 0
4408
 *
4409
 * can be replaced by
4410
 *
4411
 *  -(m div2 + div1) + m t + f' >= 0
4412
 *
4413
 * These constraint are those that we would obtain from eliminating
4414
 * div1 using Fourier-Motzkin.
4415
 *
4416
 * After all constraints have been modified, we drop the lower and upper
4417
 * bound and then drop div1.
4418
 */
4419
static struct isl_basic_map *coalesce_divs(struct isl_basic_map *bmap,
4420
  unsigned div1, unsigned div2, unsigned l, unsigned u)
4421
1
{
4422
1
  isl_ctx *ctx;
4423
1
  isl_int m;
4424
1
  unsigned dim, total;
4425
1
  int i;
4426
1
4427
1
  ctx = isl_basic_map_get_ctx(bmap);
4428
1
4429
1
  dim = isl_space_dim(bmap->dim, isl_dim_all);
4430
1
  total = 1 + dim + bmap->n_div;
4431
1
4432
1
  isl_int_init(m);
4433
1
  isl_int_add(m, bmap->ineq[l][0], bmap->ineq[u][0]);
4434
1
  isl_int_add_ui(m, m, 1);
4435
1
4436
7
  for (i = 0; 
i < bmap->n_ineq7
;
++i6
)
{6
4437
6
    if (
i == l || 6
i == u5
)
4438
2
      continue;
4439
4
    
if (4
isl_int_is_zero4
(bmap->ineq[i][1 + dim + div2]))
4440
0
      continue;
4441
4
    
if (4
isl_int_is_zero4
(bmap->ineq[i][1 + dim + div1]))
{2
4442
2
      if (isl_int_is_pos(bmap->ineq[i][1 + dim + div2]))
4443
1
        isl_seq_combine(bmap->ineq[i], m, bmap->ineq[i],
4444
1
            ctx->one, bmap->ineq[l], total);
4445
2
      else
4446
1
        isl_seq_combine(bmap->ineq[i], m, bmap->ineq[i],
4447
1
            ctx->one, bmap->ineq[u], total);
4448
2
    }
4449
4
    isl_int_set(bmap->ineq[i][1 + dim + div2],
4450
4
          bmap->ineq[i][1 + dim + div1]);
4451
4
    isl_int_set_si(bmap->ineq[i][1 + dim + div1], 0);
4452
4
  }
4453
1
4454
1
  isl_int_clear(m);
4455
1
  if (
l > u1
)
{0
4456
0
    isl_basic_map_drop_inequality(bmap, l);
4457
0
    isl_basic_map_drop_inequality(bmap, u);
4458
1
  } else {
4459
1
    isl_basic_map_drop_inequality(bmap, u);
4460
1
    isl_basic_map_drop_inequality(bmap, l);
4461
1
  }
4462
1
  bmap = isl_basic_map_drop_div(bmap, div1);
4463
1
  return bmap;
4464
1
}
4465
4466
/* First check if we can coalesce any pair of divs and
4467
 * then continue with dropping more redundant divs.
4468
 *
4469
 * We loop over all pairs of lower and upper bounds on a div
4470
 * with coefficient 1 and -1, respectively, check if there
4471
 * is any other div "c" with which we can coalesce the div
4472
 * and if so, perform the coalescing.
4473
 */
4474
static struct isl_basic_map *coalesce_or_drop_more_redundant_divs(
4475
  struct isl_basic_map *bmap, int *pairs, int n)
4476
709
{
4477
709
  int i, l, u;
4478
709
  unsigned dim;
4479
709
4480
709
  dim = isl_space_dim(bmap->dim, isl_dim_all);
4481
709
4482
1.92k
  for (i = 0; 
i < bmap->n_div1.92k
;
++i1.21k
)
{1.21k
4483
1.21k
    if (!pairs[i])
4484
301
      continue;
4485
10.6k
    
for (l = 0; 918
l < bmap->n_ineq10.6k
;
++l9.70k
)
{9.70k
4486
9.70k
      if (
!9.70k
isl_int_is_one9.70k
(bmap->ineq[l][1 + dim + i]))
4487
9.23k
        continue;
4488
4.22k
      
for (u = 0; 471
u < bmap->n_ineq4.22k
;
++u3.75k
)
{3.75k
4489
3.75k
        int c;
4490
3.75k
4491
3.75k
        if (
!3.75k
isl_int_is_negone3.75k
(bmap->ineq[u][1+dim+i]))
4492
3.28k
          continue;
4493
474
        c = div_find_coalesce(bmap, pairs, i, l, u);
4494
474
        if (c < 0)
4495
473
          continue;
4496
1
        free(pairs);
4497
1
        bmap = coalesce_divs(bmap, i, c, l, u);
4498
1
        return isl_basic_map_drop_redundant_divs(bmap);
4499
474
      }
4500
471
    }
4501
918
  }
4502
709
4503
708
  
if (708
ISL_F_ISSET708
(bmap, ISL_BASIC_MAP_EMPTY))
{0
4504
0
    free(pairs);
4505
0
    return bmap;
4506
0
  }
4507
708
4508
708
  return drop_more_redundant_divs(bmap, pairs, n);
4509
708
}
4510
4511
/* Are the "n" coefficients starting at "first" of inequality constraints
4512
 * "i" and "j" of "bmap" equal to each other?
4513
 */
4514
static int is_parallel_part(__isl_keep isl_basic_map *bmap, int i, int j,
4515
  int first, int n)
4516
286
{
4517
286
  return isl_seq_eq(bmap->ineq[i] + first, bmap->ineq[j] + first, n);
4518
286
}
4519
4520
/* Are the "n" coefficients starting at "first" of inequality constraints
4521
 * "i" and "j" of "bmap" opposite to each other?
4522
 */
4523
static int is_opposite_part(__isl_keep isl_basic_map *bmap, int i, int j,
4524
  int first, int n)
4525
2.37k
{
4526
2.37k
  return isl_seq_is_neg(bmap->ineq[i] + first, bmap->ineq[j] + first, n);
4527
2.37k
}
4528
4529
/* Are inequality constraints "i" and "j" of "bmap" opposite to each other,
4530
 * apart from the constant term?
4531
 */
4532
static isl_bool is_opposite(__isl_keep isl_basic_map *bmap, int i, int j)
4533
2.04k
{
4534
2.04k
  unsigned total;
4535
2.04k
4536
2.04k
  total = isl_basic_map_dim(bmap, isl_dim_all);
4537
2.04k
  return is_opposite_part(bmap, i, j, 1, total);
4538
2.04k
}
4539
4540
/* Are inequality constraints "i" and "j" of "bmap" equal to each other,
4541
 * apart from the constant term and the coefficient at position "pos"?
4542
 */
4543
static int is_parallel_except(__isl_keep isl_basic_map *bmap, int i, int j,
4544
  int pos)
4545
269
{
4546
269
  unsigned total;
4547
269
4548
269
  total = isl_basic_map_dim(bmap, isl_dim_all);
4549
269
  return is_parallel_part(bmap, i, j, 1, pos - 1) &&
4550
17
    is_parallel_part(bmap, i, j, pos + 1, total - pos);
4551
269
}
4552
4553
/* Are inequality constraints "i" and "j" of "bmap" opposite to each other,
4554
 * apart from the constant term and the coefficient at position "pos"?
4555
 */
4556
static int is_opposite_except(__isl_keep isl_basic_map *bmap, int i, int j,
4557
  int pos)
4558
271
{
4559
271
  unsigned total;
4560
271
4561
271
  total = isl_basic_map_dim(bmap, isl_dim_all);
4562
271
  return is_opposite_part(bmap, i, j, 1, pos - 1) &&
4563
61
    is_opposite_part(bmap, i, j, pos + 1, total - pos);
4564
271
}
4565
4566
/* Restart isl_basic_map_drop_redundant_divs after "bmap" has
4567
 * been modified, simplying it if "simplify" is set.
4568
 * Free the temporary data structure "pairs" that was associated
4569
 * to the old version of "bmap".
4570
 */
4571
static __isl_give isl_basic_map *drop_redundant_divs_again(
4572
  __isl_take isl_basic_map *bmap, __isl_take int *pairs, int simplify)
4573
2.51k
{
4574
2.51k
  if (simplify)
4575
902
    bmap = isl_basic_map_simplify(bmap);
4576
2.51k
  free(pairs);
4577
2.51k
  return isl_basic_map_drop_redundant_divs(bmap);
4578
2.51k
}
4579
4580
/* Is "div" the single unknown existentially quantified variable
4581
 * in inequality constraint "ineq" of "bmap"?
4582
 * "div" is known to have a non-zero coefficient in "ineq".
4583
 */
4584
static isl_bool single_unknown(__isl_keep isl_basic_map *bmap, int ineq,
4585
  int div)
4586
1.09k
{
4587
1.09k
  int i;
4588
1.09k
  unsigned n_div, o_div;
4589
1.09k
  isl_bool known;
4590
1.09k
4591
1.09k
  known = isl_basic_map_div_is_known(bmap, div);
4592
1.09k
  if (
known < 0 || 1.09k
known1.09k
)
4593
61
    return isl_bool_not(known);
4594
1.03k
  n_div = isl_basic_map_dim(bmap, isl_dim_div);
4595
1.03k
  if (n_div == 1)
4596
729
    return isl_bool_true;
4597
301
  o_div = isl_basic_map_offset(bmap, isl_dim_div);
4598
981
  for (i = 0; 
i < n_div981
;
++i680
)
{714
4599
714
    isl_bool known;
4600
714
4601
714
    if (i == div)
4602
284
      continue;
4603
430
    
if (430
isl_int_is_zero430
(bmap->ineq[ineq][o_div + i]))
4604
384
      continue;
4605
46
    known = isl_basic_map_div_is_known(bmap, i);
4606
46
    if (
known < 0 || 46
!known46
)
4607
34
      return known;
4608
46
  }
4609
301
4610
267
  return isl_bool_true;
4611
301
}
4612
4613
/* Does integer division "div" have coefficient 1 in inequality constraint
4614
 * "ineq" of "map"?
4615
 */
4616
static isl_bool has_coef_one(__isl_keep isl_basic_map *bmap, int div, int ineq)
4617
996
{
4618
996
  unsigned o_div;
4619
996
4620
996
  o_div = isl_basic_map_offset(bmap, isl_dim_div);
4621
996
  if (isl_int_is_one(bmap->ineq[ineq][o_div + div]))
4622
899
    return isl_bool_true;
4623
996
4624
97
  return isl_bool_false;
4625
996
}
4626
4627
/* Turn inequality constraint "ineq" of "bmap" into an equality and
4628
 * then try and drop redundant divs again,
4629
 * freeing the temporary data structure "pairs" that was associated
4630
 * to the old version of "bmap".
4631
 */
4632
static __isl_give isl_basic_map *set_eq_and_try_again(
4633
  __isl_take isl_basic_map *bmap, int ineq, __isl_take int *pairs)
4634
899
{
4635
899
  bmap = isl_basic_map_cow(bmap);
4636
899
  isl_basic_map_inequality_to_equality(bmap, ineq);
4637
899
  return drop_redundant_divs_again(bmap, pairs, 1);
4638
899
}
4639
4640
/* Drop the integer division at position "div", along with the two
4641
 * inequality constraints "ineq1" and "ineq2" in which it appears
4642
 * from "bmap" and then try and drop redundant divs again,
4643
 * freeing the temporary data structure "pairs" that was associated
4644
 * to the old version of "bmap".
4645
 */
4646
static __isl_give isl_basic_map *drop_div_and_try_again(
4647
  __isl_take isl_basic_map *bmap, int div, int ineq1, int ineq2,
4648
  __isl_take int *pairs)
4649
518
{
4650
518
  if (
ineq1 > ineq2518
)
{371
4651
371
    isl_basic_map_drop_inequality(bmap, ineq1);
4652
371
    isl_basic_map_drop_inequality(bmap, ineq2);
4653
147
  } else {
4654
147
    isl_basic_map_drop_inequality(bmap, ineq2);
4655
147
    isl_basic_map_drop_inequality(bmap, ineq1);
4656
147
  }
4657
518
  bmap = isl_basic_map_drop_div(bmap, div);
4658
518
  return drop_redundant_divs_again(bmap, pairs, 0);
4659
518
}
4660
4661
/* Given two inequality constraints
4662
 *
4663
 *  f(x) + n d + c >= 0,    (ineq)
4664
 *
4665
 * with d the variable at position "pos", and
4666
 *
4667
 *  f(x) + c0 >= 0,     (lower)
4668
 *
4669
 * compute the maximal value of the lower bound ceil((-f(x) - c)/n)
4670
 * determined by the first constraint.
4671
 * That is, store
4672
 *
4673
 *  ceil((c0 - c)/n)
4674
 *
4675
 * in *l.
4676
 */
4677
static void lower_bound_from_parallel(__isl_keep isl_basic_map *bmap,
4678
  int ineq, int lower, int pos, isl_int *l)
4679
9
{
4680
9
  isl_int_neg(*l, bmap->ineq[ineq][0]);
4681
9
  isl_int_add(*l, *l, bmap->ineq[lower][0]);
4682
9
  isl_int_cdiv_q(*l, *l, bmap->ineq[ineq][pos]);
4683
9
}
4684
4685
/* Given two inequality constraints
4686
 *
4687
 *  f(x) + n d + c >= 0,    (ineq)
4688
 *
4689
 * with d the variable at position "pos", and
4690
 *
4691
 *  -f(x) - c0 >= 0,    (upper)
4692
 *
4693
 * compute the minimal value of the lower bound ceil((-f(x) - c)/n)
4694
 * determined by the first constraint.
4695
 * That is, store
4696
 *
4697
 *  ceil((-c1 - c)/n)
4698
 *
4699
 * in *u.
4700
 */
4701
static void lower_bound_from_opposite(__isl_keep isl_basic_map *bmap,
4702
  int ineq, int upper, int pos, isl_int *u)
4703
7
{
4704
7
  isl_int_neg(*u, bmap->ineq[ineq][0]);
4705
7
  isl_int_sub(*u, *u, bmap->ineq[upper][0]);
4706
7
  isl_int_cdiv_q(*u, *u, bmap->ineq[ineq][pos]);
4707
7
}
4708
4709
/* Given a lower bound constraint "ineq" on "div" in "bmap",
4710
 * does the corresponding lower bound have a fixed value in "bmap"?
4711
 *
4712
 * In particular, "ineq" is of the form
4713
 *
4714
 *  f(x) + n d + c >= 0
4715
 *
4716
 * with n > 0, c the constant term and
4717
 * d the existentially quantified variable "div".
4718
 * That is, the lower bound is
4719
 *
4720
 *  ceil((-f(x) - c)/n)
4721
 *
4722
 * Look for a pair of constraints
4723
 *
4724
 *  f(x) + c0 >= 0
4725
 *  -f(x) + c1 >= 0
4726
 *
4727
 * i.e., -c1 <= -f(x) <= c0, that fix ceil((-f(x) - c)/n) to a constant value.
4728
 * That is, check that
4729
 *
4730
 *  ceil((-c1 - c)/n) = ceil((c0 - c)/n)
4731
 *
4732
 * If so, return the index of inequality f(x) + c0 >= 0.
4733
 * Otherwise, return -1.
4734
 */
4735
static int lower_bound_is_cst(__isl_keep isl_basic_map *bmap, int div, int ineq)
4736
97
{
4737
97
  int i;
4738
97
  int lower = -1, upper = -1;
4739
97
  unsigned o_div;
4740
97
  isl_int l, u;
4741
97
  int equal;
4742
97
4743
97
  o_div = isl_basic_map_offset(bmap, isl_dim_div);
4744
604
  for (i = 0; 
i < bmap->n_ineq && 604
(lower < 0 || 512
upper < 049
);
++i507
)
{507
4745
507
    if (i == ineq)
4746
92
      continue;
4747
415
    
if (415
!415
isl_int_is_zero415
(bmap->ineq[i][o_div + div]))
4748
118
      continue;
4749
297
    
if (297
lower < 0 &&297
4750
269
        
is_parallel_except(bmap, ineq, i, o_div + div)269
)
{13
4751
13
      lower = i;
4752
13
      continue;
4753
13
    }
4754
284
    
if (284
upper < 0 &&284
4755
271
        
is_opposite_except(bmap, ineq, i, o_div + div)271
)
{58
4756
58
      upper = i;
4757
58
    }
4758
284
  }
4759
97
4760
97
  if (
lower < 0 || 97
upper < 013
)
4761
90
    return -1;
4762
97
4763
7
  
isl_int_init7
(l);7
4764
7
  isl_int_init(u);
4765
7
4766
7
  lower_bound_from_parallel(bmap, ineq, lower, o_div + div, &l);
4767
7
  lower_bound_from_opposite(bmap, ineq, upper, o_div + div, &u);
4768
7
4769
7
  equal = isl_int_eq(l, u);
4770
7
4771
7
  isl_int_clear(l);
4772
7
  isl_int_clear(u);
4773
7
4774
5
  return equal ? 
lower2
:
-15
;
4775
97
}
4776
4777
/* Given a lower bound constraint "ineq" on the existentially quantified
4778
 * variable "div", such that the corresponding lower bound has
4779
 * a fixed value in "bmap", assign this fixed value to the variable and
4780
 * then try and drop redundant divs again,
4781
 * freeing the temporary data structure "pairs" that was associated
4782
 * to the old version of "bmap".
4783
 * "lower" determines the constant value for the lower bound.
4784
 *
4785
 * In particular, "ineq" is of the form
4786
 *
4787
 *  f(x) + n d + c >= 0,
4788
 *
4789
 * while "lower" is of the form
4790
 *
4791
 *  f(x) + c0 >= 0
4792
 *
4793
 * The lower bound is ceil((-f(x) - c)/n) and its constant value
4794
 * is ceil((c0 - c)/n).
4795
 */
4796
static __isl_give isl_basic_map *fix_cst_lower(__isl_take isl_basic_map *bmap,
4797
  int div, int ineq, int lower, int *pairs)
4798
2
{
4799
2
  isl_int c;
4800
2
  unsigned o_div;
4801
2
4802
2
  isl_int_init(c);
4803
2
4804
2
  o_div = isl_basic_map_offset(bmap, isl_dim_div);
4805
2
  lower_bound_from_parallel(bmap, ineq, lower, o_div + div, &c);
4806
2
  bmap = isl_basic_map_fix(bmap, isl_dim_div, div, c);
4807
2
  free(pairs);
4808
2
4809
2
  isl_int_clear(c);
4810
2
4811
2
  return isl_basic_map_drop_redundant_divs(bmap);
4812
2
}
4813
4814
/* Remove divs that are not strictly needed based on the inequality
4815
 * constraints.
4816
 * In particular, if a div only occurs positively (or negatively)
4817
 * in constraints, then it can simply be dropped.