add isl_val_abs_eq
[isl.git] / isl_map_simplify.c
blob30e2172dd7cd5fef4c55de99d096f415c6c82114
1 /*
2 * Copyright 2008-2009 Katholieke Universiteit Leuven
3 * Copyright 2012 Ecole Normale Superieure
4 * Copyright 2014 INRIA Rocquencourt
6 * Use of this software is governed by the MIT license
8 * Written by Sven Verdoolaege, K.U.Leuven, Departement
9 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
10 * and Ecole Normale Superieure, 45 rue d’Ulm, 75230 Paris, France
11 * and Inria Paris - Rocquencourt, Domaine de Voluceau - Rocquencourt,
12 * B.P. 105 - 78153 Le Chesnay, France
15 #include <strings.h>
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>
26 static void swap_equality(struct isl_basic_map *bmap, int a, int b)
28 isl_int *t = bmap->eq[a];
29 bmap->eq[a] = bmap->eq[b];
30 bmap->eq[b] = t;
33 static void swap_inequality(struct isl_basic_map *bmap, int a, int b)
35 if (a != b) {
36 isl_int *t = bmap->ineq[a];
37 bmap->ineq[a] = bmap->ineq[b];
38 bmap->ineq[b] = t;
42 static void constraint_drop_vars(isl_int *c, unsigned n, unsigned rem)
44 isl_seq_cpy(c, c + n, rem);
45 isl_seq_clr(c + rem, n);
48 /* Drop n dimensions starting at first.
50 * In principle, this frees up some extra variables as the number
51 * of columns remains constant, but we would have to extend
52 * the div array too as the number of rows in this array is assumed
53 * to be equal to extra.
55 struct isl_basic_set *isl_basic_set_drop_dims(
56 struct isl_basic_set *bset, unsigned first, unsigned n)
58 int i;
60 if (!bset)
61 goto error;
63 isl_assert(bset->ctx, first + n <= bset->dim->n_out, goto error);
65 if (n == 0 && !isl_space_get_tuple_name(bset->dim, isl_dim_set))
66 return bset;
68 bset = isl_basic_set_cow(bset);
69 if (!bset)
70 return NULL;
72 for (i = 0; i < bset->n_eq; ++i)
73 constraint_drop_vars(bset->eq[i]+1+bset->dim->nparam+first, n,
74 (bset->dim->n_out-first-n)+bset->extra);
76 for (i = 0; i < bset->n_ineq; ++i)
77 constraint_drop_vars(bset->ineq[i]+1+bset->dim->nparam+first, n,
78 (bset->dim->n_out-first-n)+bset->extra);
80 for (i = 0; i < bset->n_div; ++i)
81 constraint_drop_vars(bset->div[i]+1+1+bset->dim->nparam+first, n,
82 (bset->dim->n_out-first-n)+bset->extra);
84 bset->dim = isl_space_drop_outputs(bset->dim, first, n);
85 if (!bset->dim)
86 goto error;
88 ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED);
89 bset = isl_basic_set_simplify(bset);
90 return isl_basic_set_finalize(bset);
91 error:
92 isl_basic_set_free(bset);
93 return NULL;
96 struct isl_set *isl_set_drop_dims(
97 struct isl_set *set, unsigned first, unsigned n)
99 int i;
101 if (!set)
102 goto error;
104 isl_assert(set->ctx, first + n <= set->dim->n_out, goto error);
106 if (n == 0 && !isl_space_get_tuple_name(set->dim, isl_dim_set))
107 return set;
108 set = isl_set_cow(set);
109 if (!set)
110 goto error;
111 set->dim = isl_space_drop_outputs(set->dim, first, n);
112 if (!set->dim)
113 goto error;
115 for (i = 0; i < set->n; ++i) {
116 set->p[i] = isl_basic_set_drop_dims(set->p[i], first, n);
117 if (!set->p[i])
118 goto error;
121 ISL_F_CLR(set, ISL_SET_NORMALIZED);
122 return set;
123 error:
124 isl_set_free(set);
125 return NULL;
128 /* Move "n" divs starting at "first" to the end of the list of divs.
130 static struct isl_basic_map *move_divs_last(struct isl_basic_map *bmap,
131 unsigned first, unsigned n)
133 isl_int **div;
134 int i;
136 if (first + n == bmap->n_div)
137 return bmap;
139 div = isl_alloc_array(bmap->ctx, isl_int *, n);
140 if (!div)
141 goto error;
142 for (i = 0; i < n; ++i)
143 div[i] = bmap->div[first + i];
144 for (i = 0; i < bmap->n_div - first - n; ++i)
145 bmap->div[first + i] = bmap->div[first + n + i];
146 for (i = 0; i < n; ++i)
147 bmap->div[bmap->n_div - n + i] = div[i];
148 free(div);
149 return bmap;
150 error:
151 isl_basic_map_free(bmap);
152 return NULL;
155 /* Drop "n" dimensions of type "type" starting at "first".
157 * In principle, this frees up some extra variables as the number
158 * of columns remains constant, but we would have to extend
159 * the div array too as the number of rows in this array is assumed
160 * to be equal to extra.
162 struct isl_basic_map *isl_basic_map_drop(struct isl_basic_map *bmap,
163 enum isl_dim_type type, unsigned first, unsigned n)
165 int i;
166 unsigned dim;
167 unsigned offset;
168 unsigned left;
170 if (!bmap)
171 goto error;
173 dim = isl_basic_map_dim(bmap, type);
174 isl_assert(bmap->ctx, first + n <= dim, goto error);
176 if (n == 0 && !isl_space_is_named_or_nested(bmap->dim, type))
177 return bmap;
179 bmap = isl_basic_map_cow(bmap);
180 if (!bmap)
181 return NULL;
183 offset = isl_basic_map_offset(bmap, type) + first;
184 left = isl_basic_map_total_dim(bmap) - (offset - 1) - n;
185 for (i = 0; i < bmap->n_eq; ++i)
186 constraint_drop_vars(bmap->eq[i]+offset, n, left);
188 for (i = 0; i < bmap->n_ineq; ++i)
189 constraint_drop_vars(bmap->ineq[i]+offset, n, left);
191 for (i = 0; i < bmap->n_div; ++i)
192 constraint_drop_vars(bmap->div[i]+1+offset, n, left);
194 if (type == isl_dim_div) {
195 bmap = move_divs_last(bmap, first, n);
196 if (!bmap)
197 goto error;
198 isl_basic_map_free_div(bmap, n);
199 } else
200 bmap->dim = isl_space_drop_dims(bmap->dim, type, first, n);
201 if (!bmap->dim)
202 goto error;
204 ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
205 bmap = isl_basic_map_simplify(bmap);
206 return isl_basic_map_finalize(bmap);
207 error:
208 isl_basic_map_free(bmap);
209 return NULL;
212 __isl_give isl_basic_set *isl_basic_set_drop(__isl_take isl_basic_set *bset,
213 enum isl_dim_type type, unsigned first, unsigned n)
215 return (isl_basic_set *)isl_basic_map_drop((isl_basic_map *)bset,
216 type, first, n);
219 struct isl_basic_map *isl_basic_map_drop_inputs(
220 struct isl_basic_map *bmap, unsigned first, unsigned n)
222 return isl_basic_map_drop(bmap, isl_dim_in, first, n);
225 struct isl_map *isl_map_drop(struct isl_map *map,
226 enum isl_dim_type type, unsigned first, unsigned n)
228 int i;
230 if (!map)
231 goto error;
233 isl_assert(map->ctx, first + n <= isl_map_dim(map, type), goto error);
235 if (n == 0 && !isl_space_get_tuple_name(map->dim, type))
236 return map;
237 map = isl_map_cow(map);
238 if (!map)
239 goto error;
240 map->dim = isl_space_drop_dims(map->dim, type, first, n);
241 if (!map->dim)
242 goto error;
244 for (i = 0; i < map->n; ++i) {
245 map->p[i] = isl_basic_map_drop(map->p[i], type, first, n);
246 if (!map->p[i])
247 goto error;
249 ISL_F_CLR(map, ISL_MAP_NORMALIZED);
251 return map;
252 error:
253 isl_map_free(map);
254 return NULL;
257 struct isl_set *isl_set_drop(struct isl_set *set,
258 enum isl_dim_type type, unsigned first, unsigned n)
260 return (isl_set *)isl_map_drop((isl_map *)set, type, first, n);
263 struct isl_map *isl_map_drop_inputs(
264 struct isl_map *map, unsigned first, unsigned n)
266 return isl_map_drop(map, isl_dim_in, first, n);
270 * We don't cow, as the div is assumed to be redundant.
272 static struct isl_basic_map *isl_basic_map_drop_div(
273 struct isl_basic_map *bmap, unsigned div)
275 int i;
276 unsigned pos;
278 if (!bmap)
279 goto error;
281 pos = 1 + isl_space_dim(bmap->dim, isl_dim_all) + div;
283 isl_assert(bmap->ctx, div < bmap->n_div, goto error);
285 for (i = 0; i < bmap->n_eq; ++i)
286 constraint_drop_vars(bmap->eq[i]+pos, 1, bmap->extra-div-1);
288 for (i = 0; i < bmap->n_ineq; ++i) {
289 if (!isl_int_is_zero(bmap->ineq[i][pos])) {
290 isl_basic_map_drop_inequality(bmap, i);
291 --i;
292 continue;
294 constraint_drop_vars(bmap->ineq[i]+pos, 1, bmap->extra-div-1);
297 for (i = 0; i < bmap->n_div; ++i)
298 constraint_drop_vars(bmap->div[i]+1+pos, 1, bmap->extra-div-1);
300 if (div != bmap->n_div - 1) {
301 int j;
302 isl_int *t = bmap->div[div];
304 for (j = div; j < bmap->n_div - 1; ++j)
305 bmap->div[j] = bmap->div[j+1];
307 bmap->div[bmap->n_div - 1] = t;
309 ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
310 isl_basic_map_free_div(bmap, 1);
312 return bmap;
313 error:
314 isl_basic_map_free(bmap);
315 return NULL;
318 struct isl_basic_map *isl_basic_map_normalize_constraints(
319 struct isl_basic_map *bmap)
321 int i;
322 isl_int gcd;
323 unsigned total = isl_basic_map_total_dim(bmap);
325 if (!bmap)
326 return NULL;
328 isl_int_init(gcd);
329 for (i = bmap->n_eq - 1; i >= 0; --i) {
330 isl_seq_gcd(bmap->eq[i]+1, total, &gcd);
331 if (isl_int_is_zero(gcd)) {
332 if (!isl_int_is_zero(bmap->eq[i][0])) {
333 bmap = isl_basic_map_set_to_empty(bmap);
334 break;
336 isl_basic_map_drop_equality(bmap, i);
337 continue;
339 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
340 isl_int_gcd(gcd, gcd, bmap->eq[i][0]);
341 if (isl_int_is_one(gcd))
342 continue;
343 if (!isl_int_is_divisible_by(bmap->eq[i][0], gcd)) {
344 bmap = isl_basic_map_set_to_empty(bmap);
345 break;
347 isl_seq_scale_down(bmap->eq[i], bmap->eq[i], gcd, 1+total);
350 for (i = bmap->n_ineq - 1; i >= 0; --i) {
351 isl_seq_gcd(bmap->ineq[i]+1, total, &gcd);
352 if (isl_int_is_zero(gcd)) {
353 if (isl_int_is_neg(bmap->ineq[i][0])) {
354 bmap = isl_basic_map_set_to_empty(bmap);
355 break;
357 isl_basic_map_drop_inequality(bmap, i);
358 continue;
360 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
361 isl_int_gcd(gcd, gcd, bmap->ineq[i][0]);
362 if (isl_int_is_one(gcd))
363 continue;
364 isl_int_fdiv_q(bmap->ineq[i][0], bmap->ineq[i][0], gcd);
365 isl_seq_scale_down(bmap->ineq[i]+1, bmap->ineq[i]+1, gcd, total);
367 isl_int_clear(gcd);
369 return bmap;
372 struct isl_basic_set *isl_basic_set_normalize_constraints(
373 struct isl_basic_set *bset)
375 return (struct isl_basic_set *)isl_basic_map_normalize_constraints(
376 (struct isl_basic_map *)bset);
379 /* Remove any common factor in numerator and denominator of the div expression,
380 * not taking into account the constant term.
381 * That is, if the div is of the form
383 * floor((a + m f(x))/(m d))
385 * then replace it by
387 * floor((floor(a/m) + f(x))/d)
389 * The difference {a/m}/d in the argument satisfies 0 <= {a/m}/d < 1/d
390 * and can therefore not influence the result of the floor.
392 static void normalize_div_expression(__isl_keep isl_basic_map *bmap, int div)
394 unsigned total = isl_basic_map_total_dim(bmap);
395 isl_ctx *ctx = bmap->ctx;
397 if (isl_int_is_zero(bmap->div[div][0]))
398 return;
399 isl_seq_gcd(bmap->div[div] + 2, total, &ctx->normalize_gcd);
400 isl_int_gcd(ctx->normalize_gcd, ctx->normalize_gcd, bmap->div[div][0]);
401 if (isl_int_is_one(ctx->normalize_gcd))
402 return;
403 isl_int_fdiv_q(bmap->div[div][1], bmap->div[div][1],
404 ctx->normalize_gcd);
405 isl_int_divexact(bmap->div[div][0], bmap->div[div][0],
406 ctx->normalize_gcd);
407 isl_seq_scale_down(bmap->div[div] + 2, bmap->div[div] + 2,
408 ctx->normalize_gcd, total);
411 /* Remove any common factor in numerator and denominator of a div expression,
412 * not taking into account the constant term.
413 * That is, look for any div of the form
415 * floor((a + m f(x))/(m d))
417 * and replace it by
419 * floor((floor(a/m) + f(x))/d)
421 * The difference {a/m}/d in the argument satisfies 0 <= {a/m}/d < 1/d
422 * and can therefore not influence the result of the floor.
424 static __isl_give isl_basic_map *normalize_div_expressions(
425 __isl_take isl_basic_map *bmap)
427 int i;
429 if (!bmap)
430 return NULL;
431 if (bmap->n_div == 0)
432 return bmap;
434 for (i = 0; i < bmap->n_div; ++i)
435 normalize_div_expression(bmap, i);
437 return bmap;
440 /* Assumes divs have been ordered if keep_divs is set.
442 static void eliminate_var_using_equality(struct isl_basic_map *bmap,
443 unsigned pos, isl_int *eq, int keep_divs, int *progress)
445 unsigned total;
446 unsigned space_total;
447 int k;
448 int last_div;
450 total = isl_basic_map_total_dim(bmap);
451 space_total = isl_space_dim(bmap->dim, isl_dim_all);
452 last_div = isl_seq_last_non_zero(eq + 1 + space_total, bmap->n_div);
453 for (k = 0; k < bmap->n_eq; ++k) {
454 if (bmap->eq[k] == eq)
455 continue;
456 if (isl_int_is_zero(bmap->eq[k][1+pos]))
457 continue;
458 if (progress)
459 *progress = 1;
460 isl_seq_elim(bmap->eq[k], eq, 1+pos, 1+total, NULL);
461 isl_seq_normalize(bmap->ctx, bmap->eq[k], 1 + total);
464 for (k = 0; k < bmap->n_ineq; ++k) {
465 if (isl_int_is_zero(bmap->ineq[k][1+pos]))
466 continue;
467 if (progress)
468 *progress = 1;
469 isl_seq_elim(bmap->ineq[k], eq, 1+pos, 1+total, NULL);
470 isl_seq_normalize(bmap->ctx, bmap->ineq[k], 1 + total);
471 ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
474 for (k = 0; k < bmap->n_div; ++k) {
475 if (isl_int_is_zero(bmap->div[k][0]))
476 continue;
477 if (isl_int_is_zero(bmap->div[k][1+1+pos]))
478 continue;
479 if (progress)
480 *progress = 1;
481 /* We need to be careful about circular definitions,
482 * so for now we just remove the definition of div k
483 * if the equality contains any divs.
484 * If keep_divs is set, then the divs have been ordered
485 * and we can keep the definition as long as the result
486 * is still ordered.
488 if (last_div == -1 || (keep_divs && last_div < k)) {
489 isl_seq_elim(bmap->div[k]+1, eq,
490 1+pos, 1+total, &bmap->div[k][0]);
491 normalize_div_expression(bmap, k);
492 } else
493 isl_seq_clr(bmap->div[k], 1 + total);
494 ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
498 /* Assumes divs have been ordered if keep_divs is set.
500 static void eliminate_div(struct isl_basic_map *bmap, isl_int *eq,
501 unsigned div, int keep_divs)
503 unsigned pos = isl_space_dim(bmap->dim, isl_dim_all) + div;
505 eliminate_var_using_equality(bmap, pos, eq, keep_divs, NULL);
507 isl_basic_map_drop_div(bmap, div);
510 /* Check if elimination of div "div" using equality "eq" would not
511 * result in a div depending on a later div.
513 static int ok_to_eliminate_div(struct isl_basic_map *bmap, isl_int *eq,
514 unsigned div)
516 int k;
517 int last_div;
518 unsigned space_total = isl_space_dim(bmap->dim, isl_dim_all);
519 unsigned pos = space_total + div;
521 last_div = isl_seq_last_non_zero(eq + 1 + space_total, bmap->n_div);
522 if (last_div < 0 || last_div <= div)
523 return 1;
525 for (k = 0; k <= last_div; ++k) {
526 if (isl_int_is_zero(bmap->div[k][0]))
527 return 1;
528 if (!isl_int_is_zero(bmap->div[k][1 + 1 + pos]))
529 return 0;
532 return 1;
535 /* Elimininate divs based on equalities
537 static struct isl_basic_map *eliminate_divs_eq(
538 struct isl_basic_map *bmap, int *progress)
540 int d;
541 int i;
542 int modified = 0;
543 unsigned off;
545 bmap = isl_basic_map_order_divs(bmap);
547 if (!bmap)
548 return NULL;
550 off = 1 + isl_space_dim(bmap->dim, isl_dim_all);
552 for (d = bmap->n_div - 1; d >= 0 ; --d) {
553 for (i = 0; i < bmap->n_eq; ++i) {
554 if (!isl_int_is_one(bmap->eq[i][off + d]) &&
555 !isl_int_is_negone(bmap->eq[i][off + d]))
556 continue;
557 if (!ok_to_eliminate_div(bmap, bmap->eq[i], d))
558 continue;
559 modified = 1;
560 *progress = 1;
561 eliminate_div(bmap, bmap->eq[i], d, 1);
562 isl_basic_map_drop_equality(bmap, i);
563 break;
566 if (modified)
567 return eliminate_divs_eq(bmap, progress);
568 return bmap;
571 /* Elimininate divs based on inequalities
573 static struct isl_basic_map *eliminate_divs_ineq(
574 struct isl_basic_map *bmap, int *progress)
576 int d;
577 int i;
578 unsigned off;
579 struct isl_ctx *ctx;
581 if (!bmap)
582 return NULL;
584 ctx = bmap->ctx;
585 off = 1 + isl_space_dim(bmap->dim, isl_dim_all);
587 for (d = bmap->n_div - 1; d >= 0 ; --d) {
588 for (i = 0; i < bmap->n_eq; ++i)
589 if (!isl_int_is_zero(bmap->eq[i][off + d]))
590 break;
591 if (i < bmap->n_eq)
592 continue;
593 for (i = 0; i < bmap->n_ineq; ++i)
594 if (isl_int_abs_gt(bmap->ineq[i][off + d], ctx->one))
595 break;
596 if (i < bmap->n_ineq)
597 continue;
598 *progress = 1;
599 bmap = isl_basic_map_eliminate_vars(bmap, (off-1)+d, 1);
600 if (!bmap || ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
601 break;
602 bmap = isl_basic_map_drop_div(bmap, d);
603 if (!bmap)
604 break;
606 return bmap;
609 struct isl_basic_map *isl_basic_map_gauss(
610 struct isl_basic_map *bmap, int *progress)
612 int k;
613 int done;
614 int last_var;
615 unsigned total_var;
616 unsigned total;
618 bmap = isl_basic_map_order_divs(bmap);
620 if (!bmap)
621 return NULL;
623 total = isl_basic_map_total_dim(bmap);
624 total_var = total - bmap->n_div;
626 last_var = total - 1;
627 for (done = 0; done < bmap->n_eq; ++done) {
628 for (; last_var >= 0; --last_var) {
629 for (k = done; k < bmap->n_eq; ++k)
630 if (!isl_int_is_zero(bmap->eq[k][1+last_var]))
631 break;
632 if (k < bmap->n_eq)
633 break;
635 if (last_var < 0)
636 break;
637 if (k != done)
638 swap_equality(bmap, k, done);
639 if (isl_int_is_neg(bmap->eq[done][1+last_var]))
640 isl_seq_neg(bmap->eq[done], bmap->eq[done], 1+total);
642 eliminate_var_using_equality(bmap, last_var, bmap->eq[done], 1,
643 progress);
645 if (last_var >= total_var &&
646 isl_int_is_zero(bmap->div[last_var - total_var][0])) {
647 unsigned div = last_var - total_var;
648 isl_seq_neg(bmap->div[div]+1, bmap->eq[done], 1+total);
649 isl_int_set_si(bmap->div[div][1+1+last_var], 0);
650 isl_int_set(bmap->div[div][0],
651 bmap->eq[done][1+last_var]);
652 if (progress)
653 *progress = 1;
654 ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
657 if (done == bmap->n_eq)
658 return bmap;
659 for (k = done; k < bmap->n_eq; ++k) {
660 if (isl_int_is_zero(bmap->eq[k][0]))
661 continue;
662 return isl_basic_map_set_to_empty(bmap);
664 isl_basic_map_free_equality(bmap, bmap->n_eq-done);
665 return bmap;
668 struct isl_basic_set *isl_basic_set_gauss(
669 struct isl_basic_set *bset, int *progress)
671 return (struct isl_basic_set*)isl_basic_map_gauss(
672 (struct isl_basic_map *)bset, progress);
676 static unsigned int round_up(unsigned int v)
678 int old_v = v;
680 while (v) {
681 old_v = v;
682 v ^= v & -v;
684 return old_v << 1;
687 static int hash_index(isl_int ***index, unsigned int size, int bits,
688 struct isl_basic_map *bmap, int k)
690 int h;
691 unsigned total = isl_basic_map_total_dim(bmap);
692 uint32_t hash = isl_seq_get_hash_bits(bmap->ineq[k]+1, total, bits);
693 for (h = hash; index[h]; h = (h+1) % size)
694 if (&bmap->ineq[k] != index[h] &&
695 isl_seq_eq(bmap->ineq[k]+1, index[h][0]+1, total))
696 break;
697 return h;
700 static int set_hash_index(isl_int ***index, unsigned int size, int bits,
701 struct isl_basic_set *bset, int k)
703 return hash_index(index, size, bits, (struct isl_basic_map *)bset, k);
706 /* If we can eliminate more than one div, then we need to make
707 * sure we do it from last div to first div, in order not to
708 * change the position of the other divs that still need to
709 * be removed.
711 static struct isl_basic_map *remove_duplicate_divs(
712 struct isl_basic_map *bmap, int *progress)
714 unsigned int size;
715 int *index;
716 int *elim_for;
717 int k, l, h;
718 int bits;
719 struct isl_blk eq;
720 unsigned total_var;
721 unsigned total;
722 struct isl_ctx *ctx;
724 bmap = isl_basic_map_order_divs(bmap);
725 if (!bmap || bmap->n_div <= 1)
726 return bmap;
728 total_var = isl_space_dim(bmap->dim, isl_dim_all);
729 total = total_var + bmap->n_div;
731 ctx = bmap->ctx;
732 for (k = bmap->n_div - 1; k >= 0; --k)
733 if (!isl_int_is_zero(bmap->div[k][0]))
734 break;
735 if (k <= 0)
736 return bmap;
738 elim_for = isl_calloc_array(ctx, int, bmap->n_div);
739 size = round_up(4 * bmap->n_div / 3 - 1);
740 bits = ffs(size) - 1;
741 index = isl_calloc_array(ctx, int, size);
742 if (!index)
743 return bmap;
744 eq = isl_blk_alloc(ctx, 1+total);
745 if (isl_blk_is_error(eq))
746 goto out;
748 isl_seq_clr(eq.data, 1+total);
749 index[isl_seq_get_hash_bits(bmap->div[k], 2+total, bits)] = k + 1;
750 for (--k; k >= 0; --k) {
751 uint32_t hash;
753 if (isl_int_is_zero(bmap->div[k][0]))
754 continue;
756 hash = isl_seq_get_hash_bits(bmap->div[k], 2+total, bits);
757 for (h = hash; index[h]; h = (h+1) % size)
758 if (isl_seq_eq(bmap->div[k],
759 bmap->div[index[h]-1], 2+total))
760 break;
761 if (index[h]) {
762 *progress = 1;
763 l = index[h] - 1;
764 elim_for[l] = k + 1;
766 index[h] = k+1;
768 for (l = bmap->n_div - 1; l >= 0; --l) {
769 if (!elim_for[l])
770 continue;
771 k = elim_for[l] - 1;
772 isl_int_set_si(eq.data[1+total_var+k], -1);
773 isl_int_set_si(eq.data[1+total_var+l], 1);
774 eliminate_div(bmap, eq.data, l, 1);
775 isl_int_set_si(eq.data[1+total_var+k], 0);
776 isl_int_set_si(eq.data[1+total_var+l], 0);
779 isl_blk_free(ctx, eq);
780 out:
781 free(index);
782 free(elim_for);
783 return bmap;
786 static int n_pure_div_eq(struct isl_basic_map *bmap)
788 int i, j;
789 unsigned total;
791 total = isl_space_dim(bmap->dim, isl_dim_all);
792 for (i = 0, j = bmap->n_div-1; i < bmap->n_eq; ++i) {
793 while (j >= 0 && isl_int_is_zero(bmap->eq[i][1 + total + j]))
794 --j;
795 if (j < 0)
796 break;
797 if (isl_seq_first_non_zero(bmap->eq[i] + 1 + total, j) != -1)
798 return 0;
800 return i;
803 /* Normalize divs that appear in equalities.
805 * In particular, we assume that bmap contains some equalities
806 * of the form
808 * a x = m * e_i
810 * and we want to replace the set of e_i by a minimal set and
811 * such that the new e_i have a canonical representation in terms
812 * of the vector x.
813 * If any of the equalities involves more than one divs, then
814 * we currently simply bail out.
816 * Let us first additionally assume that all equalities involve
817 * a div. The equalities then express modulo constraints on the
818 * remaining variables and we can use "parameter compression"
819 * to find a minimal set of constraints. The result is a transformation
821 * x = T(x') = x_0 + G x'
823 * with G a lower-triangular matrix with all elements below the diagonal
824 * non-negative and smaller than the diagonal element on the same row.
825 * We first normalize x_0 by making the same property hold in the affine
826 * T matrix.
827 * The rows i of G with a 1 on the diagonal do not impose any modulo
828 * constraint and simply express x_i = x'_i.
829 * For each of the remaining rows i, we introduce a div and a corresponding
830 * equality. In particular
832 * g_ii e_j = x_i - g_i(x')
834 * where each x'_k is replaced either by x_k (if g_kk = 1) or the
835 * corresponding div (if g_kk != 1).
837 * If there are any equalities not involving any div, then we
838 * first apply a variable compression on the variables x:
840 * x = C x'' x'' = C_2 x
842 * and perform the above parameter compression on A C instead of on A.
843 * The resulting compression is then of the form
845 * x'' = T(x') = x_0 + G x'
847 * and in constructing the new divs and the corresponding equalities,
848 * we have to replace each x'', i.e., the x'_k with (g_kk = 1),
849 * by the corresponding row from C_2.
851 static struct isl_basic_map *normalize_divs(
852 struct isl_basic_map *bmap, int *progress)
854 int i, j, k;
855 int total;
856 int div_eq;
857 struct isl_mat *B;
858 struct isl_vec *d;
859 struct isl_mat *T = NULL;
860 struct isl_mat *C = NULL;
861 struct isl_mat *C2 = NULL;
862 isl_int v;
863 int *pos;
864 int dropped, needed;
866 if (!bmap)
867 return NULL;
869 if (bmap->n_div == 0)
870 return bmap;
872 if (bmap->n_eq == 0)
873 return bmap;
875 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_NORMALIZED_DIVS))
876 return bmap;
878 total = isl_space_dim(bmap->dim, isl_dim_all);
879 div_eq = n_pure_div_eq(bmap);
880 if (div_eq == 0)
881 return bmap;
883 if (div_eq < bmap->n_eq) {
884 B = isl_mat_sub_alloc6(bmap->ctx, bmap->eq, div_eq,
885 bmap->n_eq - div_eq, 0, 1 + total);
886 C = isl_mat_variable_compression(B, &C2);
887 if (!C || !C2)
888 goto error;
889 if (C->n_col == 0) {
890 bmap = isl_basic_map_set_to_empty(bmap);
891 isl_mat_free(C);
892 isl_mat_free(C2);
893 goto done;
897 d = isl_vec_alloc(bmap->ctx, div_eq);
898 if (!d)
899 goto error;
900 for (i = 0, j = bmap->n_div-1; i < div_eq; ++i) {
901 while (j >= 0 && isl_int_is_zero(bmap->eq[i][1 + total + j]))
902 --j;
903 isl_int_set(d->block.data[i], bmap->eq[i][1 + total + j]);
905 B = isl_mat_sub_alloc6(bmap->ctx, bmap->eq, 0, div_eq, 0, 1 + total);
907 if (C) {
908 B = isl_mat_product(B, C);
909 C = NULL;
912 T = isl_mat_parameter_compression(B, d);
913 if (!T)
914 goto error;
915 if (T->n_col == 0) {
916 bmap = isl_basic_map_set_to_empty(bmap);
917 isl_mat_free(C2);
918 isl_mat_free(T);
919 goto done;
921 isl_int_init(v);
922 for (i = 0; i < T->n_row - 1; ++i) {
923 isl_int_fdiv_q(v, T->row[1 + i][0], T->row[1 + i][1 + i]);
924 if (isl_int_is_zero(v))
925 continue;
926 isl_mat_col_submul(T, 0, v, 1 + i);
928 isl_int_clear(v);
929 pos = isl_alloc_array(bmap->ctx, int, T->n_row);
930 if (!pos)
931 goto error;
932 /* We have to be careful because dropping equalities may reorder them */
933 dropped = 0;
934 for (j = bmap->n_div - 1; j >= 0; --j) {
935 for (i = 0; i < bmap->n_eq; ++i)
936 if (!isl_int_is_zero(bmap->eq[i][1 + total + j]))
937 break;
938 if (i < bmap->n_eq) {
939 bmap = isl_basic_map_drop_div(bmap, j);
940 isl_basic_map_drop_equality(bmap, i);
941 ++dropped;
944 pos[0] = 0;
945 needed = 0;
946 for (i = 1; i < T->n_row; ++i) {
947 if (isl_int_is_one(T->row[i][i]))
948 pos[i] = i;
949 else
950 needed++;
952 if (needed > dropped) {
953 bmap = isl_basic_map_extend_space(bmap, isl_space_copy(bmap->dim),
954 needed, needed, 0);
955 if (!bmap)
956 goto error;
958 for (i = 1; i < T->n_row; ++i) {
959 if (isl_int_is_one(T->row[i][i]))
960 continue;
961 k = isl_basic_map_alloc_div(bmap);
962 pos[i] = 1 + total + k;
963 isl_seq_clr(bmap->div[k] + 1, 1 + total + bmap->n_div);
964 isl_int_set(bmap->div[k][0], T->row[i][i]);
965 if (C2)
966 isl_seq_cpy(bmap->div[k] + 1, C2->row[i], 1 + total);
967 else
968 isl_int_set_si(bmap->div[k][1 + i], 1);
969 for (j = 0; j < i; ++j) {
970 if (isl_int_is_zero(T->row[i][j]))
971 continue;
972 if (pos[j] < T->n_row && C2)
973 isl_seq_submul(bmap->div[k] + 1, T->row[i][j],
974 C2->row[pos[j]], 1 + total);
975 else
976 isl_int_neg(bmap->div[k][1 + pos[j]],
977 T->row[i][j]);
979 j = isl_basic_map_alloc_equality(bmap);
980 isl_seq_neg(bmap->eq[j], bmap->div[k]+1, 1+total+bmap->n_div);
981 isl_int_set(bmap->eq[j][pos[i]], bmap->div[k][0]);
983 free(pos);
984 isl_mat_free(C2);
985 isl_mat_free(T);
987 if (progress)
988 *progress = 1;
989 done:
990 ISL_F_SET(bmap, ISL_BASIC_MAP_NORMALIZED_DIVS);
992 return bmap;
993 error:
994 isl_mat_free(C);
995 isl_mat_free(C2);
996 isl_mat_free(T);
997 return bmap;
1000 static struct isl_basic_map *set_div_from_lower_bound(
1001 struct isl_basic_map *bmap, int div, int ineq)
1003 unsigned total = 1 + isl_space_dim(bmap->dim, isl_dim_all);
1005 isl_seq_neg(bmap->div[div] + 1, bmap->ineq[ineq], total + bmap->n_div);
1006 isl_int_set(bmap->div[div][0], bmap->ineq[ineq][total + div]);
1007 isl_int_add(bmap->div[div][1], bmap->div[div][1], bmap->div[div][0]);
1008 isl_int_sub_ui(bmap->div[div][1], bmap->div[div][1], 1);
1009 isl_int_set_si(bmap->div[div][1 + total + div], 0);
1011 return bmap;
1014 /* Check whether it is ok to define a div based on an inequality.
1015 * To avoid the introduction of circular definitions of divs, we
1016 * do not allow such a definition if the resulting expression would refer to
1017 * any other undefined divs or if any known div is defined in
1018 * terms of the unknown div.
1020 static int ok_to_set_div_from_bound(struct isl_basic_map *bmap,
1021 int div, int ineq)
1023 int j;
1024 unsigned total = 1 + isl_space_dim(bmap->dim, isl_dim_all);
1026 /* Not defined in terms of unknown divs */
1027 for (j = 0; j < bmap->n_div; ++j) {
1028 if (div == j)
1029 continue;
1030 if (isl_int_is_zero(bmap->ineq[ineq][total + j]))
1031 continue;
1032 if (isl_int_is_zero(bmap->div[j][0]))
1033 return 0;
1036 /* No other div defined in terms of this one => avoid loops */
1037 for (j = 0; j < bmap->n_div; ++j) {
1038 if (div == j)
1039 continue;
1040 if (isl_int_is_zero(bmap->div[j][0]))
1041 continue;
1042 if (!isl_int_is_zero(bmap->div[j][1 + total + div]))
1043 return 0;
1046 return 1;
1049 /* Would an expression for div "div" based on inequality "ineq" of "bmap"
1050 * be a better expression than the current one?
1052 * If we do not have any expression yet, then any expression would be better.
1053 * Otherwise we check if the last variable involved in the inequality
1054 * (disregarding the div that it would define) is in an earlier position
1055 * than the last variable involved in the current div expression.
1057 static int better_div_constraint(__isl_keep isl_basic_map *bmap,
1058 int div, int ineq)
1060 unsigned total = 1 + isl_space_dim(bmap->dim, isl_dim_all);
1061 int last_div;
1062 int last_ineq;
1064 if (isl_int_is_zero(bmap->div[div][0]))
1065 return 1;
1067 if (isl_seq_last_non_zero(bmap->ineq[ineq] + total + div + 1,
1068 bmap->n_div - (div + 1)) >= 0)
1069 return 0;
1071 last_ineq = isl_seq_last_non_zero(bmap->ineq[ineq], total + div);
1072 last_div = isl_seq_last_non_zero(bmap->div[div] + 1,
1073 total + bmap->n_div);
1075 return last_ineq < last_div;
1078 /* Given two constraints "k" and "l" that are opposite to each other,
1079 * except for the constant term, check if we can use them
1080 * to obtain an expression for one of the hitherto unknown divs or
1081 * a "better" expression for a div for which we already have an expression.
1082 * "sum" is the sum of the constant terms of the constraints.
1083 * If this sum is strictly smaller than the coefficient of one
1084 * of the divs, then this pair can be used define the div.
1085 * To avoid the introduction of circular definitions of divs, we
1086 * do not use the pair if the resulting expression would refer to
1087 * any other undefined divs or if any known div is defined in
1088 * terms of the unknown div.
1090 static struct isl_basic_map *check_for_div_constraints(
1091 struct isl_basic_map *bmap, int k, int l, isl_int sum, int *progress)
1093 int i;
1094 unsigned total = 1 + isl_space_dim(bmap->dim, isl_dim_all);
1096 for (i = 0; i < bmap->n_div; ++i) {
1097 if (isl_int_is_zero(bmap->ineq[k][total + i]))
1098 continue;
1099 if (isl_int_abs_ge(sum, bmap->ineq[k][total + i]))
1100 continue;
1101 if (!better_div_constraint(bmap, i, k))
1102 continue;
1103 if (!ok_to_set_div_from_bound(bmap, i, k))
1104 break;
1105 if (isl_int_is_pos(bmap->ineq[k][total + i]))
1106 bmap = set_div_from_lower_bound(bmap, i, k);
1107 else
1108 bmap = set_div_from_lower_bound(bmap, i, l);
1109 if (progress)
1110 *progress = 1;
1111 break;
1113 return bmap;
1116 __isl_give isl_basic_map *isl_basic_map_remove_duplicate_constraints(
1117 __isl_take isl_basic_map *bmap, int *progress, int detect_divs)
1119 unsigned int size;
1120 isl_int ***index;
1121 int k, l, h;
1122 int bits;
1123 unsigned total = isl_basic_map_total_dim(bmap);
1124 isl_int sum;
1125 isl_ctx *ctx;
1127 if (!bmap || bmap->n_ineq <= 1)
1128 return bmap;
1130 size = round_up(4 * (bmap->n_ineq+1) / 3 - 1);
1131 bits = ffs(size) - 1;
1132 ctx = isl_basic_map_get_ctx(bmap);
1133 index = isl_calloc_array(ctx, isl_int **, size);
1134 if (!index)
1135 return bmap;
1137 index[isl_seq_get_hash_bits(bmap->ineq[0]+1, total, bits)] = &bmap->ineq[0];
1138 for (k = 1; k < bmap->n_ineq; ++k) {
1139 h = hash_index(index, size, bits, bmap, k);
1140 if (!index[h]) {
1141 index[h] = &bmap->ineq[k];
1142 continue;
1144 if (progress)
1145 *progress = 1;
1146 l = index[h] - &bmap->ineq[0];
1147 if (isl_int_lt(bmap->ineq[k][0], bmap->ineq[l][0]))
1148 swap_inequality(bmap, k, l);
1149 isl_basic_map_drop_inequality(bmap, k);
1150 --k;
1152 isl_int_init(sum);
1153 for (k = 0; k < bmap->n_ineq-1; ++k) {
1154 isl_seq_neg(bmap->ineq[k]+1, bmap->ineq[k]+1, total);
1155 h = hash_index(index, size, bits, bmap, k);
1156 isl_seq_neg(bmap->ineq[k]+1, bmap->ineq[k]+1, total);
1157 if (!index[h])
1158 continue;
1159 l = index[h] - &bmap->ineq[0];
1160 isl_int_add(sum, bmap->ineq[k][0], bmap->ineq[l][0]);
1161 if (isl_int_is_pos(sum)) {
1162 if (detect_divs)
1163 bmap = check_for_div_constraints(bmap, k, l,
1164 sum, progress);
1165 continue;
1167 if (isl_int_is_zero(sum)) {
1168 /* We need to break out of the loop after these
1169 * changes since the contents of the hash
1170 * will no longer be valid.
1171 * Plus, we probably we want to regauss first.
1173 if (progress)
1174 *progress = 1;
1175 isl_basic_map_drop_inequality(bmap, l);
1176 isl_basic_map_inequality_to_equality(bmap, k);
1177 } else
1178 bmap = isl_basic_map_set_to_empty(bmap);
1179 break;
1181 isl_int_clear(sum);
1183 free(index);
1184 return bmap;
1188 /* Eliminate knowns divs from constraints where they appear with
1189 * a (positive or negative) unit coefficient.
1191 * That is, replace
1193 * floor(e/m) + f >= 0
1195 * by
1197 * e + m f >= 0
1199 * and
1201 * -floor(e/m) + f >= 0
1203 * by
1205 * -e + m f + m - 1 >= 0
1207 * The first conversion is valid because floor(e/m) >= -f is equivalent
1208 * to e/m >= -f because -f is an integral expression.
1209 * The second conversion follows from the fact that
1211 * -floor(e/m) = ceil(-e/m) = floor((-e + m - 1)/m)
1214 * Note that one of the div constraints may have been eliminated
1215 * due to being redundant with respect to the constraint that is
1216 * being modified by this function. The modified constraint may
1217 * no longer imply this div constraint, so we add it back to make
1218 * sure we do not lose any information.
1220 * We skip integral divs, i.e., those with denominator 1, as we would
1221 * risk eliminating the div from the div constraints. We do not need
1222 * to handle those divs here anyway since the div constraints will turn
1223 * out to form an equality and this equality can then be use to eliminate
1224 * the div from all constraints.
1226 static __isl_give isl_basic_map *eliminate_unit_divs(
1227 __isl_take isl_basic_map *bmap, int *progress)
1229 int i, j;
1230 isl_ctx *ctx;
1231 unsigned total;
1233 if (!bmap)
1234 return NULL;
1236 ctx = isl_basic_map_get_ctx(bmap);
1237 total = 1 + isl_space_dim(bmap->dim, isl_dim_all);
1239 for (i = 0; i < bmap->n_div; ++i) {
1240 if (isl_int_is_zero(bmap->div[i][0]))
1241 continue;
1242 if (isl_int_is_one(bmap->div[i][0]))
1243 continue;
1244 for (j = 0; j < bmap->n_ineq; ++j) {
1245 int s;
1247 if (!isl_int_is_one(bmap->ineq[j][total + i]) &&
1248 !isl_int_is_negone(bmap->ineq[j][total + i]))
1249 continue;
1251 *progress = 1;
1253 s = isl_int_sgn(bmap->ineq[j][total + i]);
1254 isl_int_set_si(bmap->ineq[j][total + i], 0);
1255 if (s < 0)
1256 isl_seq_combine(bmap->ineq[j],
1257 ctx->negone, bmap->div[i] + 1,
1258 bmap->div[i][0], bmap->ineq[j],
1259 total + bmap->n_div);
1260 else
1261 isl_seq_combine(bmap->ineq[j],
1262 ctx->one, bmap->div[i] + 1,
1263 bmap->div[i][0], bmap->ineq[j],
1264 total + bmap->n_div);
1265 if (s < 0) {
1266 isl_int_add(bmap->ineq[j][0],
1267 bmap->ineq[j][0], bmap->div[i][0]);
1268 isl_int_sub_ui(bmap->ineq[j][0],
1269 bmap->ineq[j][0], 1);
1272 bmap = isl_basic_map_extend_constraints(bmap, 0, 1);
1273 if (isl_basic_map_add_div_constraint(bmap, i, s) < 0)
1274 return isl_basic_map_free(bmap);
1278 return bmap;
1281 struct isl_basic_map *isl_basic_map_simplify(struct isl_basic_map *bmap)
1283 int progress = 1;
1284 if (!bmap)
1285 return NULL;
1286 while (progress) {
1287 progress = 0;
1288 if (!bmap)
1289 break;
1290 if (isl_basic_map_plain_is_empty(bmap))
1291 break;
1292 bmap = isl_basic_map_normalize_constraints(bmap);
1293 bmap = normalize_div_expressions(bmap);
1294 bmap = remove_duplicate_divs(bmap, &progress);
1295 bmap = eliminate_unit_divs(bmap, &progress);
1296 bmap = eliminate_divs_eq(bmap, &progress);
1297 bmap = eliminate_divs_ineq(bmap, &progress);
1298 bmap = isl_basic_map_gauss(bmap, &progress);
1299 /* requires equalities in normal form */
1300 bmap = normalize_divs(bmap, &progress);
1301 bmap = isl_basic_map_remove_duplicate_constraints(bmap,
1302 &progress, 1);
1304 return bmap;
1307 struct isl_basic_set *isl_basic_set_simplify(struct isl_basic_set *bset)
1309 return (struct isl_basic_set *)
1310 isl_basic_map_simplify((struct isl_basic_map *)bset);
1314 int isl_basic_map_is_div_constraint(__isl_keep isl_basic_map *bmap,
1315 isl_int *constraint, unsigned div)
1317 unsigned pos;
1319 if (!bmap)
1320 return -1;
1322 pos = 1 + isl_space_dim(bmap->dim, isl_dim_all) + div;
1324 if (isl_int_eq(constraint[pos], bmap->div[div][0])) {
1325 int neg;
1326 isl_int_sub(bmap->div[div][1],
1327 bmap->div[div][1], bmap->div[div][0]);
1328 isl_int_add_ui(bmap->div[div][1], bmap->div[div][1], 1);
1329 neg = isl_seq_is_neg(constraint, bmap->div[div]+1, pos);
1330 isl_int_sub_ui(bmap->div[div][1], bmap->div[div][1], 1);
1331 isl_int_add(bmap->div[div][1],
1332 bmap->div[div][1], bmap->div[div][0]);
1333 if (!neg)
1334 return 0;
1335 if (isl_seq_first_non_zero(constraint+pos+1,
1336 bmap->n_div-div-1) != -1)
1337 return 0;
1338 } else if (isl_int_abs_eq(constraint[pos], bmap->div[div][0])) {
1339 if (!isl_seq_eq(constraint, bmap->div[div]+1, pos))
1340 return 0;
1341 if (isl_seq_first_non_zero(constraint+pos+1,
1342 bmap->n_div-div-1) != -1)
1343 return 0;
1344 } else
1345 return 0;
1347 return 1;
1350 int isl_basic_set_is_div_constraint(__isl_keep isl_basic_set *bset,
1351 isl_int *constraint, unsigned div)
1353 return isl_basic_map_is_div_constraint(bset, constraint, div);
1357 /* If the only constraints a div d=floor(f/m)
1358 * appears in are its two defining constraints
1360 * f - m d >=0
1361 * -(f - (m - 1)) + m d >= 0
1363 * then it can safely be removed.
1365 static int div_is_redundant(struct isl_basic_map *bmap, int div)
1367 int i;
1368 unsigned pos = 1 + isl_space_dim(bmap->dim, isl_dim_all) + div;
1370 for (i = 0; i < bmap->n_eq; ++i)
1371 if (!isl_int_is_zero(bmap->eq[i][pos]))
1372 return 0;
1374 for (i = 0; i < bmap->n_ineq; ++i) {
1375 if (isl_int_is_zero(bmap->ineq[i][pos]))
1376 continue;
1377 if (!isl_basic_map_is_div_constraint(bmap, bmap->ineq[i], div))
1378 return 0;
1381 for (i = 0; i < bmap->n_div; ++i) {
1382 if (isl_int_is_zero(bmap->div[i][0]))
1383 continue;
1384 if (!isl_int_is_zero(bmap->div[i][1+pos]))
1385 return 0;
1388 return 1;
1392 * Remove divs that don't occur in any of the constraints or other divs.
1393 * These can arise when dropping constraints from a basic map or
1394 * when the divs of a basic map have been temporarily aligned
1395 * with the divs of another basic map.
1397 static struct isl_basic_map *remove_redundant_divs(struct isl_basic_map *bmap)
1399 int i;
1401 if (!bmap)
1402 return NULL;
1404 for (i = bmap->n_div-1; i >= 0; --i) {
1405 if (!div_is_redundant(bmap, i))
1406 continue;
1407 bmap = isl_basic_map_drop_div(bmap, i);
1409 return bmap;
1412 struct isl_basic_map *isl_basic_map_finalize(struct isl_basic_map *bmap)
1414 bmap = remove_redundant_divs(bmap);
1415 if (!bmap)
1416 return NULL;
1417 ISL_F_SET(bmap, ISL_BASIC_SET_FINAL);
1418 return bmap;
1421 struct isl_basic_set *isl_basic_set_finalize(struct isl_basic_set *bset)
1423 return (struct isl_basic_set *)
1424 isl_basic_map_finalize((struct isl_basic_map *)bset);
1427 struct isl_set *isl_set_finalize(struct isl_set *set)
1429 int i;
1431 if (!set)
1432 return NULL;
1433 for (i = 0; i < set->n; ++i) {
1434 set->p[i] = isl_basic_set_finalize(set->p[i]);
1435 if (!set->p[i])
1436 goto error;
1438 return set;
1439 error:
1440 isl_set_free(set);
1441 return NULL;
1444 struct isl_map *isl_map_finalize(struct isl_map *map)
1446 int i;
1448 if (!map)
1449 return NULL;
1450 for (i = 0; i < map->n; ++i) {
1451 map->p[i] = isl_basic_map_finalize(map->p[i]);
1452 if (!map->p[i])
1453 goto error;
1455 ISL_F_CLR(map, ISL_MAP_NORMALIZED);
1456 return map;
1457 error:
1458 isl_map_free(map);
1459 return NULL;
1463 /* Remove definition of any div that is defined in terms of the given variable.
1464 * The div itself is not removed. Functions such as
1465 * eliminate_divs_ineq depend on the other divs remaining in place.
1467 static struct isl_basic_map *remove_dependent_vars(struct isl_basic_map *bmap,
1468 int pos)
1470 int i;
1472 if (!bmap)
1473 return NULL;
1475 for (i = 0; i < bmap->n_div; ++i) {
1476 if (isl_int_is_zero(bmap->div[i][0]))
1477 continue;
1478 if (isl_int_is_zero(bmap->div[i][1+1+pos]))
1479 continue;
1480 isl_int_set_si(bmap->div[i][0], 0);
1482 return bmap;
1485 /* Eliminate the specified variables from the constraints using
1486 * Fourier-Motzkin. The variables themselves are not removed.
1488 struct isl_basic_map *isl_basic_map_eliminate_vars(
1489 struct isl_basic_map *bmap, unsigned pos, unsigned n)
1491 int d;
1492 int i, j, k;
1493 unsigned total;
1494 int need_gauss = 0;
1496 if (n == 0)
1497 return bmap;
1498 if (!bmap)
1499 return NULL;
1500 total = isl_basic_map_total_dim(bmap);
1502 bmap = isl_basic_map_cow(bmap);
1503 for (d = pos + n - 1; d >= 0 && d >= pos; --d)
1504 bmap = remove_dependent_vars(bmap, d);
1505 if (!bmap)
1506 return NULL;
1508 for (d = pos + n - 1;
1509 d >= 0 && d >= total - bmap->n_div && d >= pos; --d)
1510 isl_seq_clr(bmap->div[d-(total-bmap->n_div)], 2+total);
1511 for (d = pos + n - 1; d >= 0 && d >= pos; --d) {
1512 int n_lower, n_upper;
1513 if (!bmap)
1514 return NULL;
1515 for (i = 0; i < bmap->n_eq; ++i) {
1516 if (isl_int_is_zero(bmap->eq[i][1+d]))
1517 continue;
1518 eliminate_var_using_equality(bmap, d, bmap->eq[i], 0, NULL);
1519 isl_basic_map_drop_equality(bmap, i);
1520 need_gauss = 1;
1521 break;
1523 if (i < bmap->n_eq)
1524 continue;
1525 n_lower = 0;
1526 n_upper = 0;
1527 for (i = 0; i < bmap->n_ineq; ++i) {
1528 if (isl_int_is_pos(bmap->ineq[i][1+d]))
1529 n_lower++;
1530 else if (isl_int_is_neg(bmap->ineq[i][1+d]))
1531 n_upper++;
1533 bmap = isl_basic_map_extend_constraints(bmap,
1534 0, n_lower * n_upper);
1535 if (!bmap)
1536 goto error;
1537 for (i = bmap->n_ineq - 1; i >= 0; --i) {
1538 int last;
1539 if (isl_int_is_zero(bmap->ineq[i][1+d]))
1540 continue;
1541 last = -1;
1542 for (j = 0; j < i; ++j) {
1543 if (isl_int_is_zero(bmap->ineq[j][1+d]))
1544 continue;
1545 last = j;
1546 if (isl_int_sgn(bmap->ineq[i][1+d]) ==
1547 isl_int_sgn(bmap->ineq[j][1+d]))
1548 continue;
1549 k = isl_basic_map_alloc_inequality(bmap);
1550 if (k < 0)
1551 goto error;
1552 isl_seq_cpy(bmap->ineq[k], bmap->ineq[i],
1553 1+total);
1554 isl_seq_elim(bmap->ineq[k], bmap->ineq[j],
1555 1+d, 1+total, NULL);
1557 isl_basic_map_drop_inequality(bmap, i);
1558 i = last + 1;
1560 if (n_lower > 0 && n_upper > 0) {
1561 bmap = isl_basic_map_normalize_constraints(bmap);
1562 bmap = isl_basic_map_remove_duplicate_constraints(bmap,
1563 NULL, 0);
1564 bmap = isl_basic_map_gauss(bmap, NULL);
1565 bmap = isl_basic_map_remove_redundancies(bmap);
1566 need_gauss = 0;
1567 if (!bmap)
1568 goto error;
1569 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
1570 break;
1573 ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
1574 if (need_gauss)
1575 bmap = isl_basic_map_gauss(bmap, NULL);
1576 return bmap;
1577 error:
1578 isl_basic_map_free(bmap);
1579 return NULL;
1582 struct isl_basic_set *isl_basic_set_eliminate_vars(
1583 struct isl_basic_set *bset, unsigned pos, unsigned n)
1585 return (struct isl_basic_set *)isl_basic_map_eliminate_vars(
1586 (struct isl_basic_map *)bset, pos, n);
1589 /* Eliminate the specified n dimensions starting at first from the
1590 * constraints, without removing the dimensions from the space.
1591 * If the set is rational, the dimensions are eliminated using Fourier-Motzkin.
1592 * Otherwise, they are projected out and the original space is restored.
1594 __isl_give isl_basic_map *isl_basic_map_eliminate(
1595 __isl_take isl_basic_map *bmap,
1596 enum isl_dim_type type, unsigned first, unsigned n)
1598 isl_space *space;
1600 if (!bmap)
1601 return NULL;
1602 if (n == 0)
1603 return bmap;
1605 if (first + n > isl_basic_map_dim(bmap, type) || first + n < first)
1606 isl_die(bmap->ctx, isl_error_invalid,
1607 "index out of bounds", goto error);
1609 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL)) {
1610 first += isl_basic_map_offset(bmap, type) - 1;
1611 bmap = isl_basic_map_eliminate_vars(bmap, first, n);
1612 return isl_basic_map_finalize(bmap);
1615 space = isl_basic_map_get_space(bmap);
1616 bmap = isl_basic_map_project_out(bmap, type, first, n);
1617 bmap = isl_basic_map_insert_dims(bmap, type, first, n);
1618 bmap = isl_basic_map_reset_space(bmap, space);
1619 return bmap;
1620 error:
1621 isl_basic_map_free(bmap);
1622 return NULL;
1625 __isl_give isl_basic_set *isl_basic_set_eliminate(
1626 __isl_take isl_basic_set *bset,
1627 enum isl_dim_type type, unsigned first, unsigned n)
1629 return isl_basic_map_eliminate(bset, type, first, n);
1632 /* Don't assume equalities are in order, because align_divs
1633 * may have changed the order of the divs.
1635 static void compute_elimination_index(struct isl_basic_map *bmap, int *elim)
1637 int d, i;
1638 unsigned total;
1640 total = isl_space_dim(bmap->dim, isl_dim_all);
1641 for (d = 0; d < total; ++d)
1642 elim[d] = -1;
1643 for (i = 0; i < bmap->n_eq; ++i) {
1644 for (d = total - 1; d >= 0; --d) {
1645 if (isl_int_is_zero(bmap->eq[i][1+d]))
1646 continue;
1647 elim[d] = i;
1648 break;
1653 static void set_compute_elimination_index(struct isl_basic_set *bset, int *elim)
1655 compute_elimination_index((struct isl_basic_map *)bset, elim);
1658 static int reduced_using_equalities(isl_int *dst, isl_int *src,
1659 struct isl_basic_map *bmap, int *elim)
1661 int d;
1662 int copied = 0;
1663 unsigned total;
1665 total = isl_space_dim(bmap->dim, isl_dim_all);
1666 for (d = total - 1; d >= 0; --d) {
1667 if (isl_int_is_zero(src[1+d]))
1668 continue;
1669 if (elim[d] == -1)
1670 continue;
1671 if (!copied) {
1672 isl_seq_cpy(dst, src, 1 + total);
1673 copied = 1;
1675 isl_seq_elim(dst, bmap->eq[elim[d]], 1 + d, 1 + total, NULL);
1677 return copied;
1680 static int set_reduced_using_equalities(isl_int *dst, isl_int *src,
1681 struct isl_basic_set *bset, int *elim)
1683 return reduced_using_equalities(dst, src,
1684 (struct isl_basic_map *)bset, elim);
1687 static struct isl_basic_set *isl_basic_set_reduce_using_equalities(
1688 struct isl_basic_set *bset, struct isl_basic_set *context)
1690 int i;
1691 int *elim;
1693 if (!bset || !context)
1694 goto error;
1696 if (context->n_eq == 0) {
1697 isl_basic_set_free(context);
1698 return bset;
1701 bset = isl_basic_set_cow(bset);
1702 if (!bset)
1703 goto error;
1705 elim = isl_alloc_array(bset->ctx, int, isl_basic_set_n_dim(bset));
1706 if (!elim)
1707 goto error;
1708 set_compute_elimination_index(context, elim);
1709 for (i = 0; i < bset->n_eq; ++i)
1710 set_reduced_using_equalities(bset->eq[i], bset->eq[i],
1711 context, elim);
1712 for (i = 0; i < bset->n_ineq; ++i)
1713 set_reduced_using_equalities(bset->ineq[i], bset->ineq[i],
1714 context, elim);
1715 isl_basic_set_free(context);
1716 free(elim);
1717 bset = isl_basic_set_simplify(bset);
1718 bset = isl_basic_set_finalize(bset);
1719 return bset;
1720 error:
1721 isl_basic_set_free(bset);
1722 isl_basic_set_free(context);
1723 return NULL;
1726 static struct isl_basic_set *remove_shifted_constraints(
1727 struct isl_basic_set *bset, struct isl_basic_set *context)
1729 unsigned int size;
1730 isl_int ***index;
1731 int bits;
1732 int k, h, l;
1733 isl_ctx *ctx;
1735 if (!bset)
1736 return NULL;
1738 size = round_up(4 * (context->n_ineq+1) / 3 - 1);
1739 bits = ffs(size) - 1;
1740 ctx = isl_basic_set_get_ctx(bset);
1741 index = isl_calloc_array(ctx, isl_int **, size);
1742 if (!index)
1743 return bset;
1745 for (k = 0; k < context->n_ineq; ++k) {
1746 h = set_hash_index(index, size, bits, context, k);
1747 index[h] = &context->ineq[k];
1749 for (k = 0; k < bset->n_ineq; ++k) {
1750 h = set_hash_index(index, size, bits, bset, k);
1751 if (!index[h])
1752 continue;
1753 l = index[h] - &context->ineq[0];
1754 if (isl_int_lt(bset->ineq[k][0], context->ineq[l][0]))
1755 continue;
1756 bset = isl_basic_set_cow(bset);
1757 if (!bset)
1758 goto error;
1759 isl_basic_set_drop_inequality(bset, k);
1760 --k;
1762 free(index);
1763 return bset;
1764 error:
1765 free(index);
1766 return bset;
1769 /* Remove constraints from "bmap" that are identical to constraints
1770 * in "context" or that are more relaxed (greater constant term).
1772 * We perform the test for shifted copies on the pure constraints
1773 * in remove_shifted_constraints.
1775 static __isl_give isl_basic_map *isl_basic_map_remove_shifted_constraints(
1776 __isl_take isl_basic_map *bmap, __isl_take isl_basic_map *context)
1778 isl_basic_set *bset, *bset_context;
1780 if (!bmap || !context)
1781 goto error;
1783 if (bmap->n_ineq == 0 || context->n_ineq == 0) {
1784 isl_basic_map_free(context);
1785 return bmap;
1788 context = isl_basic_map_align_divs(context, bmap);
1789 bmap = isl_basic_map_align_divs(bmap, context);
1791 bset = isl_basic_map_underlying_set(isl_basic_map_copy(bmap));
1792 bset_context = isl_basic_map_underlying_set(context);
1793 bset = remove_shifted_constraints(bset, bset_context);
1794 isl_basic_set_free(bset_context);
1796 bmap = isl_basic_map_overlying_set(bset, bmap);
1798 return bmap;
1799 error:
1800 isl_basic_map_free(bmap);
1801 isl_basic_map_free(context);
1802 return NULL;
1805 /* Does the (linear part of a) constraint "c" involve any of the "len"
1806 * "relevant" dimensions?
1808 static int is_related(isl_int *c, int len, int *relevant)
1810 int i;
1812 for (i = 0; i < len; ++i) {
1813 if (!relevant[i])
1814 continue;
1815 if (!isl_int_is_zero(c[i]))
1816 return 1;
1819 return 0;
1822 /* Drop constraints from "bset" that do not involve any of
1823 * the dimensions marked "relevant".
1825 static __isl_give isl_basic_set *drop_unrelated_constraints(
1826 __isl_take isl_basic_set *bset, int *relevant)
1828 int i, dim;
1830 dim = isl_basic_set_dim(bset, isl_dim_set);
1831 for (i = 0; i < dim; ++i)
1832 if (!relevant[i])
1833 break;
1834 if (i >= dim)
1835 return bset;
1837 for (i = bset->n_eq - 1; i >= 0; --i)
1838 if (!is_related(bset->eq[i] + 1, dim, relevant))
1839 isl_basic_set_drop_equality(bset, i);
1841 for (i = bset->n_ineq - 1; i >= 0; --i)
1842 if (!is_related(bset->ineq[i] + 1, dim, relevant))
1843 isl_basic_set_drop_inequality(bset, i);
1845 return bset;
1848 /* Update the groups in "group" based on the (linear part of a) constraint "c".
1850 * In particular, for any variable involved in the constraint,
1851 * find the actual group id from before and replace the group
1852 * of the corresponding variable by the minimal group of all
1853 * the variables involved in the constraint considered so far
1854 * (if this minimum is smaller) or replace the minimum by this group
1855 * (if the minimum is larger).
1857 * At the end, all the variables in "c" will (indirectly) point
1858 * to the minimal of the groups that they referred to originally.
1860 static void update_groups(int dim, int *group, isl_int *c)
1862 int j;
1863 int min = dim;
1865 for (j = 0; j < dim; ++j) {
1866 if (isl_int_is_zero(c[j]))
1867 continue;
1868 while (group[j] >= 0 && group[group[j]] != group[j])
1869 group[j] = group[group[j]];
1870 if (group[j] == min)
1871 continue;
1872 if (group[j] < min) {
1873 if (min >= 0 && min < dim)
1874 group[min] = group[j];
1875 min = group[j];
1876 } else
1877 group[group[j]] = min;
1881 /* Drop constraints from "context" that are irrelevant for computing
1882 * the gist of "bset".
1884 * In particular, drop constraints in variables that are not related
1885 * to any of the variables involved in the constraints of "bset"
1886 * in the sense that there is no sequence of constraints that connects them.
1888 * We construct groups of variables that collect variables that
1889 * (indirectly) appear in some common constraint of "context".
1890 * Each group is identified by the first variable in the group,
1891 * except for the special group of variables that appear in "bset"
1892 * (or are related to those variables), which is identified by -1.
1893 * If group[i] is equal to i (or -1), then the group of i is i (or -1),
1894 * otherwise the group of i is the group of group[i].
1896 * We first initialize the -1 group with the variables that appear in "bset".
1897 * Then we initialize groups for the remaining variables.
1898 * Then we iterate over the constraints of "context" and update the
1899 * group of the variables in the constraint by the smallest group.
1900 * Finally, we resolve indirect references to groups by running over
1901 * the variables.
1903 * After computing the groups, we drop constraints that do not involve
1904 * any variables in the -1 group.
1906 static __isl_give isl_basic_set *drop_irrelevant_constraints(
1907 __isl_take isl_basic_set *context, __isl_keep isl_basic_set *bset)
1909 isl_ctx *ctx;
1910 int *group;
1911 int dim;
1912 int i, j;
1913 int last;
1915 if (!context || !bset)
1916 return isl_basic_set_free(context);
1918 dim = isl_basic_set_dim(bset, isl_dim_set);
1919 ctx = isl_basic_set_get_ctx(bset);
1920 group = isl_calloc_array(ctx, int, dim);
1922 if (!group)
1923 goto error;
1925 for (i = 0; i < dim; ++i) {
1926 for (j = 0; j < bset->n_eq; ++j)
1927 if (!isl_int_is_zero(bset->eq[j][1 + i]))
1928 break;
1929 if (j < bset->n_eq) {
1930 group[i] = -1;
1931 continue;
1933 for (j = 0; j < bset->n_ineq; ++j)
1934 if (!isl_int_is_zero(bset->ineq[j][1 + i]))
1935 break;
1936 if (j < bset->n_ineq)
1937 group[i] = -1;
1940 last = -1;
1941 for (i = 0; i < dim; ++i)
1942 if (group[i] >= 0)
1943 last = group[i] = i;
1944 if (last < 0) {
1945 free(group);
1946 return context;
1949 for (i = 0; i < context->n_eq; ++i)
1950 update_groups(dim, group, context->eq[i] + 1);
1951 for (i = 0; i < context->n_ineq; ++i)
1952 update_groups(dim, group, context->ineq[i] + 1);
1954 for (i = 0; i < dim; ++i)
1955 if (group[i] >= 0)
1956 group[i] = group[group[i]];
1958 for (i = 0; i < dim; ++i)
1959 group[i] = group[i] == -1;
1961 context = drop_unrelated_constraints(context, group);
1963 free(group);
1964 return context;
1965 error:
1966 free(group);
1967 return isl_basic_set_free(context);
1970 /* Remove all information from bset that is redundant in the context
1971 * of context. Both bset and context are assumed to be full-dimensional.
1973 * We first remove the inequalities from "bset"
1974 * that are obviously redundant with respect to some inequality in "context".
1975 * Then we remove those constraints from "context" that have become
1976 * irrelevant for computing the gist of "bset".
1977 * Note that this removal of constraints cannot be replaced by
1978 * a factorization because factors in "bset" may still be connected
1979 * to each other through constraints in "context".
1981 * If there are any inequalities left, we construct a tableau for
1982 * the context and then add the inequalities of "bset".
1983 * Before adding these inequalities, we freeze all constraints such that
1984 * they won't be considered redundant in terms of the constraints of "bset".
1985 * Then we detect all redundant constraints (among the
1986 * constraints that weren't frozen), first by checking for redundancy in the
1987 * the tableau and then by checking if replacing a constraint by its negation
1988 * would lead to an empty set. This last step is fairly expensive
1989 * and could be optimized by more reuse of the tableau.
1990 * Finally, we update bset according to the results.
1992 static __isl_give isl_basic_set *uset_gist_full(__isl_take isl_basic_set *bset,
1993 __isl_take isl_basic_set *context)
1995 int i, k;
1996 isl_basic_set *combined = NULL;
1997 struct isl_tab *tab = NULL;
1998 unsigned context_ineq;
1999 unsigned total;
2001 if (!bset || !context)
2002 goto error;
2004 if (isl_basic_set_is_universe(bset)) {
2005 isl_basic_set_free(context);
2006 return bset;
2009 if (isl_basic_set_is_universe(context)) {
2010 isl_basic_set_free(context);
2011 return bset;
2014 bset = remove_shifted_constraints(bset, context);
2015 if (!bset)
2016 goto error;
2017 if (bset->n_ineq == 0)
2018 goto done;
2020 context = drop_irrelevant_constraints(context, bset);
2021 if (!context)
2022 goto error;
2023 if (isl_basic_set_is_universe(context)) {
2024 isl_basic_set_free(context);
2025 return bset;
2028 context_ineq = context->n_ineq;
2029 combined = isl_basic_set_cow(isl_basic_set_copy(context));
2030 combined = isl_basic_set_extend_constraints(combined, 0, bset->n_ineq);
2031 tab = isl_tab_from_basic_set(combined, 0);
2032 for (i = 0; i < context_ineq; ++i)
2033 if (isl_tab_freeze_constraint(tab, i) < 0)
2034 goto error;
2035 if (isl_tab_extend_cons(tab, bset->n_ineq) < 0)
2036 goto error;
2037 for (i = 0; i < bset->n_ineq; ++i)
2038 if (isl_tab_add_ineq(tab, bset->ineq[i]) < 0)
2039 goto error;
2040 bset = isl_basic_set_add_constraints(combined, bset, 0);
2041 combined = NULL;
2042 if (!bset)
2043 goto error;
2044 if (isl_tab_detect_redundant(tab) < 0)
2045 goto error;
2046 total = isl_basic_set_total_dim(bset);
2047 for (i = context_ineq; i < bset->n_ineq; ++i) {
2048 int is_empty;
2049 if (tab->con[i].is_redundant)
2050 continue;
2051 tab->con[i].is_redundant = 1;
2052 combined = isl_basic_set_dup(bset);
2053 combined = isl_basic_set_update_from_tab(combined, tab);
2054 combined = isl_basic_set_extend_constraints(combined, 0, 1);
2055 k = isl_basic_set_alloc_inequality(combined);
2056 if (k < 0)
2057 goto error;
2058 isl_seq_neg(combined->ineq[k], bset->ineq[i], 1 + total);
2059 isl_int_sub_ui(combined->ineq[k][0], combined->ineq[k][0], 1);
2060 is_empty = isl_basic_set_is_empty(combined);
2061 if (is_empty < 0)
2062 goto error;
2063 isl_basic_set_free(combined);
2064 combined = NULL;
2065 if (!is_empty)
2066 tab->con[i].is_redundant = 0;
2068 for (i = 0; i < context_ineq; ++i)
2069 tab->con[i].is_redundant = 1;
2070 bset = isl_basic_set_update_from_tab(bset, tab);
2071 if (bset) {
2072 ISL_F_SET(bset, ISL_BASIC_SET_NO_IMPLICIT);
2073 ISL_F_SET(bset, ISL_BASIC_SET_NO_REDUNDANT);
2076 isl_tab_free(tab);
2077 done:
2078 bset = isl_basic_set_simplify(bset);
2079 bset = isl_basic_set_finalize(bset);
2080 isl_basic_set_free(context);
2081 return bset;
2082 error:
2083 isl_tab_free(tab);
2084 isl_basic_set_free(combined);
2085 isl_basic_set_free(context);
2086 isl_basic_set_free(bset);
2087 return NULL;
2090 /* Remove all information from bset that is redundant in the context
2091 * of context. In particular, equalities that are linear combinations
2092 * of those in context are removed. Then the inequalities that are
2093 * redundant in the context of the equalities and inequalities of
2094 * context are removed.
2096 * First of all, we drop those constraints from "context"
2097 * that are irrelevant for computing the gist of "bset".
2098 * Alternatively, we could factorize the intersection of "context" and "bset".
2100 * We first compute the integer affine hull of the intersection,
2101 * compute the gist inside this affine hull and then add back
2102 * those equalities that are not implied by the context.
2104 * If two constraints are mutually redundant, then uset_gist_full
2105 * will remove the second of those constraints. We therefore first
2106 * sort the constraints so that constraints not involving existentially
2107 * quantified variables are given precedence over those that do.
2108 * We have to perform this sorting before the variable compression,
2109 * because that may effect the order of the variables.
2111 static __isl_give isl_basic_set *uset_gist(__isl_take isl_basic_set *bset,
2112 __isl_take isl_basic_set *context)
2114 isl_mat *eq;
2115 isl_mat *T, *T2;
2116 isl_basic_set *aff;
2117 isl_basic_set *aff_context;
2118 unsigned total;
2120 if (!bset || !context)
2121 goto error;
2123 context = drop_irrelevant_constraints(context, bset);
2125 aff = isl_basic_set_copy(bset);
2126 aff = isl_basic_set_intersect(aff, isl_basic_set_copy(context));
2127 aff = isl_basic_set_affine_hull(aff);
2128 if (!aff)
2129 goto error;
2130 if (isl_basic_set_plain_is_empty(aff)) {
2131 isl_basic_set_free(bset);
2132 isl_basic_set_free(context);
2133 return aff;
2135 bset = isl_basic_set_sort_constraints(bset);
2136 if (aff->n_eq == 0) {
2137 isl_basic_set_free(aff);
2138 return uset_gist_full(bset, context);
2140 total = isl_basic_set_total_dim(bset);
2141 eq = isl_mat_sub_alloc6(bset->ctx, aff->eq, 0, aff->n_eq, 0, 1 + total);
2142 eq = isl_mat_cow(eq);
2143 T = isl_mat_variable_compression(eq, &T2);
2144 if (T && T->n_col == 0) {
2145 isl_mat_free(T);
2146 isl_mat_free(T2);
2147 isl_basic_set_free(context);
2148 isl_basic_set_free(aff);
2149 return isl_basic_set_set_to_empty(bset);
2152 aff_context = isl_basic_set_affine_hull(isl_basic_set_copy(context));
2154 bset = isl_basic_set_preimage(bset, isl_mat_copy(T));
2155 context = isl_basic_set_preimage(context, T);
2157 bset = uset_gist_full(bset, context);
2158 bset = isl_basic_set_preimage(bset, T2);
2159 bset = isl_basic_set_intersect(bset, aff);
2160 bset = isl_basic_set_reduce_using_equalities(bset, aff_context);
2162 if (bset) {
2163 ISL_F_SET(bset, ISL_BASIC_SET_NO_IMPLICIT);
2164 ISL_F_SET(bset, ISL_BASIC_SET_NO_REDUNDANT);
2167 return bset;
2168 error:
2169 isl_basic_set_free(bset);
2170 isl_basic_set_free(context);
2171 return NULL;
2174 /* Normalize the divs in "bmap" in the context of the equalities in "context".
2175 * We simply add the equalities in context to bmap and then do a regular
2176 * div normalizations. Better results can be obtained by normalizing
2177 * only the divs in bmap than do not also appear in context.
2178 * We need to be careful to reduce the divs using the equalities
2179 * so that later calls to isl_basic_map_overlying_set wouldn't introduce
2180 * spurious constraints.
2182 static struct isl_basic_map *normalize_divs_in_context(
2183 struct isl_basic_map *bmap, struct isl_basic_map *context)
2185 int i;
2186 unsigned total_context;
2187 int div_eq;
2189 div_eq = n_pure_div_eq(bmap);
2190 if (div_eq == 0)
2191 return bmap;
2193 if (context->n_div > 0)
2194 bmap = isl_basic_map_align_divs(bmap, context);
2196 total_context = isl_basic_map_total_dim(context);
2197 bmap = isl_basic_map_extend_constraints(bmap, context->n_eq, 0);
2198 for (i = 0; i < context->n_eq; ++i) {
2199 int k;
2200 k = isl_basic_map_alloc_equality(bmap);
2201 if (k < 0)
2202 return isl_basic_map_free(bmap);
2203 isl_seq_cpy(bmap->eq[k], context->eq[i], 1 + total_context);
2204 isl_seq_clr(bmap->eq[k] + 1 + total_context,
2205 isl_basic_map_total_dim(bmap) - total_context);
2207 bmap = isl_basic_map_gauss(bmap, NULL);
2208 bmap = normalize_divs(bmap, NULL);
2209 bmap = isl_basic_map_gauss(bmap, NULL);
2210 return bmap;
2213 /* Return a basic map that has the same intersection with "context" as "bmap"
2214 * and that is as "simple" as possible.
2216 * The core computation is performed on the pure constraints.
2217 * When we add back the meaning of the integer divisions, we need
2218 * to (re)introduce the div constraints. If we happen to have
2219 * discovered that some of these integer divisions are equal to
2220 * some affine combination of other variables, then these div
2221 * constraints may end up getting simplified in terms of the equalities,
2222 * resulting in extra inequalities on the other variables that
2223 * may have been removed already or that may not even have been
2224 * part of the input. We try and remove those constraints of
2225 * this form that are most obviously redundant with respect to
2226 * the context. We also remove those div constraints that are
2227 * redundant with respect to the other constraints in the result.
2229 struct isl_basic_map *isl_basic_map_gist(struct isl_basic_map *bmap,
2230 struct isl_basic_map *context)
2232 isl_basic_set *bset, *eq;
2233 isl_basic_map *eq_bmap;
2234 unsigned n_div, n_eq, n_ineq;
2236 if (!bmap || !context)
2237 goto error;
2239 if (isl_basic_map_is_universe(bmap)) {
2240 isl_basic_map_free(context);
2241 return bmap;
2243 if (isl_basic_map_plain_is_empty(context)) {
2244 isl_space *space = isl_basic_map_get_space(bmap);
2245 isl_basic_map_free(bmap);
2246 isl_basic_map_free(context);
2247 return isl_basic_map_universe(space);
2249 if (isl_basic_map_plain_is_empty(bmap)) {
2250 isl_basic_map_free(context);
2251 return bmap;
2254 bmap = isl_basic_map_remove_redundancies(bmap);
2255 context = isl_basic_map_remove_redundancies(context);
2256 if (!context)
2257 goto error;
2259 if (context->n_eq)
2260 bmap = normalize_divs_in_context(bmap, context);
2262 context = isl_basic_map_align_divs(context, bmap);
2263 bmap = isl_basic_map_align_divs(bmap, context);
2264 n_div = isl_basic_map_dim(bmap, isl_dim_div);
2266 bset = uset_gist(isl_basic_map_underlying_set(isl_basic_map_copy(bmap)),
2267 isl_basic_map_underlying_set(isl_basic_map_copy(context)));
2269 if (!bset || bset->n_eq == 0 || n_div == 0 ||
2270 isl_basic_set_plain_is_empty(bset)) {
2271 isl_basic_map_free(context);
2272 return isl_basic_map_overlying_set(bset, bmap);
2275 n_eq = bset->n_eq;
2276 n_ineq = bset->n_ineq;
2277 eq = isl_basic_set_copy(bset);
2278 eq = isl_basic_set_cow(bset);
2279 if (isl_basic_set_free_inequality(eq, n_ineq) < 0)
2280 eq = isl_basic_set_free(eq);
2281 if (isl_basic_set_free_equality(bset, n_eq) < 0)
2282 bset = isl_basic_set_free(bset);
2284 eq_bmap = isl_basic_map_overlying_set(eq, isl_basic_map_copy(bmap));
2285 eq_bmap = isl_basic_map_remove_shifted_constraints(eq_bmap, context);
2286 bmap = isl_basic_map_overlying_set(bset, bmap);
2287 bmap = isl_basic_map_intersect(bmap, eq_bmap);
2288 bmap = isl_basic_map_remove_redundancies(bmap);
2290 return bmap;
2291 error:
2292 isl_basic_map_free(bmap);
2293 isl_basic_map_free(context);
2294 return NULL;
2298 * Assumes context has no implicit divs.
2300 __isl_give isl_map *isl_map_gist_basic_map(__isl_take isl_map *map,
2301 __isl_take isl_basic_map *context)
2303 int i;
2305 if (!map || !context)
2306 goto error;
2308 if (isl_basic_map_plain_is_empty(context)) {
2309 isl_space *space = isl_map_get_space(map);
2310 isl_map_free(map);
2311 isl_basic_map_free(context);
2312 return isl_map_universe(space);
2315 context = isl_basic_map_remove_redundancies(context);
2316 map = isl_map_cow(map);
2317 if (!map || !context)
2318 goto error;
2319 isl_assert(map->ctx, isl_space_is_equal(map->dim, context->dim), goto error);
2320 map = isl_map_compute_divs(map);
2321 if (!map)
2322 goto error;
2323 for (i = map->n - 1; i >= 0; --i) {
2324 map->p[i] = isl_basic_map_gist(map->p[i],
2325 isl_basic_map_copy(context));
2326 if (!map->p[i])
2327 goto error;
2328 if (isl_basic_map_plain_is_empty(map->p[i])) {
2329 isl_basic_map_free(map->p[i]);
2330 if (i != map->n - 1)
2331 map->p[i] = map->p[map->n - 1];
2332 map->n--;
2335 isl_basic_map_free(context);
2336 ISL_F_CLR(map, ISL_MAP_NORMALIZED);
2337 return map;
2338 error:
2339 isl_map_free(map);
2340 isl_basic_map_free(context);
2341 return NULL;
2344 /* Return a map that has the same intersection with "context" as "map"
2345 * and that is as "simple" as possible.
2347 * If "map" is already the universe, then we cannot make it any simpler.
2348 * Similarly, if "context" is the universe, then we cannot exploit it
2349 * to simplify "map"
2350 * If "map" and "context" are identical to each other, then we can
2351 * return the corresponding universe.
2353 * If none of these cases apply, we have to work a bit harder.
2355 static __isl_give isl_map *map_gist(__isl_take isl_map *map,
2356 __isl_take isl_map *context)
2358 int equal;
2359 int is_universe;
2361 is_universe = isl_map_plain_is_universe(map);
2362 if (is_universe >= 0 && !is_universe)
2363 is_universe = isl_map_plain_is_universe(context);
2364 if (is_universe < 0)
2365 goto error;
2366 if (is_universe) {
2367 isl_map_free(context);
2368 return map;
2371 equal = isl_map_plain_is_equal(map, context);
2372 if (equal < 0)
2373 goto error;
2374 if (equal) {
2375 isl_map *res = isl_map_universe(isl_map_get_space(map));
2376 isl_map_free(map);
2377 isl_map_free(context);
2378 return res;
2381 context = isl_map_compute_divs(context);
2382 return isl_map_gist_basic_map(map, isl_map_simple_hull(context));
2383 error:
2384 isl_map_free(map);
2385 isl_map_free(context);
2386 return NULL;
2389 __isl_give isl_map *isl_map_gist(__isl_take isl_map *map,
2390 __isl_take isl_map *context)
2392 return isl_map_align_params_map_map_and(map, context, &map_gist);
2395 struct isl_basic_set *isl_basic_set_gist(struct isl_basic_set *bset,
2396 struct isl_basic_set *context)
2398 return (struct isl_basic_set *)isl_basic_map_gist(
2399 (struct isl_basic_map *)bset, (struct isl_basic_map *)context);
2402 __isl_give isl_set *isl_set_gist_basic_set(__isl_take isl_set *set,
2403 __isl_take isl_basic_set *context)
2405 return (struct isl_set *)isl_map_gist_basic_map((struct isl_map *)set,
2406 (struct isl_basic_map *)context);
2409 __isl_give isl_set *isl_set_gist_params_basic_set(__isl_take isl_set *set,
2410 __isl_take isl_basic_set *context)
2412 isl_space *space = isl_set_get_space(set);
2413 isl_basic_set *dom_context = isl_basic_set_universe(space);
2414 dom_context = isl_basic_set_intersect_params(dom_context, context);
2415 return isl_set_gist_basic_set(set, dom_context);
2418 __isl_give isl_set *isl_set_gist(__isl_take isl_set *set,
2419 __isl_take isl_set *context)
2421 return (struct isl_set *)isl_map_gist((struct isl_map *)set,
2422 (struct isl_map *)context);
2425 __isl_give isl_map *isl_map_gist_domain(__isl_take isl_map *map,
2426 __isl_take isl_set *context)
2428 isl_map *map_context = isl_map_universe(isl_map_get_space(map));
2429 map_context = isl_map_intersect_domain(map_context, context);
2430 return isl_map_gist(map, map_context);
2433 __isl_give isl_map *isl_map_gist_range(__isl_take isl_map *map,
2434 __isl_take isl_set *context)
2436 isl_map *map_context = isl_map_universe(isl_map_get_space(map));
2437 map_context = isl_map_intersect_range(map_context, context);
2438 return isl_map_gist(map, map_context);
2441 __isl_give isl_map *isl_map_gist_params(__isl_take isl_map *map,
2442 __isl_take isl_set *context)
2444 isl_map *map_context = isl_map_universe(isl_map_get_space(map));
2445 map_context = isl_map_intersect_params(map_context, context);
2446 return isl_map_gist(map, map_context);
2449 __isl_give isl_set *isl_set_gist_params(__isl_take isl_set *set,
2450 __isl_take isl_set *context)
2452 return isl_map_gist_params(set, context);
2455 /* Quick check to see if two basic maps are disjoint.
2456 * In particular, we reduce the equalities and inequalities of
2457 * one basic map in the context of the equalities of the other
2458 * basic map and check if we get a contradiction.
2460 int isl_basic_map_plain_is_disjoint(__isl_keep isl_basic_map *bmap1,
2461 __isl_keep isl_basic_map *bmap2)
2463 struct isl_vec *v = NULL;
2464 int *elim = NULL;
2465 unsigned total;
2466 int i;
2468 if (!bmap1 || !bmap2)
2469 return -1;
2470 isl_assert(bmap1->ctx, isl_space_is_equal(bmap1->dim, bmap2->dim),
2471 return -1);
2472 if (bmap1->n_div || bmap2->n_div)
2473 return 0;
2474 if (!bmap1->n_eq && !bmap2->n_eq)
2475 return 0;
2477 total = isl_space_dim(bmap1->dim, isl_dim_all);
2478 if (total == 0)
2479 return 0;
2480 v = isl_vec_alloc(bmap1->ctx, 1 + total);
2481 if (!v)
2482 goto error;
2483 elim = isl_alloc_array(bmap1->ctx, int, total);
2484 if (!elim)
2485 goto error;
2486 compute_elimination_index(bmap1, elim);
2487 for (i = 0; i < bmap2->n_eq; ++i) {
2488 int reduced;
2489 reduced = reduced_using_equalities(v->block.data, bmap2->eq[i],
2490 bmap1, elim);
2491 if (reduced && !isl_int_is_zero(v->block.data[0]) &&
2492 isl_seq_first_non_zero(v->block.data + 1, total) == -1)
2493 goto disjoint;
2495 for (i = 0; i < bmap2->n_ineq; ++i) {
2496 int reduced;
2497 reduced = reduced_using_equalities(v->block.data,
2498 bmap2->ineq[i], bmap1, elim);
2499 if (reduced && isl_int_is_neg(v->block.data[0]) &&
2500 isl_seq_first_non_zero(v->block.data + 1, total) == -1)
2501 goto disjoint;
2503 compute_elimination_index(bmap2, elim);
2504 for (i = 0; i < bmap1->n_ineq; ++i) {
2505 int reduced;
2506 reduced = reduced_using_equalities(v->block.data,
2507 bmap1->ineq[i], bmap2, elim);
2508 if (reduced && isl_int_is_neg(v->block.data[0]) &&
2509 isl_seq_first_non_zero(v->block.data + 1, total) == -1)
2510 goto disjoint;
2512 isl_vec_free(v);
2513 free(elim);
2514 return 0;
2515 disjoint:
2516 isl_vec_free(v);
2517 free(elim);
2518 return 1;
2519 error:
2520 isl_vec_free(v);
2521 free(elim);
2522 return -1;
2525 int isl_basic_set_plain_is_disjoint(__isl_keep isl_basic_set *bset1,
2526 __isl_keep isl_basic_set *bset2)
2528 return isl_basic_map_plain_is_disjoint((struct isl_basic_map *)bset1,
2529 (struct isl_basic_map *)bset2);
2532 /* Are "map1" and "map2" obviously disjoint?
2534 * If one of them is empty or if they live in different spaces (ignoring
2535 * parameters), then they are clearly disjoint.
2537 * If they have different parameters, then we skip any further tests.
2539 * If they are obviously equal, but not obviously empty, then we will
2540 * not be able to detect if they are disjoint.
2542 * Otherwise we check if each basic map in "map1" is obviously disjoint
2543 * from each basic map in "map2".
2545 int isl_map_plain_is_disjoint(__isl_keep isl_map *map1,
2546 __isl_keep isl_map *map2)
2548 int i, j;
2549 int disjoint;
2550 int intersect;
2551 int match;
2553 if (!map1 || !map2)
2554 return -1;
2556 disjoint = isl_map_plain_is_empty(map1);
2557 if (disjoint < 0 || disjoint)
2558 return disjoint;
2560 disjoint = isl_map_plain_is_empty(map2);
2561 if (disjoint < 0 || disjoint)
2562 return disjoint;
2564 match = isl_space_tuple_is_equal(map1->dim, isl_dim_in,
2565 map2->dim, isl_dim_in);
2566 if (match < 0 || !match)
2567 return match < 0 ? -1 : 1;
2569 match = isl_space_tuple_is_equal(map1->dim, isl_dim_out,
2570 map2->dim, isl_dim_out);
2571 if (match < 0 || !match)
2572 return match < 0 ? -1 : 1;
2574 match = isl_space_match(map1->dim, isl_dim_param,
2575 map2->dim, isl_dim_param);
2576 if (match < 0 || !match)
2577 return match < 0 ? -1 : 0;
2579 intersect = isl_map_plain_is_equal(map1, map2);
2580 if (intersect < 0 || intersect)
2581 return intersect < 0 ? -1 : 0;
2583 for (i = 0; i < map1->n; ++i) {
2584 for (j = 0; j < map2->n; ++j) {
2585 int d = isl_basic_map_plain_is_disjoint(map1->p[i],
2586 map2->p[j]);
2587 if (d != 1)
2588 return d;
2591 return 1;
2594 /* Are "map1" and "map2" disjoint?
2596 * They are disjoint if they are "obviously disjoint" or if one of them
2597 * is empty. Otherwise, they are not disjoint if one of them is universal.
2598 * If none of these cases apply, we compute the intersection and see if
2599 * the result is empty.
2601 int isl_map_is_disjoint(__isl_keep isl_map *map1, __isl_keep isl_map *map2)
2603 int disjoint;
2604 int intersect;
2605 isl_map *test;
2607 disjoint = isl_map_plain_is_disjoint(map1, map2);
2608 if (disjoint < 0 || disjoint)
2609 return disjoint;
2611 disjoint = isl_map_is_empty(map1);
2612 if (disjoint < 0 || disjoint)
2613 return disjoint;
2615 disjoint = isl_map_is_empty(map2);
2616 if (disjoint < 0 || disjoint)
2617 return disjoint;
2619 intersect = isl_map_plain_is_universe(map1);
2620 if (intersect < 0 || intersect)
2621 return intersect < 0 ? -1 : 0;
2623 intersect = isl_map_plain_is_universe(map2);
2624 if (intersect < 0 || intersect)
2625 return intersect < 0 ? -1 : 0;
2627 test = isl_map_intersect(isl_map_copy(map1), isl_map_copy(map2));
2628 disjoint = isl_map_is_empty(test);
2629 isl_map_free(test);
2631 return disjoint;
2634 /* Are "bmap1" and "bmap2" disjoint?
2636 * They are disjoint if they are "obviously disjoint" or if one of them
2637 * is empty. Otherwise, they are not disjoint if one of them is universal.
2638 * If none of these cases apply, we compute the intersection and see if
2639 * the result is empty.
2641 int isl_basic_map_is_disjoint(__isl_keep isl_basic_map *bmap1,
2642 __isl_keep isl_basic_map *bmap2)
2644 int disjoint;
2645 int intersect;
2646 isl_basic_map *test;
2648 disjoint = isl_basic_map_plain_is_disjoint(bmap1, bmap2);
2649 if (disjoint < 0 || disjoint)
2650 return disjoint;
2652 disjoint = isl_basic_map_is_empty(bmap1);
2653 if (disjoint < 0 || disjoint)
2654 return disjoint;
2656 disjoint = isl_basic_map_is_empty(bmap2);
2657 if (disjoint < 0 || disjoint)
2658 return disjoint;
2660 intersect = isl_basic_map_is_universe(bmap1);
2661 if (intersect < 0 || intersect)
2662 return intersect < 0 ? -1 : 0;
2664 intersect = isl_basic_map_is_universe(bmap2);
2665 if (intersect < 0 || intersect)
2666 return intersect < 0 ? -1 : 0;
2668 test = isl_basic_map_intersect(isl_basic_map_copy(bmap1),
2669 isl_basic_map_copy(bmap2));
2670 disjoint = isl_basic_map_is_empty(test);
2671 isl_basic_map_free(test);
2673 return disjoint;
2676 /* Are "bset1" and "bset2" disjoint?
2678 int isl_basic_set_is_disjoint(__isl_keep isl_basic_set *bset1,
2679 __isl_keep isl_basic_set *bset2)
2681 return isl_basic_map_is_disjoint(bset1, bset2);
2684 int isl_set_plain_is_disjoint(__isl_keep isl_set *set1,
2685 __isl_keep isl_set *set2)
2687 return isl_map_plain_is_disjoint((struct isl_map *)set1,
2688 (struct isl_map *)set2);
2691 /* Are "set1" and "set2" disjoint?
2693 int isl_set_is_disjoint(__isl_keep isl_set *set1, __isl_keep isl_set *set2)
2695 return isl_map_is_disjoint(set1, set2);
2698 int isl_set_fast_is_disjoint(__isl_keep isl_set *set1, __isl_keep isl_set *set2)
2700 return isl_set_plain_is_disjoint(set1, set2);
2703 /* Check if we can combine a given div with lower bound l and upper
2704 * bound u with some other div and if so return that other div.
2705 * Otherwise return -1.
2707 * We first check that
2708 * - the bounds are opposites of each other (except for the constant
2709 * term)
2710 * - the bounds do not reference any other div
2711 * - no div is defined in terms of this div
2713 * Let m be the size of the range allowed on the div by the bounds.
2714 * That is, the bounds are of the form
2716 * e <= a <= e + m - 1
2718 * with e some expression in the other variables.
2719 * We look for another div b such that no third div is defined in terms
2720 * of this second div b and such that in any constraint that contains
2721 * a (except for the given lower and upper bound), also contains b
2722 * with a coefficient that is m times that of b.
2723 * That is, all constraints (execpt for the lower and upper bound)
2724 * are of the form
2726 * e + f (a + m b) >= 0
2728 * If so, we return b so that "a + m b" can be replaced by
2729 * a single div "c = a + m b".
2731 static int div_find_coalesce(struct isl_basic_map *bmap, int *pairs,
2732 unsigned div, unsigned l, unsigned u)
2734 int i, j;
2735 unsigned dim;
2736 int coalesce = -1;
2738 if (bmap->n_div <= 1)
2739 return -1;
2740 dim = isl_space_dim(bmap->dim, isl_dim_all);
2741 if (isl_seq_first_non_zero(bmap->ineq[l] + 1 + dim, div) != -1)
2742 return -1;
2743 if (isl_seq_first_non_zero(bmap->ineq[l] + 1 + dim + div + 1,
2744 bmap->n_div - div - 1) != -1)
2745 return -1;
2746 if (!isl_seq_is_neg(bmap->ineq[l] + 1, bmap->ineq[u] + 1,
2747 dim + bmap->n_div))
2748 return -1;
2750 for (i = 0; i < bmap->n_div; ++i) {
2751 if (isl_int_is_zero(bmap->div[i][0]))
2752 continue;
2753 if (!isl_int_is_zero(bmap->div[i][1 + 1 + dim + div]))
2754 return -1;
2757 isl_int_add(bmap->ineq[l][0], bmap->ineq[l][0], bmap->ineq[u][0]);
2758 if (isl_int_is_neg(bmap->ineq[l][0])) {
2759 isl_int_sub(bmap->ineq[l][0],
2760 bmap->ineq[l][0], bmap->ineq[u][0]);
2761 bmap = isl_basic_map_copy(bmap);
2762 bmap = isl_basic_map_set_to_empty(bmap);
2763 isl_basic_map_free(bmap);
2764 return -1;
2766 isl_int_add_ui(bmap->ineq[l][0], bmap->ineq[l][0], 1);
2767 for (i = 0; i < bmap->n_div; ++i) {
2768 if (i == div)
2769 continue;
2770 if (!pairs[i])
2771 continue;
2772 for (j = 0; j < bmap->n_div; ++j) {
2773 if (isl_int_is_zero(bmap->div[j][0]))
2774 continue;
2775 if (!isl_int_is_zero(bmap->div[j][1 + 1 + dim + i]))
2776 break;
2778 if (j < bmap->n_div)
2779 continue;
2780 for (j = 0; j < bmap->n_ineq; ++j) {
2781 int valid;
2782 if (j == l || j == u)
2783 continue;
2784 if (isl_int_is_zero(bmap->ineq[j][1 + dim + div]))
2785 continue;
2786 if (isl_int_is_zero(bmap->ineq[j][1 + dim + i]))
2787 break;
2788 isl_int_mul(bmap->ineq[j][1 + dim + div],
2789 bmap->ineq[j][1 + dim + div],
2790 bmap->ineq[l][0]);
2791 valid = isl_int_eq(bmap->ineq[j][1 + dim + div],
2792 bmap->ineq[j][1 + dim + i]);
2793 isl_int_divexact(bmap->ineq[j][1 + dim + div],
2794 bmap->ineq[j][1 + dim + div],
2795 bmap->ineq[l][0]);
2796 if (!valid)
2797 break;
2799 if (j < bmap->n_ineq)
2800 continue;
2801 coalesce = i;
2802 break;
2804 isl_int_sub_ui(bmap->ineq[l][0], bmap->ineq[l][0], 1);
2805 isl_int_sub(bmap->ineq[l][0], bmap->ineq[l][0], bmap->ineq[u][0]);
2806 return coalesce;
2809 /* Given a lower and an upper bound on div i, construct an inequality
2810 * that when nonnegative ensures that this pair of bounds always allows
2811 * for an integer value of the given div.
2812 * The lower bound is inequality l, while the upper bound is inequality u.
2813 * The constructed inequality is stored in ineq.
2814 * g, fl, fu are temporary scalars.
2816 * Let the upper bound be
2818 * -n_u a + e_u >= 0
2820 * and the lower bound
2822 * n_l a + e_l >= 0
2824 * Let n_u = f_u g and n_l = f_l g, with g = gcd(n_u, n_l).
2825 * We have
2827 * - f_u e_l <= f_u f_l g a <= f_l e_u
2829 * Since all variables are integer valued, this is equivalent to
2831 * - f_u e_l - (f_u - 1) <= f_u f_l g a <= f_l e_u + (f_l - 1)
2833 * If this interval is at least f_u f_l g, then it contains at least
2834 * one integer value for a.
2835 * That is, the test constraint is
2837 * f_l e_u + f_u e_l + f_l - 1 + f_u - 1 + 1 >= f_u f_l g
2839 static void construct_test_ineq(struct isl_basic_map *bmap, int i,
2840 int l, int u, isl_int *ineq, isl_int g, isl_int fl, isl_int fu)
2842 unsigned dim;
2843 dim = isl_space_dim(bmap->dim, isl_dim_all);
2845 isl_int_gcd(g, bmap->ineq[l][1 + dim + i], bmap->ineq[u][1 + dim + i]);
2846 isl_int_divexact(fl, bmap->ineq[l][1 + dim + i], g);
2847 isl_int_divexact(fu, bmap->ineq[u][1 + dim + i], g);
2848 isl_int_neg(fu, fu);
2849 isl_seq_combine(ineq, fl, bmap->ineq[u], fu, bmap->ineq[l],
2850 1 + dim + bmap->n_div);
2851 isl_int_add(ineq[0], ineq[0], fl);
2852 isl_int_add(ineq[0], ineq[0], fu);
2853 isl_int_sub_ui(ineq[0], ineq[0], 1);
2854 isl_int_mul(g, g, fl);
2855 isl_int_mul(g, g, fu);
2856 isl_int_sub(ineq[0], ineq[0], g);
2859 /* Remove more kinds of divs that are not strictly needed.
2860 * In particular, if all pairs of lower and upper bounds on a div
2861 * are such that they allow at least one integer value of the div,
2862 * the we can eliminate the div using Fourier-Motzkin without
2863 * introducing any spurious solutions.
2865 static struct isl_basic_map *drop_more_redundant_divs(
2866 struct isl_basic_map *bmap, int *pairs, int n)
2868 struct isl_tab *tab = NULL;
2869 struct isl_vec *vec = NULL;
2870 unsigned dim;
2871 int remove = -1;
2872 isl_int g, fl, fu;
2874 isl_int_init(g);
2875 isl_int_init(fl);
2876 isl_int_init(fu);
2878 if (!bmap)
2879 goto error;
2881 dim = isl_space_dim(bmap->dim, isl_dim_all);
2882 vec = isl_vec_alloc(bmap->ctx, 1 + dim + bmap->n_div);
2883 if (!vec)
2884 goto error;
2886 tab = isl_tab_from_basic_map(bmap, 0);
2888 while (n > 0) {
2889 int i, l, u;
2890 int best = -1;
2891 enum isl_lp_result res;
2893 for (i = 0; i < bmap->n_div; ++i) {
2894 if (!pairs[i])
2895 continue;
2896 if (best >= 0 && pairs[best] <= pairs[i])
2897 continue;
2898 best = i;
2901 i = best;
2902 for (l = 0; l < bmap->n_ineq; ++l) {
2903 if (!isl_int_is_pos(bmap->ineq[l][1 + dim + i]))
2904 continue;
2905 for (u = 0; u < bmap->n_ineq; ++u) {
2906 if (!isl_int_is_neg(bmap->ineq[u][1 + dim + i]))
2907 continue;
2908 construct_test_ineq(bmap, i, l, u,
2909 vec->el, g, fl, fu);
2910 res = isl_tab_min(tab, vec->el,
2911 bmap->ctx->one, &g, NULL, 0);
2912 if (res == isl_lp_error)
2913 goto error;
2914 if (res == isl_lp_empty) {
2915 bmap = isl_basic_map_set_to_empty(bmap);
2916 break;
2918 if (res != isl_lp_ok || isl_int_is_neg(g))
2919 break;
2921 if (u < bmap->n_ineq)
2922 break;
2924 if (l == bmap->n_ineq) {
2925 remove = i;
2926 break;
2928 pairs[i] = 0;
2929 --n;
2932 isl_tab_free(tab);
2933 isl_vec_free(vec);
2935 isl_int_clear(g);
2936 isl_int_clear(fl);
2937 isl_int_clear(fu);
2939 free(pairs);
2941 if (remove < 0)
2942 return bmap;
2944 bmap = isl_basic_map_remove_dims(bmap, isl_dim_div, remove, 1);
2945 return isl_basic_map_drop_redundant_divs(bmap);
2946 error:
2947 free(pairs);
2948 isl_basic_map_free(bmap);
2949 isl_tab_free(tab);
2950 isl_vec_free(vec);
2951 isl_int_clear(g);
2952 isl_int_clear(fl);
2953 isl_int_clear(fu);
2954 return NULL;
2957 /* Given a pair of divs div1 and div2 such that, expect for the lower bound l
2958 * and the upper bound u, div1 always occurs together with div2 in the form
2959 * (div1 + m div2), where m is the constant range on the variable div1
2960 * allowed by l and u, replace the pair div1 and div2 by a single
2961 * div that is equal to div1 + m div2.
2963 * The new div will appear in the location that contains div2.
2964 * We need to modify all constraints that contain
2965 * div2 = (div - div1) / m
2966 * (If a constraint does not contain div2, it will also not contain div1.)
2967 * If the constraint also contains div1, then we know they appear
2968 * as f (div1 + m div2) and we can simply replace (div1 + m div2) by div,
2969 * i.e., the coefficient of div is f.
2971 * Otherwise, we first need to introduce div1 into the constraint.
2972 * Let the l be
2974 * div1 + f >=0
2976 * and u
2978 * -div1 + f' >= 0
2980 * A lower bound on div2
2982 * n div2 + t >= 0
2984 * can be replaced by
2986 * (n * (m div 2 + div1) + m t + n f)/g >= 0
2988 * with g = gcd(m,n).
2989 * An upper bound
2991 * -n div2 + t >= 0
2993 * can be replaced by
2995 * (-n * (m div2 + div1) + m t + n f')/g >= 0
2997 * These constraint are those that we would obtain from eliminating
2998 * div1 using Fourier-Motzkin.
3000 * After all constraints have been modified, we drop the lower and upper
3001 * bound and then drop div1.
3003 static struct isl_basic_map *coalesce_divs(struct isl_basic_map *bmap,
3004 unsigned div1, unsigned div2, unsigned l, unsigned u)
3006 isl_int a;
3007 isl_int b;
3008 isl_int m;
3009 unsigned dim, total;
3010 int i;
3012 dim = isl_space_dim(bmap->dim, isl_dim_all);
3013 total = 1 + dim + bmap->n_div;
3015 isl_int_init(a);
3016 isl_int_init(b);
3017 isl_int_init(m);
3018 isl_int_add(m, bmap->ineq[l][0], bmap->ineq[u][0]);
3019 isl_int_add_ui(m, m, 1);
3021 for (i = 0; i < bmap->n_ineq; ++i) {
3022 if (i == l || i == u)
3023 continue;
3024 if (isl_int_is_zero(bmap->ineq[i][1 + dim + div2]))
3025 continue;
3026 if (isl_int_is_zero(bmap->ineq[i][1 + dim + div1])) {
3027 isl_int_gcd(b, m, bmap->ineq[i][1 + dim + div2]);
3028 isl_int_divexact(a, m, b);
3029 isl_int_divexact(b, bmap->ineq[i][1 + dim + div2], b);
3030 if (isl_int_is_pos(b)) {
3031 isl_seq_combine(bmap->ineq[i], a, bmap->ineq[i],
3032 b, bmap->ineq[l], total);
3033 } else {
3034 isl_int_neg(b, b);
3035 isl_seq_combine(bmap->ineq[i], a, bmap->ineq[i],
3036 b, bmap->ineq[u], total);
3039 isl_int_set(bmap->ineq[i][1 + dim + div2],
3040 bmap->ineq[i][1 + dim + div1]);
3041 isl_int_set_si(bmap->ineq[i][1 + dim + div1], 0);
3044 isl_int_clear(a);
3045 isl_int_clear(b);
3046 isl_int_clear(m);
3047 if (l > u) {
3048 isl_basic_map_drop_inequality(bmap, l);
3049 isl_basic_map_drop_inequality(bmap, u);
3050 } else {
3051 isl_basic_map_drop_inequality(bmap, u);
3052 isl_basic_map_drop_inequality(bmap, l);
3054 bmap = isl_basic_map_drop_div(bmap, div1);
3055 return bmap;
3058 /* First check if we can coalesce any pair of divs and
3059 * then continue with dropping more redundant divs.
3061 * We loop over all pairs of lower and upper bounds on a div
3062 * with coefficient 1 and -1, respectively, check if there
3063 * is any other div "c" with which we can coalesce the div
3064 * and if so, perform the coalescing.
3066 static struct isl_basic_map *coalesce_or_drop_more_redundant_divs(
3067 struct isl_basic_map *bmap, int *pairs, int n)
3069 int i, l, u;
3070 unsigned dim;
3072 dim = isl_space_dim(bmap->dim, isl_dim_all);
3074 for (i = 0; i < bmap->n_div; ++i) {
3075 if (!pairs[i])
3076 continue;
3077 for (l = 0; l < bmap->n_ineq; ++l) {
3078 if (!isl_int_is_one(bmap->ineq[l][1 + dim + i]))
3079 continue;
3080 for (u = 0; u < bmap->n_ineq; ++u) {
3081 int c;
3083 if (!isl_int_is_negone(bmap->ineq[u][1+dim+i]))
3084 continue;
3085 c = div_find_coalesce(bmap, pairs, i, l, u);
3086 if (c < 0)
3087 continue;
3088 free(pairs);
3089 bmap = coalesce_divs(bmap, i, c, l, u);
3090 return isl_basic_map_drop_redundant_divs(bmap);
3095 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
3096 return bmap;
3098 return drop_more_redundant_divs(bmap, pairs, n);
3101 /* Remove divs that are not strictly needed.
3102 * In particular, if a div only occurs positively (or negatively)
3103 * in constraints, then it can simply be dropped.
3104 * Also, if a div occurs in only two constraints and if moreover
3105 * those two constraints are opposite to each other, except for the constant
3106 * term and if the sum of the constant terms is such that for any value
3107 * of the other values, there is always at least one integer value of the
3108 * div, i.e., if one plus this sum is greater than or equal to
3109 * the (absolute value) of the coefficent of the div in the constraints,
3110 * then we can also simply drop the div.
3112 * We skip divs that appear in equalities or in the definition of other divs.
3113 * Divs that appear in the definition of other divs usually occur in at least
3114 * 4 constraints, but the constraints may have been simplified.
3116 * If any divs are left after these simple checks then we move on
3117 * to more complicated cases in drop_more_redundant_divs.
3119 struct isl_basic_map *isl_basic_map_drop_redundant_divs(
3120 struct isl_basic_map *bmap)
3122 int i, j;
3123 unsigned off;
3124 int *pairs = NULL;
3125 int n = 0;
3127 if (!bmap)
3128 goto error;
3129 if (bmap->n_div == 0)
3130 return bmap;
3132 off = isl_space_dim(bmap->dim, isl_dim_all);
3133 pairs = isl_calloc_array(bmap->ctx, int, bmap->n_div);
3134 if (!pairs)
3135 goto error;
3137 for (i = 0; i < bmap->n_div; ++i) {
3138 int pos, neg;
3139 int last_pos, last_neg;
3140 int redundant;
3141 int defined;
3143 defined = !isl_int_is_zero(bmap->div[i][0]);
3144 for (j = i; j < bmap->n_div; ++j)
3145 if (!isl_int_is_zero(bmap->div[j][1 + 1 + off + i]))
3146 break;
3147 if (j < bmap->n_div)
3148 continue;
3149 for (j = 0; j < bmap->n_eq; ++j)
3150 if (!isl_int_is_zero(bmap->eq[j][1 + off + i]))
3151 break;
3152 if (j < bmap->n_eq)
3153 continue;
3154 ++n;
3155 pos = neg = 0;
3156 for (j = 0; j < bmap->n_ineq; ++j) {
3157 if (isl_int_is_pos(bmap->ineq[j][1 + off + i])) {
3158 last_pos = j;
3159 ++pos;
3161 if (isl_int_is_neg(bmap->ineq[j][1 + off + i])) {
3162 last_neg = j;
3163 ++neg;
3166 pairs[i] = pos * neg;
3167 if (pairs[i] == 0) {
3168 for (j = bmap->n_ineq - 1; j >= 0; --j)
3169 if (!isl_int_is_zero(bmap->ineq[j][1+off+i]))
3170 isl_basic_map_drop_inequality(bmap, j);
3171 bmap = isl_basic_map_drop_div(bmap, i);
3172 free(pairs);
3173 return isl_basic_map_drop_redundant_divs(bmap);
3175 if (pairs[i] != 1)
3176 continue;
3177 if (!isl_seq_is_neg(bmap->ineq[last_pos] + 1,
3178 bmap->ineq[last_neg] + 1,
3179 off + bmap->n_div))
3180 continue;
3182 isl_int_add(bmap->ineq[last_pos][0],
3183 bmap->ineq[last_pos][0], bmap->ineq[last_neg][0]);
3184 isl_int_add_ui(bmap->ineq[last_pos][0],
3185 bmap->ineq[last_pos][0], 1);
3186 redundant = isl_int_ge(bmap->ineq[last_pos][0],
3187 bmap->ineq[last_pos][1+off+i]);
3188 isl_int_sub_ui(bmap->ineq[last_pos][0],
3189 bmap->ineq[last_pos][0], 1);
3190 isl_int_sub(bmap->ineq[last_pos][0],
3191 bmap->ineq[last_pos][0], bmap->ineq[last_neg][0]);
3192 if (!redundant) {
3193 if (defined ||
3194 !ok_to_set_div_from_bound(bmap, i, last_pos)) {
3195 pairs[i] = 0;
3196 --n;
3197 continue;
3199 bmap = set_div_from_lower_bound(bmap, i, last_pos);
3200 bmap = isl_basic_map_simplify(bmap);
3201 free(pairs);
3202 return isl_basic_map_drop_redundant_divs(bmap);
3204 if (last_pos > last_neg) {
3205 isl_basic_map_drop_inequality(bmap, last_pos);
3206 isl_basic_map_drop_inequality(bmap, last_neg);
3207 } else {
3208 isl_basic_map_drop_inequality(bmap, last_neg);
3209 isl_basic_map_drop_inequality(bmap, last_pos);
3211 bmap = isl_basic_map_drop_div(bmap, i);
3212 free(pairs);
3213 return isl_basic_map_drop_redundant_divs(bmap);
3216 if (n > 0)
3217 return coalesce_or_drop_more_redundant_divs(bmap, pairs, n);
3219 free(pairs);
3220 return bmap;
3221 error:
3222 free(pairs);
3223 isl_basic_map_free(bmap);
3224 return NULL;
3227 struct isl_basic_set *isl_basic_set_drop_redundant_divs(
3228 struct isl_basic_set *bset)
3230 return (struct isl_basic_set *)
3231 isl_basic_map_drop_redundant_divs((struct isl_basic_map *)bset);
3234 struct isl_map *isl_map_drop_redundant_divs(struct isl_map *map)
3236 int i;
3238 if (!map)
3239 return NULL;
3240 for (i = 0; i < map->n; ++i) {
3241 map->p[i] = isl_basic_map_drop_redundant_divs(map->p[i]);
3242 if (!map->p[i])
3243 goto error;
3245 ISL_F_CLR(map, ISL_MAP_NORMALIZED);
3246 return map;
3247 error:
3248 isl_map_free(map);
3249 return NULL;
3252 struct isl_set *isl_set_drop_redundant_divs(struct isl_set *set)
3254 return (struct isl_set *)
3255 isl_map_drop_redundant_divs((struct isl_map *)set);