isl_map_simplify.c: extract out uset_gist_{un,}compressed
[isl.git] / isl_map_simplify.c
blobe24378c9336263a9fb9a78a86cb817d0e2a684d0
1 /*
2 * Copyright 2008-2009 Katholieke Universiteit Leuven
3 * Copyright 2012-2013 Ecole Normale Superieure
4 * Copyright 2014-2015 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 <isl_ctx_private.h>
16 #include <isl_map_private.h>
17 #include "isl_equalities.h"
18 #include <isl/map.h>
19 #include <isl_seq.h>
20 #include "isl_tab.h"
21 #include <isl_space_private.h>
22 #include <isl_mat_private.h>
23 #include <isl_vec_private.h>
25 static void swap_equality(struct isl_basic_map *bmap, int a, int b)
27 isl_int *t = bmap->eq[a];
28 bmap->eq[a] = bmap->eq[b];
29 bmap->eq[b] = t;
32 static void swap_inequality(struct isl_basic_map *bmap, int a, int b)
34 if (a != b) {
35 isl_int *t = bmap->ineq[a];
36 bmap->ineq[a] = bmap->ineq[b];
37 bmap->ineq[b] = t;
41 static void constraint_drop_vars(isl_int *c, unsigned n, unsigned rem)
43 isl_seq_cpy(c, c + n, rem);
44 isl_seq_clr(c + rem, n);
47 /* Drop n dimensions starting at first.
49 * In principle, this frees up some extra variables as the number
50 * of columns remains constant, but we would have to extend
51 * the div array too as the number of rows in this array is assumed
52 * to be equal to extra.
54 struct isl_basic_set *isl_basic_set_drop_dims(
55 struct isl_basic_set *bset, unsigned first, unsigned n)
57 int i;
59 if (!bset)
60 goto error;
62 isl_assert(bset->ctx, first + n <= bset->dim->n_out, goto error);
64 if (n == 0 && !isl_space_get_tuple_name(bset->dim, isl_dim_set))
65 return bset;
67 bset = isl_basic_set_cow(bset);
68 if (!bset)
69 return NULL;
71 for (i = 0; i < bset->n_eq; ++i)
72 constraint_drop_vars(bset->eq[i]+1+bset->dim->nparam+first, n,
73 (bset->dim->n_out-first-n)+bset->extra);
75 for (i = 0; i < bset->n_ineq; ++i)
76 constraint_drop_vars(bset->ineq[i]+1+bset->dim->nparam+first, n,
77 (bset->dim->n_out-first-n)+bset->extra);
79 for (i = 0; i < bset->n_div; ++i)
80 constraint_drop_vars(bset->div[i]+1+1+bset->dim->nparam+first, n,
81 (bset->dim->n_out-first-n)+bset->extra);
83 bset->dim = isl_space_drop_outputs(bset->dim, first, n);
84 if (!bset->dim)
85 goto error;
87 ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED);
88 bset = isl_basic_set_simplify(bset);
89 return isl_basic_set_finalize(bset);
90 error:
91 isl_basic_set_free(bset);
92 return NULL;
95 struct isl_set *isl_set_drop_dims(
96 struct isl_set *set, unsigned first, unsigned n)
98 int i;
100 if (!set)
101 goto error;
103 isl_assert(set->ctx, first + n <= set->dim->n_out, goto error);
105 if (n == 0 && !isl_space_get_tuple_name(set->dim, isl_dim_set))
106 return set;
107 set = isl_set_cow(set);
108 if (!set)
109 goto error;
110 set->dim = isl_space_drop_outputs(set->dim, first, n);
111 if (!set->dim)
112 goto error;
114 for (i = 0; i < set->n; ++i) {
115 set->p[i] = isl_basic_set_drop_dims(set->p[i], first, n);
116 if (!set->p[i])
117 goto error;
120 ISL_F_CLR(set, ISL_SET_NORMALIZED);
121 return set;
122 error:
123 isl_set_free(set);
124 return NULL;
127 /* Move "n" divs starting at "first" to the end of the list of divs.
129 static struct isl_basic_map *move_divs_last(struct isl_basic_map *bmap,
130 unsigned first, unsigned n)
132 isl_int **div;
133 int i;
135 if (first + n == bmap->n_div)
136 return bmap;
138 div = isl_alloc_array(bmap->ctx, isl_int *, n);
139 if (!div)
140 goto error;
141 for (i = 0; i < n; ++i)
142 div[i] = bmap->div[first + i];
143 for (i = 0; i < bmap->n_div - first - n; ++i)
144 bmap->div[first + i] = bmap->div[first + n + i];
145 for (i = 0; i < n; ++i)
146 bmap->div[bmap->n_div - n + i] = div[i];
147 free(div);
148 return bmap;
149 error:
150 isl_basic_map_free(bmap);
151 return NULL;
154 /* Drop "n" dimensions of type "type" starting at "first".
156 * In principle, this frees up some extra variables as the number
157 * of columns remains constant, but we would have to extend
158 * the div array too as the number of rows in this array is assumed
159 * to be equal to extra.
161 struct isl_basic_map *isl_basic_map_drop(struct isl_basic_map *bmap,
162 enum isl_dim_type type, unsigned first, unsigned n)
164 int i;
165 unsigned dim;
166 unsigned offset;
167 unsigned left;
169 if (!bmap)
170 goto error;
172 dim = isl_basic_map_dim(bmap, type);
173 isl_assert(bmap->ctx, first + n <= dim, goto error);
175 if (n == 0 && !isl_space_is_named_or_nested(bmap->dim, type))
176 return bmap;
178 bmap = isl_basic_map_cow(bmap);
179 if (!bmap)
180 return NULL;
182 offset = isl_basic_map_offset(bmap, type) + first;
183 left = isl_basic_map_total_dim(bmap) - (offset - 1) - n;
184 for (i = 0; i < bmap->n_eq; ++i)
185 constraint_drop_vars(bmap->eq[i]+offset, n, left);
187 for (i = 0; i < bmap->n_ineq; ++i)
188 constraint_drop_vars(bmap->ineq[i]+offset, n, left);
190 for (i = 0; i < bmap->n_div; ++i)
191 constraint_drop_vars(bmap->div[i]+1+offset, n, left);
193 if (type == isl_dim_div) {
194 bmap = move_divs_last(bmap, first, n);
195 if (!bmap)
196 goto error;
197 isl_basic_map_free_div(bmap, n);
198 } else
199 bmap->dim = isl_space_drop_dims(bmap->dim, type, first, n);
200 if (!bmap->dim)
201 goto error;
203 ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
204 bmap = isl_basic_map_simplify(bmap);
205 return isl_basic_map_finalize(bmap);
206 error:
207 isl_basic_map_free(bmap);
208 return NULL;
211 __isl_give isl_basic_set *isl_basic_set_drop(__isl_take isl_basic_set *bset,
212 enum isl_dim_type type, unsigned first, unsigned n)
214 return (isl_basic_set *)isl_basic_map_drop((isl_basic_map *)bset,
215 type, first, n);
218 struct isl_basic_map *isl_basic_map_drop_inputs(
219 struct isl_basic_map *bmap, unsigned first, unsigned n)
221 return isl_basic_map_drop(bmap, isl_dim_in, first, n);
224 struct isl_map *isl_map_drop(struct isl_map *map,
225 enum isl_dim_type type, unsigned first, unsigned n)
227 int i;
229 if (!map)
230 goto error;
232 isl_assert(map->ctx, first + n <= isl_map_dim(map, type), goto error);
234 if (n == 0 && !isl_space_get_tuple_name(map->dim, type))
235 return map;
236 map = isl_map_cow(map);
237 if (!map)
238 goto error;
239 map->dim = isl_space_drop_dims(map->dim, type, first, n);
240 if (!map->dim)
241 goto error;
243 for (i = 0; i < map->n; ++i) {
244 map->p[i] = isl_basic_map_drop(map->p[i], type, first, n);
245 if (!map->p[i])
246 goto error;
248 ISL_F_CLR(map, ISL_MAP_NORMALIZED);
250 return map;
251 error:
252 isl_map_free(map);
253 return NULL;
256 struct isl_set *isl_set_drop(struct isl_set *set,
257 enum isl_dim_type type, unsigned first, unsigned n)
259 return (isl_set *)isl_map_drop((isl_map *)set, type, first, n);
262 struct isl_map *isl_map_drop_inputs(
263 struct isl_map *map, unsigned first, unsigned n)
265 return isl_map_drop(map, isl_dim_in, first, n);
269 * We don't cow, as the div is assumed to be redundant.
271 static struct isl_basic_map *isl_basic_map_drop_div(
272 struct isl_basic_map *bmap, unsigned div)
274 int i;
275 unsigned pos;
277 if (!bmap)
278 goto error;
280 pos = 1 + isl_space_dim(bmap->dim, isl_dim_all) + div;
282 isl_assert(bmap->ctx, div < bmap->n_div, goto error);
284 for (i = 0; i < bmap->n_eq; ++i)
285 constraint_drop_vars(bmap->eq[i]+pos, 1, bmap->extra-div-1);
287 for (i = 0; i < bmap->n_ineq; ++i) {
288 if (!isl_int_is_zero(bmap->ineq[i][pos])) {
289 isl_basic_map_drop_inequality(bmap, i);
290 --i;
291 continue;
293 constraint_drop_vars(bmap->ineq[i]+pos, 1, bmap->extra-div-1);
296 for (i = 0; i < bmap->n_div; ++i)
297 constraint_drop_vars(bmap->div[i]+1+pos, 1, bmap->extra-div-1);
299 if (div != bmap->n_div - 1) {
300 int j;
301 isl_int *t = bmap->div[div];
303 for (j = div; j < bmap->n_div - 1; ++j)
304 bmap->div[j] = bmap->div[j+1];
306 bmap->div[bmap->n_div - 1] = t;
308 ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
309 isl_basic_map_free_div(bmap, 1);
311 return bmap;
312 error:
313 isl_basic_map_free(bmap);
314 return NULL;
317 struct isl_basic_map *isl_basic_map_normalize_constraints(
318 struct isl_basic_map *bmap)
320 int i;
321 isl_int gcd;
322 unsigned total = isl_basic_map_total_dim(bmap);
324 if (!bmap)
325 return NULL;
327 isl_int_init(gcd);
328 for (i = bmap->n_eq - 1; i >= 0; --i) {
329 isl_seq_gcd(bmap->eq[i]+1, total, &gcd);
330 if (isl_int_is_zero(gcd)) {
331 if (!isl_int_is_zero(bmap->eq[i][0])) {
332 bmap = isl_basic_map_set_to_empty(bmap);
333 break;
335 isl_basic_map_drop_equality(bmap, i);
336 continue;
338 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
339 isl_int_gcd(gcd, gcd, bmap->eq[i][0]);
340 if (isl_int_is_one(gcd))
341 continue;
342 if (!isl_int_is_divisible_by(bmap->eq[i][0], gcd)) {
343 bmap = isl_basic_map_set_to_empty(bmap);
344 break;
346 isl_seq_scale_down(bmap->eq[i], bmap->eq[i], gcd, 1+total);
349 for (i = bmap->n_ineq - 1; i >= 0; --i) {
350 isl_seq_gcd(bmap->ineq[i]+1, total, &gcd);
351 if (isl_int_is_zero(gcd)) {
352 if (isl_int_is_neg(bmap->ineq[i][0])) {
353 bmap = isl_basic_map_set_to_empty(bmap);
354 break;
356 isl_basic_map_drop_inequality(bmap, i);
357 continue;
359 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
360 isl_int_gcd(gcd, gcd, bmap->ineq[i][0]);
361 if (isl_int_is_one(gcd))
362 continue;
363 isl_int_fdiv_q(bmap->ineq[i][0], bmap->ineq[i][0], gcd);
364 isl_seq_scale_down(bmap->ineq[i]+1, bmap->ineq[i]+1, gcd, total);
366 isl_int_clear(gcd);
368 return bmap;
371 struct isl_basic_set *isl_basic_set_normalize_constraints(
372 struct isl_basic_set *bset)
374 return (struct isl_basic_set *)isl_basic_map_normalize_constraints(
375 (struct isl_basic_map *)bset);
378 /* Assuming the variable at position "pos" has an integer coefficient
379 * in integer division "div", extract it from this integer division.
380 * "pos" is as determined by isl_basic_map_offset, i.e., pos == 0
381 * corresponds to the constant term.
383 * That is, the integer division is of the form
385 * floor((... + c * d * x_pos + ...)/d)
387 * Replace it by
389 * floor((... + 0 * x_pos + ...)/d) + c * x_pos
391 static __isl_give isl_basic_map *remove_var_from_div(
392 __isl_take isl_basic_map *bmap, int div, int pos)
394 isl_int shift;
396 isl_int_init(shift);
397 isl_int_divexact(shift, bmap->div[div][1 + pos], bmap->div[div][0]);
398 isl_int_neg(shift, shift);
399 bmap = isl_basic_map_shift_div(bmap, div, pos, shift);
400 isl_int_clear(shift);
402 return bmap;
405 /* Check if integer division "div" has any integral coefficient
406 * (or constant term). If so, extract them from the integer division.
408 static __isl_give isl_basic_map *remove_independent_vars_from_div(
409 __isl_take isl_basic_map *bmap, int div)
411 int i;
412 unsigned total = 1 + isl_basic_map_total_dim(bmap);
414 for (i = 0; i < total; ++i) {
415 if (isl_int_is_zero(bmap->div[div][1 + i]))
416 continue;
417 if (!isl_int_is_divisible_by(bmap->div[div][1 + i],
418 bmap->div[div][0]))
419 continue;
420 bmap = remove_var_from_div(bmap, div, i);
421 if (!bmap)
422 break;
425 return bmap;
428 /* Check if any known integer division has any integral coefficient
429 * (or constant term). If so, extract them from the integer division.
431 static __isl_give isl_basic_map *remove_independent_vars_from_divs(
432 __isl_take isl_basic_map *bmap)
434 int i;
436 if (!bmap)
437 return NULL;
438 if (bmap->n_div == 0)
439 return bmap;
441 for (i = 0; i < bmap->n_div; ++i) {
442 if (isl_int_is_zero(bmap->div[i][0]))
443 continue;
444 bmap = remove_independent_vars_from_div(bmap, i);
445 if (!bmap)
446 break;
449 return bmap;
452 /* Remove any common factor in numerator and denominator of the div expression,
453 * not taking into account the constant term.
454 * That is, if the div is of the form
456 * floor((a + m f(x))/(m d))
458 * then replace it by
460 * floor((floor(a/m) + f(x))/d)
462 * The difference {a/m}/d in the argument satisfies 0 <= {a/m}/d < 1/d
463 * and can therefore not influence the result of the floor.
465 static void normalize_div_expression(__isl_keep isl_basic_map *bmap, int div)
467 unsigned total = isl_basic_map_total_dim(bmap);
468 isl_ctx *ctx = bmap->ctx;
470 if (isl_int_is_zero(bmap->div[div][0]))
471 return;
472 isl_seq_gcd(bmap->div[div] + 2, total, &ctx->normalize_gcd);
473 isl_int_gcd(ctx->normalize_gcd, ctx->normalize_gcd, bmap->div[div][0]);
474 if (isl_int_is_one(ctx->normalize_gcd))
475 return;
476 isl_int_fdiv_q(bmap->div[div][1], bmap->div[div][1],
477 ctx->normalize_gcd);
478 isl_int_divexact(bmap->div[div][0], bmap->div[div][0],
479 ctx->normalize_gcd);
480 isl_seq_scale_down(bmap->div[div] + 2, bmap->div[div] + 2,
481 ctx->normalize_gcd, total);
484 /* Remove any common factor in numerator and denominator of a div expression,
485 * not taking into account the constant term.
486 * That is, look for any div of the form
488 * floor((a + m f(x))/(m d))
490 * and replace it by
492 * floor((floor(a/m) + f(x))/d)
494 * The difference {a/m}/d in the argument satisfies 0 <= {a/m}/d < 1/d
495 * and can therefore not influence the result of the floor.
497 static __isl_give isl_basic_map *normalize_div_expressions(
498 __isl_take isl_basic_map *bmap)
500 int i;
502 if (!bmap)
503 return NULL;
504 if (bmap->n_div == 0)
505 return bmap;
507 for (i = 0; i < bmap->n_div; ++i)
508 normalize_div_expression(bmap, i);
510 return bmap;
513 /* Assumes divs have been ordered if keep_divs is set.
515 static void eliminate_var_using_equality(struct isl_basic_map *bmap,
516 unsigned pos, isl_int *eq, int keep_divs, int *progress)
518 unsigned total;
519 unsigned space_total;
520 int k;
521 int last_div;
523 total = isl_basic_map_total_dim(bmap);
524 space_total = isl_space_dim(bmap->dim, isl_dim_all);
525 last_div = isl_seq_last_non_zero(eq + 1 + space_total, bmap->n_div);
526 for (k = 0; k < bmap->n_eq; ++k) {
527 if (bmap->eq[k] == eq)
528 continue;
529 if (isl_int_is_zero(bmap->eq[k][1+pos]))
530 continue;
531 if (progress)
532 *progress = 1;
533 isl_seq_elim(bmap->eq[k], eq, 1+pos, 1+total, NULL);
534 isl_seq_normalize(bmap->ctx, bmap->eq[k], 1 + total);
537 for (k = 0; k < bmap->n_ineq; ++k) {
538 if (isl_int_is_zero(bmap->ineq[k][1+pos]))
539 continue;
540 if (progress)
541 *progress = 1;
542 isl_seq_elim(bmap->ineq[k], eq, 1+pos, 1+total, NULL);
543 isl_seq_normalize(bmap->ctx, bmap->ineq[k], 1 + total);
544 ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
547 for (k = 0; k < bmap->n_div; ++k) {
548 if (isl_int_is_zero(bmap->div[k][0]))
549 continue;
550 if (isl_int_is_zero(bmap->div[k][1+1+pos]))
551 continue;
552 if (progress)
553 *progress = 1;
554 /* We need to be careful about circular definitions,
555 * so for now we just remove the definition of div k
556 * if the equality contains any divs.
557 * If keep_divs is set, then the divs have been ordered
558 * and we can keep the definition as long as the result
559 * is still ordered.
561 if (last_div == -1 || (keep_divs && last_div < k)) {
562 isl_seq_elim(bmap->div[k]+1, eq,
563 1+pos, 1+total, &bmap->div[k][0]);
564 normalize_div_expression(bmap, k);
565 } else
566 isl_seq_clr(bmap->div[k], 1 + total);
567 ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
571 /* Assumes divs have been ordered if keep_divs is set.
573 static __isl_give isl_basic_map *eliminate_div(__isl_take isl_basic_map *bmap,
574 isl_int *eq, unsigned div, int keep_divs)
576 unsigned pos = isl_space_dim(bmap->dim, isl_dim_all) + div;
578 eliminate_var_using_equality(bmap, pos, eq, keep_divs, NULL);
580 bmap = isl_basic_map_drop_div(bmap, div);
582 return bmap;
585 /* Check if elimination of div "div" using equality "eq" would not
586 * result in a div depending on a later div.
588 static int ok_to_eliminate_div(struct isl_basic_map *bmap, isl_int *eq,
589 unsigned div)
591 int k;
592 int last_div;
593 unsigned space_total = isl_space_dim(bmap->dim, isl_dim_all);
594 unsigned pos = space_total + div;
596 last_div = isl_seq_last_non_zero(eq + 1 + space_total, bmap->n_div);
597 if (last_div < 0 || last_div <= div)
598 return 1;
600 for (k = 0; k <= last_div; ++k) {
601 if (isl_int_is_zero(bmap->div[k][0]))
602 return 1;
603 if (!isl_int_is_zero(bmap->div[k][1 + 1 + pos]))
604 return 0;
607 return 1;
610 /* Elimininate divs based on equalities
612 static struct isl_basic_map *eliminate_divs_eq(
613 struct isl_basic_map *bmap, int *progress)
615 int d;
616 int i;
617 int modified = 0;
618 unsigned off;
620 bmap = isl_basic_map_order_divs(bmap);
622 if (!bmap)
623 return NULL;
625 off = 1 + isl_space_dim(bmap->dim, isl_dim_all);
627 for (d = bmap->n_div - 1; d >= 0 ; --d) {
628 for (i = 0; i < bmap->n_eq; ++i) {
629 if (!isl_int_is_one(bmap->eq[i][off + d]) &&
630 !isl_int_is_negone(bmap->eq[i][off + d]))
631 continue;
632 if (!ok_to_eliminate_div(bmap, bmap->eq[i], d))
633 continue;
634 modified = 1;
635 *progress = 1;
636 bmap = eliminate_div(bmap, bmap->eq[i], d, 1);
637 if (isl_basic_map_drop_equality(bmap, i) < 0)
638 return isl_basic_map_free(bmap);
639 break;
642 if (modified)
643 return eliminate_divs_eq(bmap, progress);
644 return bmap;
647 /* Elimininate divs based on inequalities
649 static struct isl_basic_map *eliminate_divs_ineq(
650 struct isl_basic_map *bmap, int *progress)
652 int d;
653 int i;
654 unsigned off;
655 struct isl_ctx *ctx;
657 if (!bmap)
658 return NULL;
660 ctx = bmap->ctx;
661 off = 1 + isl_space_dim(bmap->dim, isl_dim_all);
663 for (d = bmap->n_div - 1; d >= 0 ; --d) {
664 for (i = 0; i < bmap->n_eq; ++i)
665 if (!isl_int_is_zero(bmap->eq[i][off + d]))
666 break;
667 if (i < bmap->n_eq)
668 continue;
669 for (i = 0; i < bmap->n_ineq; ++i)
670 if (isl_int_abs_gt(bmap->ineq[i][off + d], ctx->one))
671 break;
672 if (i < bmap->n_ineq)
673 continue;
674 *progress = 1;
675 bmap = isl_basic_map_eliminate_vars(bmap, (off-1)+d, 1);
676 if (!bmap || ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
677 break;
678 bmap = isl_basic_map_drop_div(bmap, d);
679 if (!bmap)
680 break;
682 return bmap;
685 struct isl_basic_map *isl_basic_map_gauss(
686 struct isl_basic_map *bmap, int *progress)
688 int k;
689 int done;
690 int last_var;
691 unsigned total_var;
692 unsigned total;
694 bmap = isl_basic_map_order_divs(bmap);
696 if (!bmap)
697 return NULL;
699 total = isl_basic_map_total_dim(bmap);
700 total_var = total - bmap->n_div;
702 last_var = total - 1;
703 for (done = 0; done < bmap->n_eq; ++done) {
704 for (; last_var >= 0; --last_var) {
705 for (k = done; k < bmap->n_eq; ++k)
706 if (!isl_int_is_zero(bmap->eq[k][1+last_var]))
707 break;
708 if (k < bmap->n_eq)
709 break;
711 if (last_var < 0)
712 break;
713 if (k != done)
714 swap_equality(bmap, k, done);
715 if (isl_int_is_neg(bmap->eq[done][1+last_var]))
716 isl_seq_neg(bmap->eq[done], bmap->eq[done], 1+total);
718 eliminate_var_using_equality(bmap, last_var, bmap->eq[done], 1,
719 progress);
721 if (last_var >= total_var &&
722 isl_int_is_zero(bmap->div[last_var - total_var][0])) {
723 unsigned div = last_var - total_var;
724 isl_seq_neg(bmap->div[div]+1, bmap->eq[done], 1+total);
725 isl_int_set_si(bmap->div[div][1+1+last_var], 0);
726 isl_int_set(bmap->div[div][0],
727 bmap->eq[done][1+last_var]);
728 if (progress)
729 *progress = 1;
730 ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
733 if (done == bmap->n_eq)
734 return bmap;
735 for (k = done; k < bmap->n_eq; ++k) {
736 if (isl_int_is_zero(bmap->eq[k][0]))
737 continue;
738 return isl_basic_map_set_to_empty(bmap);
740 isl_basic_map_free_equality(bmap, bmap->n_eq-done);
741 return bmap;
744 struct isl_basic_set *isl_basic_set_gauss(
745 struct isl_basic_set *bset, int *progress)
747 return (struct isl_basic_set*)isl_basic_map_gauss(
748 (struct isl_basic_map *)bset, progress);
752 static unsigned int round_up(unsigned int v)
754 int old_v = v;
756 while (v) {
757 old_v = v;
758 v ^= v & -v;
760 return old_v << 1;
763 /* Hash table of inequalities in a basic map.
764 * "index" is an array of addresses of inequalities in the basic map, some
765 * of which are NULL. The inequalities are hashed on the coefficients
766 * except the constant term.
767 * "size" is the number of elements in the array and is always a power of two
768 * "bits" is the number of bits need to represent an index into the array.
769 * "total" is the total dimension of the basic map.
771 struct isl_constraint_index {
772 unsigned int size;
773 int bits;
774 isl_int ***index;
775 unsigned total;
778 /* Fill in the "ci" data structure for holding the inequalities of "bmap".
780 static isl_stat create_constraint_index(struct isl_constraint_index *ci,
781 __isl_keep isl_basic_map *bmap)
783 isl_ctx *ctx;
785 ci->index = NULL;
786 if (!bmap)
787 return isl_stat_error;
788 ci->total = isl_basic_set_total_dim(bmap);
789 if (bmap->n_ineq == 0)
790 return isl_stat_ok;
791 ci->size = round_up(4 * (bmap->n_ineq + 1) / 3 - 1);
792 ci->bits = ffs(ci->size) - 1;
793 ctx = isl_basic_map_get_ctx(bmap);
794 ci->index = isl_calloc_array(ctx, isl_int **, ci->size);
795 if (!ci->index)
796 return isl_stat_error;
798 return isl_stat_ok;
801 /* Free the memory allocated by create_constraint_index.
803 static void constraint_index_free(struct isl_constraint_index *ci)
805 free(ci->index);
808 /* Return the position in ci->index that contains the address of
809 * an inequality that is equal to *ineq up to the constant term,
810 * provided this address is not identical to "ineq".
811 * If there is no such inequality, then return the position where
812 * such an inequality should be inserted.
814 static int hash_index_ineq(struct isl_constraint_index *ci, isl_int **ineq)
816 int h;
817 uint32_t hash = isl_seq_get_hash_bits((*ineq) + 1, ci->total, ci->bits);
818 for (h = hash; ci->index[h]; h = (h+1) % ci->size)
819 if (ineq != ci->index[h] &&
820 isl_seq_eq((*ineq) + 1, ci->index[h][0]+1, ci->total))
821 break;
822 return h;
825 /* Return the position in ci->index that contains the address of
826 * an inequality that is equal to the k'th inequality of "bmap"
827 * up to the constant term, provided it does not point to the very
828 * same inequality.
829 * If there is no such inequality, then return the position where
830 * such an inequality should be inserted.
832 static int hash_index(struct isl_constraint_index *ci,
833 __isl_keep isl_basic_map *bmap, int k)
835 return hash_index_ineq(ci, &bmap->ineq[k]);
838 static int set_hash_index(struct isl_constraint_index *ci,
839 struct isl_basic_set *bset, int k)
841 return hash_index(ci, bset, k);
844 /* Fill in the "ci" data structure with the inequalities of "bset".
846 static isl_stat setup_constraint_index(struct isl_constraint_index *ci,
847 __isl_keep isl_basic_set *bset)
849 int k, h;
851 if (create_constraint_index(ci, bset) < 0)
852 return isl_stat_error;
854 for (k = 0; k < bset->n_ineq; ++k) {
855 h = set_hash_index(ci, bset, k);
856 ci->index[h] = &bset->ineq[k];
859 return isl_stat_ok;
862 /* Is the inequality ineq (obviously) redundant with respect
863 * to the constraints in "ci"?
865 * Look for an inequality in "ci" with the same coefficients and then
866 * check if the contant term of "ineq" is greater than or equal
867 * to the constant term of that inequality. If so, "ineq" is clearly
868 * redundant.
870 * Note that hash_index_ineq ignores a stored constraint if it has
871 * the same address as the passed inequality. It is ok to pass
872 * the address of a local variable here since it will never be
873 * the same as the address of a constraint in "ci".
875 static isl_bool constraint_index_is_redundant(struct isl_constraint_index *ci,
876 isl_int *ineq)
878 int h;
880 h = hash_index_ineq(ci, &ineq);
881 if (!ci->index[h])
882 return isl_bool_false;
883 return isl_int_ge(ineq[0], (*ci->index[h])[0]);
886 /* If we can eliminate more than one div, then we need to make
887 * sure we do it from last div to first div, in order not to
888 * change the position of the other divs that still need to
889 * be removed.
891 static struct isl_basic_map *remove_duplicate_divs(
892 struct isl_basic_map *bmap, int *progress)
894 unsigned int size;
895 int *index;
896 int *elim_for;
897 int k, l, h;
898 int bits;
899 struct isl_blk eq;
900 unsigned total_var;
901 unsigned total;
902 struct isl_ctx *ctx;
904 bmap = isl_basic_map_order_divs(bmap);
905 if (!bmap || bmap->n_div <= 1)
906 return bmap;
908 total_var = isl_space_dim(bmap->dim, isl_dim_all);
909 total = total_var + bmap->n_div;
911 ctx = bmap->ctx;
912 for (k = bmap->n_div - 1; k >= 0; --k)
913 if (!isl_int_is_zero(bmap->div[k][0]))
914 break;
915 if (k <= 0)
916 return bmap;
918 size = round_up(4 * bmap->n_div / 3 - 1);
919 if (size == 0)
920 return bmap;
921 elim_for = isl_calloc_array(ctx, int, bmap->n_div);
922 bits = ffs(size) - 1;
923 index = isl_calloc_array(ctx, int, size);
924 if (!elim_for || !index)
925 goto out;
926 eq = isl_blk_alloc(ctx, 1+total);
927 if (isl_blk_is_error(eq))
928 goto out;
930 isl_seq_clr(eq.data, 1+total);
931 index[isl_seq_get_hash_bits(bmap->div[k], 2+total, bits)] = k + 1;
932 for (--k; k >= 0; --k) {
933 uint32_t hash;
935 if (isl_int_is_zero(bmap->div[k][0]))
936 continue;
938 hash = isl_seq_get_hash_bits(bmap->div[k], 2+total, bits);
939 for (h = hash; index[h]; h = (h+1) % size)
940 if (isl_seq_eq(bmap->div[k],
941 bmap->div[index[h]-1], 2+total))
942 break;
943 if (index[h]) {
944 *progress = 1;
945 l = index[h] - 1;
946 elim_for[l] = k + 1;
948 index[h] = k+1;
950 for (l = bmap->n_div - 1; l >= 0; --l) {
951 if (!elim_for[l])
952 continue;
953 k = elim_for[l] - 1;
954 isl_int_set_si(eq.data[1+total_var+k], -1);
955 isl_int_set_si(eq.data[1+total_var+l], 1);
956 bmap = eliminate_div(bmap, eq.data, l, 1);
957 if (!bmap)
958 break;
959 isl_int_set_si(eq.data[1+total_var+k], 0);
960 isl_int_set_si(eq.data[1+total_var+l], 0);
963 isl_blk_free(ctx, eq);
964 out:
965 free(index);
966 free(elim_for);
967 return bmap;
970 static int n_pure_div_eq(struct isl_basic_map *bmap)
972 int i, j;
973 unsigned total;
975 total = isl_space_dim(bmap->dim, isl_dim_all);
976 for (i = 0, j = bmap->n_div-1; i < bmap->n_eq; ++i) {
977 while (j >= 0 && isl_int_is_zero(bmap->eq[i][1 + total + j]))
978 --j;
979 if (j < 0)
980 break;
981 if (isl_seq_first_non_zero(bmap->eq[i] + 1 + total, j) != -1)
982 return 0;
984 return i;
987 /* Normalize divs that appear in equalities.
989 * In particular, we assume that bmap contains some equalities
990 * of the form
992 * a x = m * e_i
994 * and we want to replace the set of e_i by a minimal set and
995 * such that the new e_i have a canonical representation in terms
996 * of the vector x.
997 * If any of the equalities involves more than one divs, then
998 * we currently simply bail out.
1000 * Let us first additionally assume that all equalities involve
1001 * a div. The equalities then express modulo constraints on the
1002 * remaining variables and we can use "parameter compression"
1003 * to find a minimal set of constraints. The result is a transformation
1005 * x = T(x') = x_0 + G x'
1007 * with G a lower-triangular matrix with all elements below the diagonal
1008 * non-negative and smaller than the diagonal element on the same row.
1009 * We first normalize x_0 by making the same property hold in the affine
1010 * T matrix.
1011 * The rows i of G with a 1 on the diagonal do not impose any modulo
1012 * constraint and simply express x_i = x'_i.
1013 * For each of the remaining rows i, we introduce a div and a corresponding
1014 * equality. In particular
1016 * g_ii e_j = x_i - g_i(x')
1018 * where each x'_k is replaced either by x_k (if g_kk = 1) or the
1019 * corresponding div (if g_kk != 1).
1021 * If there are any equalities not involving any div, then we
1022 * first apply a variable compression on the variables x:
1024 * x = C x'' x'' = C_2 x
1026 * and perform the above parameter compression on A C instead of on A.
1027 * The resulting compression is then of the form
1029 * x'' = T(x') = x_0 + G x'
1031 * and in constructing the new divs and the corresponding equalities,
1032 * we have to replace each x'', i.e., the x'_k with (g_kk = 1),
1033 * by the corresponding row from C_2.
1035 static struct isl_basic_map *normalize_divs(
1036 struct isl_basic_map *bmap, int *progress)
1038 int i, j, k;
1039 int total;
1040 int div_eq;
1041 struct isl_mat *B;
1042 struct isl_vec *d;
1043 struct isl_mat *T = NULL;
1044 struct isl_mat *C = NULL;
1045 struct isl_mat *C2 = NULL;
1046 isl_int v;
1047 int *pos;
1048 int dropped, needed;
1050 if (!bmap)
1051 return NULL;
1053 if (bmap->n_div == 0)
1054 return bmap;
1056 if (bmap->n_eq == 0)
1057 return bmap;
1059 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_NORMALIZED_DIVS))
1060 return bmap;
1062 total = isl_space_dim(bmap->dim, isl_dim_all);
1063 div_eq = n_pure_div_eq(bmap);
1064 if (div_eq == 0)
1065 return bmap;
1067 if (div_eq < bmap->n_eq) {
1068 B = isl_mat_sub_alloc6(bmap->ctx, bmap->eq, div_eq,
1069 bmap->n_eq - div_eq, 0, 1 + total);
1070 C = isl_mat_variable_compression(B, &C2);
1071 if (!C || !C2)
1072 goto error;
1073 if (C->n_col == 0) {
1074 bmap = isl_basic_map_set_to_empty(bmap);
1075 isl_mat_free(C);
1076 isl_mat_free(C2);
1077 goto done;
1081 d = isl_vec_alloc(bmap->ctx, div_eq);
1082 if (!d)
1083 goto error;
1084 for (i = 0, j = bmap->n_div-1; i < div_eq; ++i) {
1085 while (j >= 0 && isl_int_is_zero(bmap->eq[i][1 + total + j]))
1086 --j;
1087 isl_int_set(d->block.data[i], bmap->eq[i][1 + total + j]);
1089 B = isl_mat_sub_alloc6(bmap->ctx, bmap->eq, 0, div_eq, 0, 1 + total);
1091 if (C) {
1092 B = isl_mat_product(B, C);
1093 C = NULL;
1096 T = isl_mat_parameter_compression(B, d);
1097 if (!T)
1098 goto error;
1099 if (T->n_col == 0) {
1100 bmap = isl_basic_map_set_to_empty(bmap);
1101 isl_mat_free(C2);
1102 isl_mat_free(T);
1103 goto done;
1105 isl_int_init(v);
1106 for (i = 0; i < T->n_row - 1; ++i) {
1107 isl_int_fdiv_q(v, T->row[1 + i][0], T->row[1 + i][1 + i]);
1108 if (isl_int_is_zero(v))
1109 continue;
1110 isl_mat_col_submul(T, 0, v, 1 + i);
1112 isl_int_clear(v);
1113 pos = isl_alloc_array(bmap->ctx, int, T->n_row);
1114 if (!pos)
1115 goto error;
1116 /* We have to be careful because dropping equalities may reorder them */
1117 dropped = 0;
1118 for (j = bmap->n_div - 1; j >= 0; --j) {
1119 for (i = 0; i < bmap->n_eq; ++i)
1120 if (!isl_int_is_zero(bmap->eq[i][1 + total + j]))
1121 break;
1122 if (i < bmap->n_eq) {
1123 bmap = isl_basic_map_drop_div(bmap, j);
1124 isl_basic_map_drop_equality(bmap, i);
1125 ++dropped;
1128 pos[0] = 0;
1129 needed = 0;
1130 for (i = 1; i < T->n_row; ++i) {
1131 if (isl_int_is_one(T->row[i][i]))
1132 pos[i] = i;
1133 else
1134 needed++;
1136 if (needed > dropped) {
1137 bmap = isl_basic_map_extend_space(bmap, isl_space_copy(bmap->dim),
1138 needed, needed, 0);
1139 if (!bmap)
1140 goto error;
1142 for (i = 1; i < T->n_row; ++i) {
1143 if (isl_int_is_one(T->row[i][i]))
1144 continue;
1145 k = isl_basic_map_alloc_div(bmap);
1146 pos[i] = 1 + total + k;
1147 isl_seq_clr(bmap->div[k] + 1, 1 + total + bmap->n_div);
1148 isl_int_set(bmap->div[k][0], T->row[i][i]);
1149 if (C2)
1150 isl_seq_cpy(bmap->div[k] + 1, C2->row[i], 1 + total);
1151 else
1152 isl_int_set_si(bmap->div[k][1 + i], 1);
1153 for (j = 0; j < i; ++j) {
1154 if (isl_int_is_zero(T->row[i][j]))
1155 continue;
1156 if (pos[j] < T->n_row && C2)
1157 isl_seq_submul(bmap->div[k] + 1, T->row[i][j],
1158 C2->row[pos[j]], 1 + total);
1159 else
1160 isl_int_neg(bmap->div[k][1 + pos[j]],
1161 T->row[i][j]);
1163 j = isl_basic_map_alloc_equality(bmap);
1164 isl_seq_neg(bmap->eq[j], bmap->div[k]+1, 1+total+bmap->n_div);
1165 isl_int_set(bmap->eq[j][pos[i]], bmap->div[k][0]);
1167 free(pos);
1168 isl_mat_free(C2);
1169 isl_mat_free(T);
1171 if (progress)
1172 *progress = 1;
1173 done:
1174 ISL_F_SET(bmap, ISL_BASIC_MAP_NORMALIZED_DIVS);
1176 return bmap;
1177 error:
1178 isl_mat_free(C);
1179 isl_mat_free(C2);
1180 isl_mat_free(T);
1181 return bmap;
1184 static struct isl_basic_map *set_div_from_lower_bound(
1185 struct isl_basic_map *bmap, int div, int ineq)
1187 unsigned total = 1 + isl_space_dim(bmap->dim, isl_dim_all);
1189 isl_seq_neg(bmap->div[div] + 1, bmap->ineq[ineq], total + bmap->n_div);
1190 isl_int_set(bmap->div[div][0], bmap->ineq[ineq][total + div]);
1191 isl_int_add(bmap->div[div][1], bmap->div[div][1], bmap->div[div][0]);
1192 isl_int_sub_ui(bmap->div[div][1], bmap->div[div][1], 1);
1193 isl_int_set_si(bmap->div[div][1 + total + div], 0);
1195 return bmap;
1198 /* Check whether it is ok to define a div based on an inequality.
1199 * To avoid the introduction of circular definitions of divs, we
1200 * do not allow such a definition if the resulting expression would refer to
1201 * any other undefined divs or if any known div is defined in
1202 * terms of the unknown div.
1204 static int ok_to_set_div_from_bound(struct isl_basic_map *bmap,
1205 int div, int ineq)
1207 int j;
1208 unsigned total = 1 + isl_space_dim(bmap->dim, isl_dim_all);
1210 /* Not defined in terms of unknown divs */
1211 for (j = 0; j < bmap->n_div; ++j) {
1212 if (div == j)
1213 continue;
1214 if (isl_int_is_zero(bmap->ineq[ineq][total + j]))
1215 continue;
1216 if (isl_int_is_zero(bmap->div[j][0]))
1217 return 0;
1220 /* No other div defined in terms of this one => avoid loops */
1221 for (j = 0; j < bmap->n_div; ++j) {
1222 if (div == j)
1223 continue;
1224 if (isl_int_is_zero(bmap->div[j][0]))
1225 continue;
1226 if (!isl_int_is_zero(bmap->div[j][1 + total + div]))
1227 return 0;
1230 return 1;
1233 /* Would an expression for div "div" based on inequality "ineq" of "bmap"
1234 * be a better expression than the current one?
1236 * If we do not have any expression yet, then any expression would be better.
1237 * Otherwise we check if the last variable involved in the inequality
1238 * (disregarding the div that it would define) is in an earlier position
1239 * than the last variable involved in the current div expression.
1241 static int better_div_constraint(__isl_keep isl_basic_map *bmap,
1242 int div, int ineq)
1244 unsigned total = 1 + isl_space_dim(bmap->dim, isl_dim_all);
1245 int last_div;
1246 int last_ineq;
1248 if (isl_int_is_zero(bmap->div[div][0]))
1249 return 1;
1251 if (isl_seq_last_non_zero(bmap->ineq[ineq] + total + div + 1,
1252 bmap->n_div - (div + 1)) >= 0)
1253 return 0;
1255 last_ineq = isl_seq_last_non_zero(bmap->ineq[ineq], total + div);
1256 last_div = isl_seq_last_non_zero(bmap->div[div] + 1,
1257 total + bmap->n_div);
1259 return last_ineq < last_div;
1262 /* Given two constraints "k" and "l" that are opposite to each other,
1263 * except for the constant term, check if we can use them
1264 * to obtain an expression for one of the hitherto unknown divs or
1265 * a "better" expression for a div for which we already have an expression.
1266 * "sum" is the sum of the constant terms of the constraints.
1267 * If this sum is strictly smaller than the coefficient of one
1268 * of the divs, then this pair can be used define the div.
1269 * To avoid the introduction of circular definitions of divs, we
1270 * do not use the pair if the resulting expression would refer to
1271 * any other undefined divs or if any known div is defined in
1272 * terms of the unknown div.
1274 static struct isl_basic_map *check_for_div_constraints(
1275 struct isl_basic_map *bmap, int k, int l, isl_int sum, int *progress)
1277 int i;
1278 unsigned total = 1 + isl_space_dim(bmap->dim, isl_dim_all);
1280 for (i = 0; i < bmap->n_div; ++i) {
1281 if (isl_int_is_zero(bmap->ineq[k][total + i]))
1282 continue;
1283 if (isl_int_abs_ge(sum, bmap->ineq[k][total + i]))
1284 continue;
1285 if (!better_div_constraint(bmap, i, k))
1286 continue;
1287 if (!ok_to_set_div_from_bound(bmap, i, k))
1288 break;
1289 if (isl_int_is_pos(bmap->ineq[k][total + i]))
1290 bmap = set_div_from_lower_bound(bmap, i, k);
1291 else
1292 bmap = set_div_from_lower_bound(bmap, i, l);
1293 if (progress)
1294 *progress = 1;
1295 break;
1297 return bmap;
1300 __isl_give isl_basic_map *isl_basic_map_remove_duplicate_constraints(
1301 __isl_take isl_basic_map *bmap, int *progress, int detect_divs)
1303 struct isl_constraint_index ci;
1304 int k, l, h;
1305 unsigned total = isl_basic_map_total_dim(bmap);
1306 isl_int sum;
1308 if (!bmap || bmap->n_ineq <= 1)
1309 return bmap;
1311 if (create_constraint_index(&ci, bmap) < 0)
1312 return bmap;
1314 h = isl_seq_get_hash_bits(bmap->ineq[0] + 1, total, ci.bits);
1315 ci.index[h] = &bmap->ineq[0];
1316 for (k = 1; k < bmap->n_ineq; ++k) {
1317 h = hash_index(&ci, bmap, k);
1318 if (!ci.index[h]) {
1319 ci.index[h] = &bmap->ineq[k];
1320 continue;
1322 if (progress)
1323 *progress = 1;
1324 l = ci.index[h] - &bmap->ineq[0];
1325 if (isl_int_lt(bmap->ineq[k][0], bmap->ineq[l][0]))
1326 swap_inequality(bmap, k, l);
1327 isl_basic_map_drop_inequality(bmap, k);
1328 --k;
1330 isl_int_init(sum);
1331 for (k = 0; k < bmap->n_ineq-1; ++k) {
1332 isl_seq_neg(bmap->ineq[k]+1, bmap->ineq[k]+1, total);
1333 h = hash_index(&ci, bmap, k);
1334 isl_seq_neg(bmap->ineq[k]+1, bmap->ineq[k]+1, total);
1335 if (!ci.index[h])
1336 continue;
1337 l = ci.index[h] - &bmap->ineq[0];
1338 isl_int_add(sum, bmap->ineq[k][0], bmap->ineq[l][0]);
1339 if (isl_int_is_pos(sum)) {
1340 if (detect_divs)
1341 bmap = check_for_div_constraints(bmap, k, l,
1342 sum, progress);
1343 continue;
1345 if (isl_int_is_zero(sum)) {
1346 /* We need to break out of the loop after these
1347 * changes since the contents of the hash
1348 * will no longer be valid.
1349 * Plus, we probably we want to regauss first.
1351 if (progress)
1352 *progress = 1;
1353 isl_basic_map_drop_inequality(bmap, l);
1354 isl_basic_map_inequality_to_equality(bmap, k);
1355 } else
1356 bmap = isl_basic_map_set_to_empty(bmap);
1357 break;
1359 isl_int_clear(sum);
1361 constraint_index_free(&ci);
1362 return bmap;
1365 /* Detect all pairs of inequalities that form an equality.
1367 * isl_basic_map_remove_duplicate_constraints detects at most one such pair.
1368 * Call it repeatedly while it is making progress.
1370 __isl_give isl_basic_map *isl_basic_map_detect_inequality_pairs(
1371 __isl_take isl_basic_map *bmap, int *progress)
1373 int duplicate;
1375 do {
1376 duplicate = 0;
1377 bmap = isl_basic_map_remove_duplicate_constraints(bmap,
1378 &duplicate, 0);
1379 if (progress && duplicate)
1380 *progress = 1;
1381 } while (duplicate);
1383 return bmap;
1386 /* Eliminate knowns divs from constraints where they appear with
1387 * a (positive or negative) unit coefficient.
1389 * That is, replace
1391 * floor(e/m) + f >= 0
1393 * by
1395 * e + m f >= 0
1397 * and
1399 * -floor(e/m) + f >= 0
1401 * by
1403 * -e + m f + m - 1 >= 0
1405 * The first conversion is valid because floor(e/m) >= -f is equivalent
1406 * to e/m >= -f because -f is an integral expression.
1407 * The second conversion follows from the fact that
1409 * -floor(e/m) = ceil(-e/m) = floor((-e + m - 1)/m)
1412 * Note that one of the div constraints may have been eliminated
1413 * due to being redundant with respect to the constraint that is
1414 * being modified by this function. The modified constraint may
1415 * no longer imply this div constraint, so we add it back to make
1416 * sure we do not lose any information.
1418 * We skip integral divs, i.e., those with denominator 1, as we would
1419 * risk eliminating the div from the div constraints. We do not need
1420 * to handle those divs here anyway since the div constraints will turn
1421 * out to form an equality and this equality can then be use to eliminate
1422 * the div from all constraints.
1424 static __isl_give isl_basic_map *eliminate_unit_divs(
1425 __isl_take isl_basic_map *bmap, int *progress)
1427 int i, j;
1428 isl_ctx *ctx;
1429 unsigned total;
1431 if (!bmap)
1432 return NULL;
1434 ctx = isl_basic_map_get_ctx(bmap);
1435 total = 1 + isl_space_dim(bmap->dim, isl_dim_all);
1437 for (i = 0; i < bmap->n_div; ++i) {
1438 if (isl_int_is_zero(bmap->div[i][0]))
1439 continue;
1440 if (isl_int_is_one(bmap->div[i][0]))
1441 continue;
1442 for (j = 0; j < bmap->n_ineq; ++j) {
1443 int s;
1445 if (!isl_int_is_one(bmap->ineq[j][total + i]) &&
1446 !isl_int_is_negone(bmap->ineq[j][total + i]))
1447 continue;
1449 *progress = 1;
1451 s = isl_int_sgn(bmap->ineq[j][total + i]);
1452 isl_int_set_si(bmap->ineq[j][total + i], 0);
1453 if (s < 0)
1454 isl_seq_combine(bmap->ineq[j],
1455 ctx->negone, bmap->div[i] + 1,
1456 bmap->div[i][0], bmap->ineq[j],
1457 total + bmap->n_div);
1458 else
1459 isl_seq_combine(bmap->ineq[j],
1460 ctx->one, bmap->div[i] + 1,
1461 bmap->div[i][0], bmap->ineq[j],
1462 total + bmap->n_div);
1463 if (s < 0) {
1464 isl_int_add(bmap->ineq[j][0],
1465 bmap->ineq[j][0], bmap->div[i][0]);
1466 isl_int_sub_ui(bmap->ineq[j][0],
1467 bmap->ineq[j][0], 1);
1470 bmap = isl_basic_map_extend_constraints(bmap, 0, 1);
1471 if (isl_basic_map_add_div_constraint(bmap, i, s) < 0)
1472 return isl_basic_map_free(bmap);
1476 return bmap;
1479 struct isl_basic_map *isl_basic_map_simplify(struct isl_basic_map *bmap)
1481 int progress = 1;
1482 if (!bmap)
1483 return NULL;
1484 while (progress) {
1485 progress = 0;
1486 if (!bmap)
1487 break;
1488 if (isl_basic_map_plain_is_empty(bmap))
1489 break;
1490 bmap = isl_basic_map_normalize_constraints(bmap);
1491 bmap = remove_independent_vars_from_divs(bmap);
1492 bmap = normalize_div_expressions(bmap);
1493 bmap = remove_duplicate_divs(bmap, &progress);
1494 bmap = eliminate_unit_divs(bmap, &progress);
1495 bmap = eliminate_divs_eq(bmap, &progress);
1496 bmap = eliminate_divs_ineq(bmap, &progress);
1497 bmap = isl_basic_map_gauss(bmap, &progress);
1498 /* requires equalities in normal form */
1499 bmap = normalize_divs(bmap, &progress);
1500 bmap = isl_basic_map_remove_duplicate_constraints(bmap,
1501 &progress, 1);
1502 if (bmap && progress)
1503 ISL_F_CLR(bmap, ISL_BASIC_MAP_REDUCED_COEFFICIENTS);
1505 return bmap;
1508 struct isl_basic_set *isl_basic_set_simplify(struct isl_basic_set *bset)
1510 return (struct isl_basic_set *)
1511 isl_basic_map_simplify((struct isl_basic_map *)bset);
1515 int isl_basic_map_is_div_constraint(__isl_keep isl_basic_map *bmap,
1516 isl_int *constraint, unsigned div)
1518 unsigned pos;
1520 if (!bmap)
1521 return -1;
1523 pos = 1 + isl_space_dim(bmap->dim, isl_dim_all) + div;
1525 if (isl_int_eq(constraint[pos], bmap->div[div][0])) {
1526 int neg;
1527 isl_int_sub(bmap->div[div][1],
1528 bmap->div[div][1], bmap->div[div][0]);
1529 isl_int_add_ui(bmap->div[div][1], bmap->div[div][1], 1);
1530 neg = isl_seq_is_neg(constraint, bmap->div[div]+1, pos);
1531 isl_int_sub_ui(bmap->div[div][1], bmap->div[div][1], 1);
1532 isl_int_add(bmap->div[div][1],
1533 bmap->div[div][1], bmap->div[div][0]);
1534 if (!neg)
1535 return 0;
1536 if (isl_seq_first_non_zero(constraint+pos+1,
1537 bmap->n_div-div-1) != -1)
1538 return 0;
1539 } else if (isl_int_abs_eq(constraint[pos], bmap->div[div][0])) {
1540 if (!isl_seq_eq(constraint, bmap->div[div]+1, pos))
1541 return 0;
1542 if (isl_seq_first_non_zero(constraint+pos+1,
1543 bmap->n_div-div-1) != -1)
1544 return 0;
1545 } else
1546 return 0;
1548 return 1;
1551 int isl_basic_set_is_div_constraint(__isl_keep isl_basic_set *bset,
1552 isl_int *constraint, unsigned div)
1554 return isl_basic_map_is_div_constraint(bset, constraint, div);
1558 /* If the only constraints a div d=floor(f/m)
1559 * appears in are its two defining constraints
1561 * f - m d >=0
1562 * -(f - (m - 1)) + m d >= 0
1564 * then it can safely be removed.
1566 static int div_is_redundant(struct isl_basic_map *bmap, int div)
1568 int i;
1569 unsigned pos = 1 + isl_space_dim(bmap->dim, isl_dim_all) + div;
1571 for (i = 0; i < bmap->n_eq; ++i)
1572 if (!isl_int_is_zero(bmap->eq[i][pos]))
1573 return 0;
1575 for (i = 0; i < bmap->n_ineq; ++i) {
1576 if (isl_int_is_zero(bmap->ineq[i][pos]))
1577 continue;
1578 if (!isl_basic_map_is_div_constraint(bmap, bmap->ineq[i], div))
1579 return 0;
1582 for (i = 0; i < bmap->n_div; ++i) {
1583 if (isl_int_is_zero(bmap->div[i][0]))
1584 continue;
1585 if (!isl_int_is_zero(bmap->div[i][1+pos]))
1586 return 0;
1589 return 1;
1593 * Remove divs that don't occur in any of the constraints or other divs.
1594 * These can arise when dropping constraints from a basic map or
1595 * when the divs of a basic map have been temporarily aligned
1596 * with the divs of another basic map.
1598 static struct isl_basic_map *remove_redundant_divs(struct isl_basic_map *bmap)
1600 int i;
1602 if (!bmap)
1603 return NULL;
1605 for (i = bmap->n_div-1; i >= 0; --i) {
1606 if (!div_is_redundant(bmap, i))
1607 continue;
1608 bmap = isl_basic_map_drop_div(bmap, i);
1610 return bmap;
1613 struct isl_basic_map *isl_basic_map_finalize(struct isl_basic_map *bmap)
1615 bmap = remove_redundant_divs(bmap);
1616 if (!bmap)
1617 return NULL;
1618 ISL_F_SET(bmap, ISL_BASIC_SET_FINAL);
1619 return bmap;
1622 struct isl_basic_set *isl_basic_set_finalize(struct isl_basic_set *bset)
1624 return (struct isl_basic_set *)
1625 isl_basic_map_finalize((struct isl_basic_map *)bset);
1628 struct isl_set *isl_set_finalize(struct isl_set *set)
1630 int i;
1632 if (!set)
1633 return NULL;
1634 for (i = 0; i < set->n; ++i) {
1635 set->p[i] = isl_basic_set_finalize(set->p[i]);
1636 if (!set->p[i])
1637 goto error;
1639 return set;
1640 error:
1641 isl_set_free(set);
1642 return NULL;
1645 struct isl_map *isl_map_finalize(struct isl_map *map)
1647 int i;
1649 if (!map)
1650 return NULL;
1651 for (i = 0; i < map->n; ++i) {
1652 map->p[i] = isl_basic_map_finalize(map->p[i]);
1653 if (!map->p[i])
1654 goto error;
1656 ISL_F_CLR(map, ISL_MAP_NORMALIZED);
1657 return map;
1658 error:
1659 isl_map_free(map);
1660 return NULL;
1664 /* Remove definition of any div that is defined in terms of the given variable.
1665 * The div itself is not removed. Functions such as
1666 * eliminate_divs_ineq depend on the other divs remaining in place.
1668 static struct isl_basic_map *remove_dependent_vars(struct isl_basic_map *bmap,
1669 int pos)
1671 int i;
1673 if (!bmap)
1674 return NULL;
1676 for (i = 0; i < bmap->n_div; ++i) {
1677 if (isl_int_is_zero(bmap->div[i][0]))
1678 continue;
1679 if (isl_int_is_zero(bmap->div[i][1+1+pos]))
1680 continue;
1681 isl_int_set_si(bmap->div[i][0], 0);
1683 return bmap;
1686 /* Eliminate the specified variables from the constraints using
1687 * Fourier-Motzkin. The variables themselves are not removed.
1689 struct isl_basic_map *isl_basic_map_eliminate_vars(
1690 struct isl_basic_map *bmap, unsigned pos, unsigned n)
1692 int d;
1693 int i, j, k;
1694 unsigned total;
1695 int need_gauss = 0;
1697 if (n == 0)
1698 return bmap;
1699 if (!bmap)
1700 return NULL;
1701 total = isl_basic_map_total_dim(bmap);
1703 bmap = isl_basic_map_cow(bmap);
1704 for (d = pos + n - 1; d >= 0 && d >= pos; --d)
1705 bmap = remove_dependent_vars(bmap, d);
1706 if (!bmap)
1707 return NULL;
1709 for (d = pos + n - 1;
1710 d >= 0 && d >= total - bmap->n_div && d >= pos; --d)
1711 isl_seq_clr(bmap->div[d-(total-bmap->n_div)], 2+total);
1712 for (d = pos + n - 1; d >= 0 && d >= pos; --d) {
1713 int n_lower, n_upper;
1714 if (!bmap)
1715 return NULL;
1716 for (i = 0; i < bmap->n_eq; ++i) {
1717 if (isl_int_is_zero(bmap->eq[i][1+d]))
1718 continue;
1719 eliminate_var_using_equality(bmap, d, bmap->eq[i], 0, NULL);
1720 isl_basic_map_drop_equality(bmap, i);
1721 need_gauss = 1;
1722 break;
1724 if (i < bmap->n_eq)
1725 continue;
1726 n_lower = 0;
1727 n_upper = 0;
1728 for (i = 0; i < bmap->n_ineq; ++i) {
1729 if (isl_int_is_pos(bmap->ineq[i][1+d]))
1730 n_lower++;
1731 else if (isl_int_is_neg(bmap->ineq[i][1+d]))
1732 n_upper++;
1734 bmap = isl_basic_map_extend_constraints(bmap,
1735 0, n_lower * n_upper);
1736 if (!bmap)
1737 goto error;
1738 for (i = bmap->n_ineq - 1; i >= 0; --i) {
1739 int last;
1740 if (isl_int_is_zero(bmap->ineq[i][1+d]))
1741 continue;
1742 last = -1;
1743 for (j = 0; j < i; ++j) {
1744 if (isl_int_is_zero(bmap->ineq[j][1+d]))
1745 continue;
1746 last = j;
1747 if (isl_int_sgn(bmap->ineq[i][1+d]) ==
1748 isl_int_sgn(bmap->ineq[j][1+d]))
1749 continue;
1750 k = isl_basic_map_alloc_inequality(bmap);
1751 if (k < 0)
1752 goto error;
1753 isl_seq_cpy(bmap->ineq[k], bmap->ineq[i],
1754 1+total);
1755 isl_seq_elim(bmap->ineq[k], bmap->ineq[j],
1756 1+d, 1+total, NULL);
1758 isl_basic_map_drop_inequality(bmap, i);
1759 i = last + 1;
1761 if (n_lower > 0 && n_upper > 0) {
1762 bmap = isl_basic_map_normalize_constraints(bmap);
1763 bmap = isl_basic_map_remove_duplicate_constraints(bmap,
1764 NULL, 0);
1765 bmap = isl_basic_map_gauss(bmap, NULL);
1766 bmap = isl_basic_map_remove_redundancies(bmap);
1767 need_gauss = 0;
1768 if (!bmap)
1769 goto error;
1770 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
1771 break;
1774 ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
1775 if (need_gauss)
1776 bmap = isl_basic_map_gauss(bmap, NULL);
1777 return bmap;
1778 error:
1779 isl_basic_map_free(bmap);
1780 return NULL;
1783 struct isl_basic_set *isl_basic_set_eliminate_vars(
1784 struct isl_basic_set *bset, unsigned pos, unsigned n)
1786 return (struct isl_basic_set *)isl_basic_map_eliminate_vars(
1787 (struct isl_basic_map *)bset, pos, n);
1790 /* Eliminate the specified n dimensions starting at first from the
1791 * constraints, without removing the dimensions from the space.
1792 * If the set is rational, the dimensions are eliminated using Fourier-Motzkin.
1793 * Otherwise, they are projected out and the original space is restored.
1795 __isl_give isl_basic_map *isl_basic_map_eliminate(
1796 __isl_take isl_basic_map *bmap,
1797 enum isl_dim_type type, unsigned first, unsigned n)
1799 isl_space *space;
1801 if (!bmap)
1802 return NULL;
1803 if (n == 0)
1804 return bmap;
1806 if (first + n > isl_basic_map_dim(bmap, type) || first + n < first)
1807 isl_die(bmap->ctx, isl_error_invalid,
1808 "index out of bounds", goto error);
1810 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL)) {
1811 first += isl_basic_map_offset(bmap, type) - 1;
1812 bmap = isl_basic_map_eliminate_vars(bmap, first, n);
1813 return isl_basic_map_finalize(bmap);
1816 space = isl_basic_map_get_space(bmap);
1817 bmap = isl_basic_map_project_out(bmap, type, first, n);
1818 bmap = isl_basic_map_insert_dims(bmap, type, first, n);
1819 bmap = isl_basic_map_reset_space(bmap, space);
1820 return bmap;
1821 error:
1822 isl_basic_map_free(bmap);
1823 return NULL;
1826 __isl_give isl_basic_set *isl_basic_set_eliminate(
1827 __isl_take isl_basic_set *bset,
1828 enum isl_dim_type type, unsigned first, unsigned n)
1830 return isl_basic_map_eliminate(bset, type, first, n);
1833 /* Don't assume equalities are in order, because align_divs
1834 * may have changed the order of the divs.
1836 static void compute_elimination_index(struct isl_basic_map *bmap, int *elim)
1838 int d, i;
1839 unsigned total;
1841 total = isl_space_dim(bmap->dim, isl_dim_all);
1842 for (d = 0; d < total; ++d)
1843 elim[d] = -1;
1844 for (i = 0; i < bmap->n_eq; ++i) {
1845 for (d = total - 1; d >= 0; --d) {
1846 if (isl_int_is_zero(bmap->eq[i][1+d]))
1847 continue;
1848 elim[d] = i;
1849 break;
1854 static void set_compute_elimination_index(struct isl_basic_set *bset, int *elim)
1856 compute_elimination_index((struct isl_basic_map *)bset, elim);
1859 static int reduced_using_equalities(isl_int *dst, isl_int *src,
1860 struct isl_basic_map *bmap, int *elim)
1862 int d;
1863 int copied = 0;
1864 unsigned total;
1866 total = isl_space_dim(bmap->dim, isl_dim_all);
1867 for (d = total - 1; d >= 0; --d) {
1868 if (isl_int_is_zero(src[1+d]))
1869 continue;
1870 if (elim[d] == -1)
1871 continue;
1872 if (!copied) {
1873 isl_seq_cpy(dst, src, 1 + total);
1874 copied = 1;
1876 isl_seq_elim(dst, bmap->eq[elim[d]], 1 + d, 1 + total, NULL);
1878 return copied;
1881 static int set_reduced_using_equalities(isl_int *dst, isl_int *src,
1882 struct isl_basic_set *bset, int *elim)
1884 return reduced_using_equalities(dst, src,
1885 (struct isl_basic_map *)bset, elim);
1888 static struct isl_basic_set *isl_basic_set_reduce_using_equalities(
1889 struct isl_basic_set *bset, struct isl_basic_set *context)
1891 int i;
1892 int *elim;
1894 if (!bset || !context)
1895 goto error;
1897 if (context->n_eq == 0) {
1898 isl_basic_set_free(context);
1899 return bset;
1902 bset = isl_basic_set_cow(bset);
1903 if (!bset)
1904 goto error;
1906 elim = isl_alloc_array(bset->ctx, int, isl_basic_set_n_dim(bset));
1907 if (!elim)
1908 goto error;
1909 set_compute_elimination_index(context, elim);
1910 for (i = 0; i < bset->n_eq; ++i)
1911 set_reduced_using_equalities(bset->eq[i], bset->eq[i],
1912 context, elim);
1913 for (i = 0; i < bset->n_ineq; ++i)
1914 set_reduced_using_equalities(bset->ineq[i], bset->ineq[i],
1915 context, elim);
1916 isl_basic_set_free(context);
1917 free(elim);
1918 bset = isl_basic_set_simplify(bset);
1919 bset = isl_basic_set_finalize(bset);
1920 return bset;
1921 error:
1922 isl_basic_set_free(bset);
1923 isl_basic_set_free(context);
1924 return NULL;
1927 /* For each inequality in "bset" that is a shifted (more relaxed)
1928 * copy of an inequality in "context", mark the corresponding entry
1929 * in "row" with -1.
1931 static isl_stat mark_shifted_constraints(__isl_keep isl_basic_set *bset,
1932 __isl_keep isl_basic_set *context, int *row)
1934 struct isl_constraint_index ci;
1935 unsigned total;
1936 int k;
1938 if (!bset || !context)
1939 return isl_stat_error;
1940 if (context->n_ineq == 0)
1941 return isl_stat_ok;
1942 if (setup_constraint_index(&ci, context) < 0)
1943 return isl_stat_error;
1945 total = isl_basic_set_total_dim(bset);
1946 for (k = 0; k < bset->n_ineq; ++k) {
1947 isl_bool redundant;
1949 redundant = constraint_index_is_redundant(&ci, bset->ineq[k]);
1950 if (redundant < 0)
1951 goto error;
1952 if (!redundant)
1953 continue;
1954 row[k] = -1;
1956 constraint_index_free(&ci);
1957 return isl_stat_ok;
1958 error:
1959 constraint_index_free(&ci);
1960 return isl_stat_error;
1963 static struct isl_basic_set *remove_shifted_constraints(
1964 struct isl_basic_set *bset, struct isl_basic_set *context)
1966 struct isl_constraint_index ci;
1967 int k;
1969 if (!bset || !context)
1970 return bset;
1972 if (context->n_ineq == 0)
1973 return bset;
1974 if (setup_constraint_index(&ci, context) < 0)
1975 return bset;
1977 for (k = 0; k < bset->n_ineq; ++k) {
1978 isl_bool redundant;
1980 redundant = constraint_index_is_redundant(&ci, bset->ineq[k]);
1981 if (redundant < 0)
1982 goto error;
1983 if (!redundant)
1984 continue;
1985 bset = isl_basic_set_cow(bset);
1986 if (!bset)
1987 goto error;
1988 isl_basic_set_drop_inequality(bset, k);
1989 --k;
1991 constraint_index_free(&ci);
1992 return bset;
1993 error:
1994 constraint_index_free(&ci);
1995 return bset;
1998 /* Remove constraints from "bmap" that are identical to constraints
1999 * in "context" or that are more relaxed (greater constant term).
2001 * We perform the test for shifted copies on the pure constraints
2002 * in remove_shifted_constraints.
2004 static __isl_give isl_basic_map *isl_basic_map_remove_shifted_constraints(
2005 __isl_take isl_basic_map *bmap, __isl_take isl_basic_map *context)
2007 isl_basic_set *bset, *bset_context;
2009 if (!bmap || !context)
2010 goto error;
2012 if (bmap->n_ineq == 0 || context->n_ineq == 0) {
2013 isl_basic_map_free(context);
2014 return bmap;
2017 context = isl_basic_map_align_divs(context, bmap);
2018 bmap = isl_basic_map_align_divs(bmap, context);
2020 bset = isl_basic_map_underlying_set(isl_basic_map_copy(bmap));
2021 bset_context = isl_basic_map_underlying_set(context);
2022 bset = remove_shifted_constraints(bset, bset_context);
2023 isl_basic_set_free(bset_context);
2025 bmap = isl_basic_map_overlying_set(bset, bmap);
2027 return bmap;
2028 error:
2029 isl_basic_map_free(bmap);
2030 isl_basic_map_free(context);
2031 return NULL;
2034 /* Does the (linear part of a) constraint "c" involve any of the "len"
2035 * "relevant" dimensions?
2037 static int is_related(isl_int *c, int len, int *relevant)
2039 int i;
2041 for (i = 0; i < len; ++i) {
2042 if (!relevant[i])
2043 continue;
2044 if (!isl_int_is_zero(c[i]))
2045 return 1;
2048 return 0;
2051 /* Drop constraints from "bset" that do not involve any of
2052 * the dimensions marked "relevant".
2054 static __isl_give isl_basic_set *drop_unrelated_constraints(
2055 __isl_take isl_basic_set *bset, int *relevant)
2057 int i, dim;
2059 dim = isl_basic_set_dim(bset, isl_dim_set);
2060 for (i = 0; i < dim; ++i)
2061 if (!relevant[i])
2062 break;
2063 if (i >= dim)
2064 return bset;
2066 for (i = bset->n_eq - 1; i >= 0; --i)
2067 if (!is_related(bset->eq[i] + 1, dim, relevant))
2068 isl_basic_set_drop_equality(bset, i);
2070 for (i = bset->n_ineq - 1; i >= 0; --i)
2071 if (!is_related(bset->ineq[i] + 1, dim, relevant))
2072 isl_basic_set_drop_inequality(bset, i);
2074 return bset;
2077 /* Update the groups in "group" based on the (linear part of a) constraint "c".
2079 * In particular, for any variable involved in the constraint,
2080 * find the actual group id from before and replace the group
2081 * of the corresponding variable by the minimal group of all
2082 * the variables involved in the constraint considered so far
2083 * (if this minimum is smaller) or replace the minimum by this group
2084 * (if the minimum is larger).
2086 * At the end, all the variables in "c" will (indirectly) point
2087 * to the minimal of the groups that they referred to originally.
2089 static void update_groups(int dim, int *group, isl_int *c)
2091 int j;
2092 int min = dim;
2094 for (j = 0; j < dim; ++j) {
2095 if (isl_int_is_zero(c[j]))
2096 continue;
2097 while (group[j] >= 0 && group[group[j]] != group[j])
2098 group[j] = group[group[j]];
2099 if (group[j] == min)
2100 continue;
2101 if (group[j] < min) {
2102 if (min >= 0 && min < dim)
2103 group[min] = group[j];
2104 min = group[j];
2105 } else
2106 group[group[j]] = min;
2110 /* Allocate an array of groups of variables, one for each variable
2111 * in "context", initialized to zero.
2113 static int *alloc_groups(__isl_keep isl_basic_set *context)
2115 isl_ctx *ctx;
2116 int dim;
2118 dim = isl_basic_set_dim(context, isl_dim_set);
2119 ctx = isl_basic_set_get_ctx(context);
2120 return isl_calloc_array(ctx, int, dim);
2123 /* Drop constraints from "context" that only involve variables that are
2124 * not related to any of the variables marked with a "-1" in "group".
2126 * We construct groups of variables that collect variables that
2127 * (indirectly) appear in some common constraint of "context".
2128 * Each group is identified by the first variable in the group,
2129 * except for the special group of variables that was already identified
2130 * in the input as -1 (or are related to those variables).
2131 * If group[i] is equal to i (or -1), then the group of i is i (or -1),
2132 * otherwise the group of i is the group of group[i].
2134 * We first initialize groups for the remaining variables.
2135 * Then we iterate over the constraints of "context" and update the
2136 * group of the variables in the constraint by the smallest group.
2137 * Finally, we resolve indirect references to groups by running over
2138 * the variables.
2140 * After computing the groups, we drop constraints that do not involve
2141 * any variables in the -1 group.
2143 static __isl_give isl_basic_set *group_and_drop_irrelevant_constraints(
2144 __isl_take isl_basic_set *context, __isl_take int *group)
2146 int dim;
2147 int i;
2148 int last;
2150 dim = isl_basic_set_dim(context, isl_dim_set);
2152 last = -1;
2153 for (i = 0; i < dim; ++i)
2154 if (group[i] >= 0)
2155 last = group[i] = i;
2156 if (last < 0) {
2157 free(group);
2158 return context;
2161 for (i = 0; i < context->n_eq; ++i)
2162 update_groups(dim, group, context->eq[i] + 1);
2163 for (i = 0; i < context->n_ineq; ++i)
2164 update_groups(dim, group, context->ineq[i] + 1);
2166 for (i = 0; i < dim; ++i)
2167 if (group[i] >= 0)
2168 group[i] = group[group[i]];
2170 for (i = 0; i < dim; ++i)
2171 group[i] = group[i] == -1;
2173 context = drop_unrelated_constraints(context, group);
2175 free(group);
2176 return context;
2179 /* Drop constraints from "context" that are irrelevant for computing
2180 * the gist of "bset".
2182 * In particular, drop constraints in variables that are not related
2183 * to any of the variables involved in the constraints of "bset"
2184 * in the sense that there is no sequence of constraints that connects them.
2186 * We first mark all variables that appear in "bset" as belonging
2187 * to a "-1" group and then continue with group_and_drop_irrelevant_constraints.
2189 static __isl_give isl_basic_set *drop_irrelevant_constraints(
2190 __isl_take isl_basic_set *context, __isl_keep isl_basic_set *bset)
2192 int *group;
2193 int dim;
2194 int i, j;
2196 if (!context || !bset)
2197 return isl_basic_set_free(context);
2199 group = alloc_groups(context);
2201 if (!group)
2202 return isl_basic_set_free(context);
2204 dim = isl_basic_set_dim(bset, isl_dim_set);
2205 for (i = 0; i < dim; ++i) {
2206 for (j = 0; j < bset->n_eq; ++j)
2207 if (!isl_int_is_zero(bset->eq[j][1 + i]))
2208 break;
2209 if (j < bset->n_eq) {
2210 group[i] = -1;
2211 continue;
2213 for (j = 0; j < bset->n_ineq; ++j)
2214 if (!isl_int_is_zero(bset->ineq[j][1 + i]))
2215 break;
2216 if (j < bset->n_ineq)
2217 group[i] = -1;
2220 return group_and_drop_irrelevant_constraints(context, group);
2223 /* Drop constraints from "context" that are irrelevant for computing
2224 * the gist of "bset", which does not involve any equalities.
2225 * Inequalities in "bset" for which the corresponding element of row
2226 * is set to -1 have already been marked for removal and should be ignored.
2228 * In particular, drop constraints in variables that are not related
2229 * to any of the variables involved in the constraints of "bset"
2230 * in the sense that there is no sequence of constraints that connects them.
2232 * We first mark all variables that appear in "bset" as belonging
2233 * to a "-1" group and then continue with group_and_drop_irrelevant_constraints.
2235 static __isl_give isl_basic_set *drop_irrelevant_constraints_marked(
2236 __isl_take isl_basic_set *context, __isl_keep isl_basic_set *bset,
2237 int *row)
2239 int *group;
2240 int dim;
2241 int i, j;
2243 if (!context || !bset)
2244 return isl_basic_set_free(context);
2246 group = alloc_groups(context);
2248 if (!group)
2249 return isl_basic_set_free(context);
2251 dim = isl_basic_set_dim(context, isl_dim_set);
2252 for (i = 0; i < dim; ++i) {
2253 for (j = 0; j < bset->n_ineq; ++j) {
2254 if (row[j] < 0)
2255 continue;
2256 if (!isl_int_is_zero(bset->ineq[j][1 + i]))
2257 break;
2259 if (j < bset->n_ineq)
2260 group[i] = -1;
2263 return group_and_drop_irrelevant_constraints(context, group);
2266 /* Do all "n" entries of "row" contain a negative value?
2268 static int all_neg(int *row, int n)
2270 int i;
2272 for (i = 0; i < n; ++i)
2273 if (row[i] >= 0)
2274 return 0;
2276 return 1;
2279 /* Update the inequalities in "bset" based on the information in "row"
2280 * and "tab".
2282 * In particular, the array "row" contains either -1, meaning that
2283 * the corresponding inequality of "bset" is redundant, or the index
2284 * of an inequality in "tab".
2286 * If the row entry is -1, then drop the inequality.
2287 * Otherwise, if the constraint is marked redundant in the tableau,
2288 * then drop the inequality.
2290 static __isl_give isl_basic_set *update_ineq(__isl_take isl_basic_set *bset,
2291 __isl_keep int *row, struct isl_tab *tab)
2293 int i;
2294 unsigned n_ineq;
2295 unsigned n_eq;
2297 if (!bset)
2298 return NULL;
2300 n_ineq = bset->n_ineq;
2301 for (i = n_ineq - 1; i >= 0; --i) {
2302 if (row[i] < 0) {
2303 if (isl_basic_set_drop_inequality(bset, i) < 0)
2304 return isl_basic_set_free(bset);
2305 continue;
2307 if (!tab)
2308 continue;
2309 n_eq = tab->n_eq;
2310 if (isl_tab_is_redundant(tab, n_eq + row[i])) {
2311 if (isl_basic_set_drop_inequality(bset, i) < 0)
2312 return isl_basic_set_free(bset);
2316 bset = isl_basic_set_finalize(bset);
2317 return bset;
2320 /* Update the inequalities in "bset" based on the information in "row"
2321 * and "tab" and free all arguments (other than "bset").
2323 static __isl_give isl_basic_set *update_ineq_free(
2324 __isl_take isl_basic_set *bset, __isl_take isl_basic_set *context,
2325 __isl_take int *row, struct isl_tab *tab)
2327 isl_basic_set_free(context);
2329 bset = update_ineq(bset, row, tab);
2331 free(row);
2332 isl_tab_free(tab);
2333 return bset;
2336 /* Remove all information from bset that is redundant in the context
2337 * of context. Both bset and context are assumed to be full-dimensional.
2339 * "row" keeps track of the constraint index of a "bset" inequality in "tab".
2340 * A value of -1 means that the inequality is obviously redundant and may
2341 * not even appear in "tab".
2343 * We first mark the inequalities of "bset"
2344 * that are obviously redundant with respect to some inequality in "context".
2345 * Then we remove those constraints from "context" that have become
2346 * irrelevant for computing the gist of "bset".
2347 * Note that this removal of constraints cannot be replaced by
2348 * a factorization because factors in "bset" may still be connected
2349 * to each other through constraints in "context".
2351 * If there are any inequalities left, we construct a tableau for
2352 * the context and then add the inequalities of "bset".
2353 * Before adding these inequalities, we freeze all constraints such that
2354 * they won't be considered redundant in terms of the constraints of "bset".
2355 * Then we detect all redundant constraints (among the
2356 * constraints that weren't frozen), first by checking for redundancy in the
2357 * the tableau and then by checking if replacing a constraint by its negation
2358 * would lead to an empty set. This last step is fairly expensive
2359 * and could be optimized by more reuse of the tableau.
2360 * Finally, we update bset according to the results.
2362 static __isl_give isl_basic_set *uset_gist_full(__isl_take isl_basic_set *bset,
2363 __isl_take isl_basic_set *context)
2365 int i, r;
2366 int *row = NULL;
2367 isl_ctx *ctx;
2368 isl_basic_set *combined = NULL;
2369 struct isl_tab *tab = NULL;
2370 unsigned context_ineq;
2371 unsigned total;
2373 if (!bset || !context)
2374 goto error;
2376 if (isl_basic_set_is_universe(bset)) {
2377 isl_basic_set_free(context);
2378 return bset;
2381 if (isl_basic_set_is_universe(context)) {
2382 isl_basic_set_free(context);
2383 return bset;
2386 ctx = isl_basic_set_get_ctx(context);
2387 row = isl_calloc_array(ctx, int, bset->n_ineq);
2388 if (!row)
2389 goto error;
2391 if (mark_shifted_constraints(bset, context, row) < 0)
2392 goto error;
2393 if (all_neg(row, bset->n_ineq))
2394 return update_ineq_free(bset, context, row, NULL);
2396 context = drop_irrelevant_constraints_marked(context, bset, row);
2397 if (!context)
2398 goto error;
2399 if (isl_basic_set_is_universe(context))
2400 return update_ineq_free(bset, context, row, NULL);
2402 context_ineq = context->n_ineq;
2403 combined = isl_basic_set_cow(isl_basic_set_copy(context));
2404 combined = isl_basic_set_extend_constraints(combined, 0, bset->n_ineq);
2405 tab = isl_tab_from_basic_set(combined, 0);
2406 for (i = 0; i < context_ineq; ++i)
2407 if (isl_tab_freeze_constraint(tab, i) < 0)
2408 goto error;
2409 if (isl_tab_extend_cons(tab, bset->n_ineq) < 0)
2410 goto error;
2411 r = context_ineq;
2412 for (i = 0; i < bset->n_ineq; ++i) {
2413 if (row[i] < 0)
2414 continue;
2415 combined = isl_basic_set_add_ineq(combined, bset->ineq[i]);
2416 if (isl_tab_add_ineq(tab, bset->ineq[i]) < 0)
2417 goto error;
2418 row[i] = r++;
2420 if (isl_tab_detect_redundant(tab) < 0)
2421 goto error;
2422 total = isl_basic_set_total_dim(bset);
2423 for (i = bset->n_ineq - 1; i >= 0; --i) {
2424 isl_basic_set *test;
2425 int is_empty;
2427 if (row[i] < 0)
2428 continue;
2429 r = row[i];
2430 if (tab->con[r].is_redundant)
2431 continue;
2432 test = isl_basic_set_dup(combined);
2433 if (isl_inequality_negate(test, r) < 0)
2434 test = isl_basic_set_free(test);
2435 test = isl_basic_set_update_from_tab(test, tab);
2436 is_empty = isl_basic_set_is_empty(test);
2437 isl_basic_set_free(test);
2438 if (is_empty < 0)
2439 goto error;
2440 if (is_empty)
2441 tab->con[r].is_redundant = 1;
2443 bset = update_ineq_free(bset, context, row, tab);
2444 if (bset) {
2445 ISL_F_SET(bset, ISL_BASIC_SET_NO_IMPLICIT);
2446 ISL_F_SET(bset, ISL_BASIC_SET_NO_REDUNDANT);
2449 isl_basic_set_free(combined);
2450 return bset;
2451 error:
2452 free(row);
2453 isl_tab_free(tab);
2454 isl_basic_set_free(combined);
2455 isl_basic_set_free(context);
2456 isl_basic_set_free(bset);
2457 return NULL;
2460 /* Remove all information from "bset" that is redundant in the context
2461 * of "context", for the case where both "bset" and "context" are
2462 * full-dimensional.
2464 static __isl_give isl_basic_set *uset_gist_uncompressed(
2465 __isl_take isl_basic_set *bset, __isl_take isl_basic_set *context)
2467 return uset_gist_full(bset, context);
2470 /* Remove all information from "bset" that is redundant in the context
2471 * of "context", for the case where the combined equalities of
2472 * "bset" and "context" allow for a compression that can be obtained
2473 * by preapplication of "T".
2474 * Preapplication of "T2" moves back to the original space.
2476 static __isl_give isl_basic_set *uset_gist_compressed(
2477 __isl_take isl_basic_set *bset, __isl_take isl_basic_set *context,
2478 __isl_take isl_mat *T, __isl_take isl_mat *T2)
2480 bset = isl_basic_set_preimage(bset, isl_mat_copy(T));
2481 context = isl_basic_set_preimage(context, T);
2483 bset = uset_gist_full(bset, context);
2484 return isl_basic_set_preimage(bset, T2);
2487 /* Remove all information from bset that is redundant in the context
2488 * of context. In particular, equalities that are linear combinations
2489 * of those in context are removed. Then the inequalities that are
2490 * redundant in the context of the equalities and inequalities of
2491 * context are removed.
2493 * First of all, we drop those constraints from "context"
2494 * that are irrelevant for computing the gist of "bset".
2495 * Alternatively, we could factorize the intersection of "context" and "bset".
2497 * We first compute the integer affine hull of the intersection,
2498 * compute the gist inside this affine hull and then add back
2499 * those equalities that are not implied by the context.
2501 * If two constraints are mutually redundant, then uset_gist_full
2502 * will remove the second of those constraints. We therefore first
2503 * sort the constraints so that constraints not involving existentially
2504 * quantified variables are given precedence over those that do.
2505 * We have to perform this sorting before the variable compression,
2506 * because that may effect the order of the variables.
2508 static __isl_give isl_basic_set *uset_gist(__isl_take isl_basic_set *bset,
2509 __isl_take isl_basic_set *context)
2511 isl_mat *eq;
2512 isl_mat *T, *T2;
2513 isl_basic_set *aff;
2514 isl_basic_set *aff_context;
2515 unsigned total;
2517 if (!bset || !context)
2518 goto error;
2520 context = drop_irrelevant_constraints(context, bset);
2522 aff = isl_basic_set_copy(bset);
2523 aff = isl_basic_set_intersect(aff, isl_basic_set_copy(context));
2524 aff = isl_basic_set_affine_hull(aff);
2525 if (!aff)
2526 goto error;
2527 if (isl_basic_set_plain_is_empty(aff)) {
2528 isl_basic_set_free(bset);
2529 isl_basic_set_free(context);
2530 return aff;
2532 bset = isl_basic_set_sort_constraints(bset);
2533 if (aff->n_eq == 0) {
2534 isl_basic_set_free(aff);
2535 return uset_gist_uncompressed(bset, context);
2537 total = isl_basic_set_total_dim(bset);
2538 eq = isl_mat_sub_alloc6(bset->ctx, aff->eq, 0, aff->n_eq, 0, 1 + total);
2539 eq = isl_mat_cow(eq);
2540 T = isl_mat_variable_compression(eq, &T2);
2541 if (T && T->n_col == 0) {
2542 isl_mat_free(T);
2543 isl_mat_free(T2);
2544 isl_basic_set_free(context);
2545 isl_basic_set_free(aff);
2546 return isl_basic_set_set_to_empty(bset);
2549 aff_context = isl_basic_set_affine_hull(isl_basic_set_copy(context));
2551 bset = uset_gist_compressed(bset, context, T, T2);
2552 bset = isl_basic_set_intersect(bset, aff);
2553 bset = isl_basic_set_reduce_using_equalities(bset, aff_context);
2555 if (bset) {
2556 ISL_F_SET(bset, ISL_BASIC_SET_NO_IMPLICIT);
2557 ISL_F_SET(bset, ISL_BASIC_SET_NO_REDUNDANT);
2560 return bset;
2561 error:
2562 isl_basic_set_free(bset);
2563 isl_basic_set_free(context);
2564 return NULL;
2567 /* Return a basic map that has the same intersection with "context" as "bmap"
2568 * and that is as "simple" as possible.
2570 * The core computation is performed on the pure constraints.
2571 * When we add back the meaning of the integer divisions, we need
2572 * to (re)introduce the div constraints. If we happen to have
2573 * discovered that some of these integer divisions are equal to
2574 * some affine combination of other variables, then these div
2575 * constraints may end up getting simplified in terms of the equalities,
2576 * resulting in extra inequalities on the other variables that
2577 * may have been removed already or that may not even have been
2578 * part of the input. We try and remove those constraints of
2579 * this form that are most obviously redundant with respect to
2580 * the context. We also remove those div constraints that are
2581 * redundant with respect to the other constraints in the result.
2583 struct isl_basic_map *isl_basic_map_gist(struct isl_basic_map *bmap,
2584 struct isl_basic_map *context)
2586 isl_basic_set *bset, *eq;
2587 isl_basic_map *eq_bmap;
2588 unsigned n_div, n_eq, n_ineq;
2590 if (!bmap || !context)
2591 goto error;
2593 if (isl_basic_map_is_universe(bmap)) {
2594 isl_basic_map_free(context);
2595 return bmap;
2597 if (isl_basic_map_plain_is_empty(context)) {
2598 isl_space *space = isl_basic_map_get_space(bmap);
2599 isl_basic_map_free(bmap);
2600 isl_basic_map_free(context);
2601 return isl_basic_map_universe(space);
2603 if (isl_basic_map_plain_is_empty(bmap)) {
2604 isl_basic_map_free(context);
2605 return bmap;
2608 bmap = isl_basic_map_remove_redundancies(bmap);
2609 context = isl_basic_map_remove_redundancies(context);
2610 if (!context)
2611 goto error;
2613 context = isl_basic_map_align_divs(context, bmap);
2614 bmap = isl_basic_map_align_divs(bmap, context);
2615 n_div = isl_basic_map_dim(bmap, isl_dim_div);
2617 bset = uset_gist(isl_basic_map_underlying_set(isl_basic_map_copy(bmap)),
2618 isl_basic_map_underlying_set(isl_basic_map_copy(context)));
2620 if (!bset || bset->n_eq == 0 || n_div == 0 ||
2621 isl_basic_set_plain_is_empty(bset)) {
2622 isl_basic_map_free(context);
2623 return isl_basic_map_overlying_set(bset, bmap);
2626 n_eq = bset->n_eq;
2627 n_ineq = bset->n_ineq;
2628 eq = isl_basic_set_copy(bset);
2629 eq = isl_basic_set_cow(eq);
2630 if (isl_basic_set_free_inequality(eq, n_ineq) < 0)
2631 eq = isl_basic_set_free(eq);
2632 if (isl_basic_set_free_equality(bset, n_eq) < 0)
2633 bset = isl_basic_set_free(bset);
2635 eq_bmap = isl_basic_map_overlying_set(eq, isl_basic_map_copy(bmap));
2636 eq_bmap = isl_basic_map_remove_shifted_constraints(eq_bmap, context);
2637 bmap = isl_basic_map_overlying_set(bset, bmap);
2638 bmap = isl_basic_map_intersect(bmap, eq_bmap);
2639 bmap = isl_basic_map_remove_redundancies(bmap);
2641 return bmap;
2642 error:
2643 isl_basic_map_free(bmap);
2644 isl_basic_map_free(context);
2645 return NULL;
2649 * Assumes context has no implicit divs.
2651 __isl_give isl_map *isl_map_gist_basic_map(__isl_take isl_map *map,
2652 __isl_take isl_basic_map *context)
2654 int i;
2656 if (!map || !context)
2657 goto error;
2659 if (isl_basic_map_plain_is_empty(context)) {
2660 isl_space *space = isl_map_get_space(map);
2661 isl_map_free(map);
2662 isl_basic_map_free(context);
2663 return isl_map_universe(space);
2666 context = isl_basic_map_remove_redundancies(context);
2667 map = isl_map_cow(map);
2668 if (!map || !context)
2669 goto error;
2670 isl_assert(map->ctx, isl_space_is_equal(map->dim, context->dim), goto error);
2671 map = isl_map_compute_divs(map);
2672 if (!map)
2673 goto error;
2674 for (i = map->n - 1; i >= 0; --i) {
2675 map->p[i] = isl_basic_map_gist(map->p[i],
2676 isl_basic_map_copy(context));
2677 if (!map->p[i])
2678 goto error;
2679 if (isl_basic_map_plain_is_empty(map->p[i])) {
2680 isl_basic_map_free(map->p[i]);
2681 if (i != map->n - 1)
2682 map->p[i] = map->p[map->n - 1];
2683 map->n--;
2686 isl_basic_map_free(context);
2687 ISL_F_CLR(map, ISL_MAP_NORMALIZED);
2688 return map;
2689 error:
2690 isl_map_free(map);
2691 isl_basic_map_free(context);
2692 return NULL;
2695 /* Return a map that has the same intersection with "context" as "map"
2696 * and that is as "simple" as possible.
2698 * If "map" is already the universe, then we cannot make it any simpler.
2699 * Similarly, if "context" is the universe, then we cannot exploit it
2700 * to simplify "map"
2701 * If "map" and "context" are identical to each other, then we can
2702 * return the corresponding universe.
2704 * If none of these cases apply, we have to work a bit harder.
2705 * During this computation, we make use of a single disjunct context,
2706 * so if the original context consists of more than one disjunct
2707 * then we need to approximate the context by a single disjunct set.
2708 * Simply taking the simple hull may drop constraints that are
2709 * only implicitly available in each disjunct. We therefore also
2710 * look for constraints among those defining "map" that are valid
2711 * for the context. These can then be used to simplify away
2712 * the corresponding constraints in "map".
2714 static __isl_give isl_map *map_gist(__isl_take isl_map *map,
2715 __isl_take isl_map *context)
2717 int equal;
2718 int is_universe;
2719 isl_basic_map *hull;
2721 is_universe = isl_map_plain_is_universe(map);
2722 if (is_universe >= 0 && !is_universe)
2723 is_universe = isl_map_plain_is_universe(context);
2724 if (is_universe < 0)
2725 goto error;
2726 if (is_universe) {
2727 isl_map_free(context);
2728 return map;
2731 equal = isl_map_plain_is_equal(map, context);
2732 if (equal < 0)
2733 goto error;
2734 if (equal) {
2735 isl_map *res = isl_map_universe(isl_map_get_space(map));
2736 isl_map_free(map);
2737 isl_map_free(context);
2738 return res;
2741 context = isl_map_compute_divs(context);
2742 if (!context)
2743 goto error;
2744 if (isl_map_n_basic_map(context) == 1) {
2745 hull = isl_map_simple_hull(context);
2746 } else {
2747 isl_ctx *ctx;
2748 isl_map_list *list;
2750 ctx = isl_map_get_ctx(map);
2751 list = isl_map_list_alloc(ctx, 2);
2752 list = isl_map_list_add(list, isl_map_copy(context));
2753 list = isl_map_list_add(list, isl_map_copy(map));
2754 hull = isl_map_unshifted_simple_hull_from_map_list(context,
2755 list);
2757 return isl_map_gist_basic_map(map, hull);
2758 error:
2759 isl_map_free(map);
2760 isl_map_free(context);
2761 return NULL;
2764 __isl_give isl_map *isl_map_gist(__isl_take isl_map *map,
2765 __isl_take isl_map *context)
2767 return isl_map_align_params_map_map_and(map, context, &map_gist);
2770 struct isl_basic_set *isl_basic_set_gist(struct isl_basic_set *bset,
2771 struct isl_basic_set *context)
2773 return (struct isl_basic_set *)isl_basic_map_gist(
2774 (struct isl_basic_map *)bset, (struct isl_basic_map *)context);
2777 __isl_give isl_set *isl_set_gist_basic_set(__isl_take isl_set *set,
2778 __isl_take isl_basic_set *context)
2780 return (struct isl_set *)isl_map_gist_basic_map((struct isl_map *)set,
2781 (struct isl_basic_map *)context);
2784 __isl_give isl_set *isl_set_gist_params_basic_set(__isl_take isl_set *set,
2785 __isl_take isl_basic_set *context)
2787 isl_space *space = isl_set_get_space(set);
2788 isl_basic_set *dom_context = isl_basic_set_universe(space);
2789 dom_context = isl_basic_set_intersect_params(dom_context, context);
2790 return isl_set_gist_basic_set(set, dom_context);
2793 __isl_give isl_set *isl_set_gist(__isl_take isl_set *set,
2794 __isl_take isl_set *context)
2796 return (struct isl_set *)isl_map_gist((struct isl_map *)set,
2797 (struct isl_map *)context);
2800 /* Compute the gist of "bmap" with respect to the constraints "context"
2801 * on the domain.
2803 __isl_give isl_basic_map *isl_basic_map_gist_domain(
2804 __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *context)
2806 isl_space *space = isl_basic_map_get_space(bmap);
2807 isl_basic_map *bmap_context = isl_basic_map_universe(space);
2809 bmap_context = isl_basic_map_intersect_domain(bmap_context, context);
2810 return isl_basic_map_gist(bmap, bmap_context);
2813 __isl_give isl_map *isl_map_gist_domain(__isl_take isl_map *map,
2814 __isl_take isl_set *context)
2816 isl_map *map_context = isl_map_universe(isl_map_get_space(map));
2817 map_context = isl_map_intersect_domain(map_context, context);
2818 return isl_map_gist(map, map_context);
2821 __isl_give isl_map *isl_map_gist_range(__isl_take isl_map *map,
2822 __isl_take isl_set *context)
2824 isl_map *map_context = isl_map_universe(isl_map_get_space(map));
2825 map_context = isl_map_intersect_range(map_context, context);
2826 return isl_map_gist(map, map_context);
2829 __isl_give isl_map *isl_map_gist_params(__isl_take isl_map *map,
2830 __isl_take isl_set *context)
2832 isl_map *map_context = isl_map_universe(isl_map_get_space(map));
2833 map_context = isl_map_intersect_params(map_context, context);
2834 return isl_map_gist(map, map_context);
2837 __isl_give isl_set *isl_set_gist_params(__isl_take isl_set *set,
2838 __isl_take isl_set *context)
2840 return isl_map_gist_params(set, context);
2843 /* Quick check to see if two basic maps are disjoint.
2844 * In particular, we reduce the equalities and inequalities of
2845 * one basic map in the context of the equalities of the other
2846 * basic map and check if we get a contradiction.
2848 isl_bool isl_basic_map_plain_is_disjoint(__isl_keep isl_basic_map *bmap1,
2849 __isl_keep isl_basic_map *bmap2)
2851 struct isl_vec *v = NULL;
2852 int *elim = NULL;
2853 unsigned total;
2854 int i;
2856 if (!bmap1 || !bmap2)
2857 return isl_bool_error;
2858 isl_assert(bmap1->ctx, isl_space_is_equal(bmap1->dim, bmap2->dim),
2859 return isl_bool_error);
2860 if (bmap1->n_div || bmap2->n_div)
2861 return isl_bool_false;
2862 if (!bmap1->n_eq && !bmap2->n_eq)
2863 return isl_bool_false;
2865 total = isl_space_dim(bmap1->dim, isl_dim_all);
2866 if (total == 0)
2867 return isl_bool_false;
2868 v = isl_vec_alloc(bmap1->ctx, 1 + total);
2869 if (!v)
2870 goto error;
2871 elim = isl_alloc_array(bmap1->ctx, int, total);
2872 if (!elim)
2873 goto error;
2874 compute_elimination_index(bmap1, elim);
2875 for (i = 0; i < bmap2->n_eq; ++i) {
2876 int reduced;
2877 reduced = reduced_using_equalities(v->block.data, bmap2->eq[i],
2878 bmap1, elim);
2879 if (reduced && !isl_int_is_zero(v->block.data[0]) &&
2880 isl_seq_first_non_zero(v->block.data + 1, total) == -1)
2881 goto disjoint;
2883 for (i = 0; i < bmap2->n_ineq; ++i) {
2884 int reduced;
2885 reduced = reduced_using_equalities(v->block.data,
2886 bmap2->ineq[i], bmap1, elim);
2887 if (reduced && isl_int_is_neg(v->block.data[0]) &&
2888 isl_seq_first_non_zero(v->block.data + 1, total) == -1)
2889 goto disjoint;
2891 compute_elimination_index(bmap2, elim);
2892 for (i = 0; i < bmap1->n_ineq; ++i) {
2893 int reduced;
2894 reduced = reduced_using_equalities(v->block.data,
2895 bmap1->ineq[i], bmap2, elim);
2896 if (reduced && isl_int_is_neg(v->block.data[0]) &&
2897 isl_seq_first_non_zero(v->block.data + 1, total) == -1)
2898 goto disjoint;
2900 isl_vec_free(v);
2901 free(elim);
2902 return isl_bool_false;
2903 disjoint:
2904 isl_vec_free(v);
2905 free(elim);
2906 return isl_bool_true;
2907 error:
2908 isl_vec_free(v);
2909 free(elim);
2910 return isl_bool_error;
2913 int isl_basic_set_plain_is_disjoint(__isl_keep isl_basic_set *bset1,
2914 __isl_keep isl_basic_set *bset2)
2916 return isl_basic_map_plain_is_disjoint((struct isl_basic_map *)bset1,
2917 (struct isl_basic_map *)bset2);
2920 /* Are "map1" and "map2" obviously disjoint?
2922 * If one of them is empty or if they live in different spaces (ignoring
2923 * parameters), then they are clearly disjoint.
2925 * If they have different parameters, then we skip any further tests.
2927 * If they are obviously equal, but not obviously empty, then we will
2928 * not be able to detect if they are disjoint.
2930 * Otherwise we check if each basic map in "map1" is obviously disjoint
2931 * from each basic map in "map2".
2933 isl_bool isl_map_plain_is_disjoint(__isl_keep isl_map *map1,
2934 __isl_keep isl_map *map2)
2936 int i, j;
2937 isl_bool disjoint;
2938 isl_bool intersect;
2939 isl_bool match;
2941 if (!map1 || !map2)
2942 return isl_bool_error;
2944 disjoint = isl_map_plain_is_empty(map1);
2945 if (disjoint < 0 || disjoint)
2946 return disjoint;
2948 disjoint = isl_map_plain_is_empty(map2);
2949 if (disjoint < 0 || disjoint)
2950 return disjoint;
2952 match = isl_space_tuple_is_equal(map1->dim, isl_dim_in,
2953 map2->dim, isl_dim_in);
2954 if (match < 0 || !match)
2955 return match < 0 ? isl_bool_error : isl_bool_true;
2957 match = isl_space_tuple_is_equal(map1->dim, isl_dim_out,
2958 map2->dim, isl_dim_out);
2959 if (match < 0 || !match)
2960 return match < 0 ? isl_bool_error : isl_bool_true;
2962 match = isl_space_match(map1->dim, isl_dim_param,
2963 map2->dim, isl_dim_param);
2964 if (match < 0 || !match)
2965 return match < 0 ? isl_bool_error : isl_bool_false;
2967 intersect = isl_map_plain_is_equal(map1, map2);
2968 if (intersect < 0 || intersect)
2969 return intersect < 0 ? isl_bool_error : isl_bool_false;
2971 for (i = 0; i < map1->n; ++i) {
2972 for (j = 0; j < map2->n; ++j) {
2973 isl_bool d = isl_basic_map_plain_is_disjoint(map1->p[i],
2974 map2->p[j]);
2975 if (d != isl_bool_true)
2976 return d;
2979 return isl_bool_true;
2982 /* Are "map1" and "map2" disjoint?
2984 * They are disjoint if they are "obviously disjoint" or if one of them
2985 * is empty. Otherwise, they are not disjoint if one of them is universal.
2986 * If none of these cases apply, we compute the intersection and see if
2987 * the result is empty.
2989 isl_bool isl_map_is_disjoint(__isl_keep isl_map *map1, __isl_keep isl_map *map2)
2991 isl_bool disjoint;
2992 isl_bool intersect;
2993 isl_map *test;
2995 disjoint = isl_map_plain_is_disjoint(map1, map2);
2996 if (disjoint < 0 || disjoint)
2997 return disjoint;
2999 disjoint = isl_map_is_empty(map1);
3000 if (disjoint < 0 || disjoint)
3001 return disjoint;
3003 disjoint = isl_map_is_empty(map2);
3004 if (disjoint < 0 || disjoint)
3005 return disjoint;
3007 intersect = isl_map_plain_is_universe(map1);
3008 if (intersect < 0 || intersect)
3009 return intersect < 0 ? isl_bool_error : isl_bool_false;
3011 intersect = isl_map_plain_is_universe(map2);
3012 if (intersect < 0 || intersect)
3013 return intersect < 0 ? isl_bool_error : isl_bool_false;
3015 test = isl_map_intersect(isl_map_copy(map1), isl_map_copy(map2));
3016 disjoint = isl_map_is_empty(test);
3017 isl_map_free(test);
3019 return disjoint;
3022 /* Are "bmap1" and "bmap2" disjoint?
3024 * They are disjoint if they are "obviously disjoint" or if one of them
3025 * is empty. Otherwise, they are not disjoint if one of them is universal.
3026 * If none of these cases apply, we compute the intersection and see if
3027 * the result is empty.
3029 isl_bool isl_basic_map_is_disjoint(__isl_keep isl_basic_map *bmap1,
3030 __isl_keep isl_basic_map *bmap2)
3032 isl_bool disjoint;
3033 isl_bool intersect;
3034 isl_basic_map *test;
3036 disjoint = isl_basic_map_plain_is_disjoint(bmap1, bmap2);
3037 if (disjoint < 0 || disjoint)
3038 return disjoint;
3040 disjoint = isl_basic_map_is_empty(bmap1);
3041 if (disjoint < 0 || disjoint)
3042 return disjoint;
3044 disjoint = isl_basic_map_is_empty(bmap2);
3045 if (disjoint < 0 || disjoint)
3046 return disjoint;
3048 intersect = isl_basic_map_is_universe(bmap1);
3049 if (intersect < 0 || intersect)
3050 return intersect < 0 ? isl_bool_error : isl_bool_false;
3052 intersect = isl_basic_map_is_universe(bmap2);
3053 if (intersect < 0 || intersect)
3054 return intersect < 0 ? isl_bool_error : isl_bool_false;
3056 test = isl_basic_map_intersect(isl_basic_map_copy(bmap1),
3057 isl_basic_map_copy(bmap2));
3058 disjoint = isl_basic_map_is_empty(test);
3059 isl_basic_map_free(test);
3061 return disjoint;
3064 /* Are "bset1" and "bset2" disjoint?
3066 isl_bool isl_basic_set_is_disjoint(__isl_keep isl_basic_set *bset1,
3067 __isl_keep isl_basic_set *bset2)
3069 return isl_basic_map_is_disjoint(bset1, bset2);
3072 isl_bool isl_set_plain_is_disjoint(__isl_keep isl_set *set1,
3073 __isl_keep isl_set *set2)
3075 return isl_map_plain_is_disjoint((struct isl_map *)set1,
3076 (struct isl_map *)set2);
3079 /* Are "set1" and "set2" disjoint?
3081 isl_bool isl_set_is_disjoint(__isl_keep isl_set *set1, __isl_keep isl_set *set2)
3083 return isl_map_is_disjoint(set1, set2);
3086 /* Check if we can combine a given div with lower bound l and upper
3087 * bound u with some other div and if so return that other div.
3088 * Otherwise return -1.
3090 * We first check that
3091 * - the bounds are opposites of each other (except for the constant
3092 * term)
3093 * - the bounds do not reference any other div
3094 * - no div is defined in terms of this div
3096 * Let m be the size of the range allowed on the div by the bounds.
3097 * That is, the bounds are of the form
3099 * e <= a <= e + m - 1
3101 * with e some expression in the other variables.
3102 * We look for another div b such that no third div is defined in terms
3103 * of this second div b and such that in any constraint that contains
3104 * a (except for the given lower and upper bound), also contains b
3105 * with a coefficient that is m times that of b.
3106 * That is, all constraints (execpt for the lower and upper bound)
3107 * are of the form
3109 * e + f (a + m b) >= 0
3111 * If so, we return b so that "a + m b" can be replaced by
3112 * a single div "c = a + m b".
3114 static int div_find_coalesce(struct isl_basic_map *bmap, int *pairs,
3115 unsigned div, unsigned l, unsigned u)
3117 int i, j;
3118 unsigned dim;
3119 int coalesce = -1;
3121 if (bmap->n_div <= 1)
3122 return -1;
3123 dim = isl_space_dim(bmap->dim, isl_dim_all);
3124 if (isl_seq_first_non_zero(bmap->ineq[l] + 1 + dim, div) != -1)
3125 return -1;
3126 if (isl_seq_first_non_zero(bmap->ineq[l] + 1 + dim + div + 1,
3127 bmap->n_div - div - 1) != -1)
3128 return -1;
3129 if (!isl_seq_is_neg(bmap->ineq[l] + 1, bmap->ineq[u] + 1,
3130 dim + bmap->n_div))
3131 return -1;
3133 for (i = 0; i < bmap->n_div; ++i) {
3134 if (isl_int_is_zero(bmap->div[i][0]))
3135 continue;
3136 if (!isl_int_is_zero(bmap->div[i][1 + 1 + dim + div]))
3137 return -1;
3140 isl_int_add(bmap->ineq[l][0], bmap->ineq[l][0], bmap->ineq[u][0]);
3141 if (isl_int_is_neg(bmap->ineq[l][0])) {
3142 isl_int_sub(bmap->ineq[l][0],
3143 bmap->ineq[l][0], bmap->ineq[u][0]);
3144 bmap = isl_basic_map_copy(bmap);
3145 bmap = isl_basic_map_set_to_empty(bmap);
3146 isl_basic_map_free(bmap);
3147 return -1;
3149 isl_int_add_ui(bmap->ineq[l][0], bmap->ineq[l][0], 1);
3150 for (i = 0; i < bmap->n_div; ++i) {
3151 if (i == div)
3152 continue;
3153 if (!pairs[i])
3154 continue;
3155 for (j = 0; j < bmap->n_div; ++j) {
3156 if (isl_int_is_zero(bmap->div[j][0]))
3157 continue;
3158 if (!isl_int_is_zero(bmap->div[j][1 + 1 + dim + i]))
3159 break;
3161 if (j < bmap->n_div)
3162 continue;
3163 for (j = 0; j < bmap->n_ineq; ++j) {
3164 int valid;
3165 if (j == l || j == u)
3166 continue;
3167 if (isl_int_is_zero(bmap->ineq[j][1 + dim + div]))
3168 continue;
3169 if (isl_int_is_zero(bmap->ineq[j][1 + dim + i]))
3170 break;
3171 isl_int_mul(bmap->ineq[j][1 + dim + div],
3172 bmap->ineq[j][1 + dim + div],
3173 bmap->ineq[l][0]);
3174 valid = isl_int_eq(bmap->ineq[j][1 + dim + div],
3175 bmap->ineq[j][1 + dim + i]);
3176 isl_int_divexact(bmap->ineq[j][1 + dim + div],
3177 bmap->ineq[j][1 + dim + div],
3178 bmap->ineq[l][0]);
3179 if (!valid)
3180 break;
3182 if (j < bmap->n_ineq)
3183 continue;
3184 coalesce = i;
3185 break;
3187 isl_int_sub_ui(bmap->ineq[l][0], bmap->ineq[l][0], 1);
3188 isl_int_sub(bmap->ineq[l][0], bmap->ineq[l][0], bmap->ineq[u][0]);
3189 return coalesce;
3192 /* Given a lower and an upper bound on div i, construct an inequality
3193 * that when nonnegative ensures that this pair of bounds always allows
3194 * for an integer value of the given div.
3195 * The lower bound is inequality l, while the upper bound is inequality u.
3196 * The constructed inequality is stored in ineq.
3197 * g, fl, fu are temporary scalars.
3199 * Let the upper bound be
3201 * -n_u a + e_u >= 0
3203 * and the lower bound
3205 * n_l a + e_l >= 0
3207 * Let n_u = f_u g and n_l = f_l g, with g = gcd(n_u, n_l).
3208 * We have
3210 * - f_u e_l <= f_u f_l g a <= f_l e_u
3212 * Since all variables are integer valued, this is equivalent to
3214 * - f_u e_l - (f_u - 1) <= f_u f_l g a <= f_l e_u + (f_l - 1)
3216 * If this interval is at least f_u f_l g, then it contains at least
3217 * one integer value for a.
3218 * That is, the test constraint is
3220 * f_l e_u + f_u e_l + f_l - 1 + f_u - 1 + 1 >= f_u f_l g
3222 static void construct_test_ineq(struct isl_basic_map *bmap, int i,
3223 int l, int u, isl_int *ineq, isl_int g, isl_int fl, isl_int fu)
3225 unsigned dim;
3226 dim = isl_space_dim(bmap->dim, isl_dim_all);
3228 isl_int_gcd(g, bmap->ineq[l][1 + dim + i], bmap->ineq[u][1 + dim + i]);
3229 isl_int_divexact(fl, bmap->ineq[l][1 + dim + i], g);
3230 isl_int_divexact(fu, bmap->ineq[u][1 + dim + i], g);
3231 isl_int_neg(fu, fu);
3232 isl_seq_combine(ineq, fl, bmap->ineq[u], fu, bmap->ineq[l],
3233 1 + dim + bmap->n_div);
3234 isl_int_add(ineq[0], ineq[0], fl);
3235 isl_int_add(ineq[0], ineq[0], fu);
3236 isl_int_sub_ui(ineq[0], ineq[0], 1);
3237 isl_int_mul(g, g, fl);
3238 isl_int_mul(g, g, fu);
3239 isl_int_sub(ineq[0], ineq[0], g);
3242 /* Remove more kinds of divs that are not strictly needed.
3243 * In particular, if all pairs of lower and upper bounds on a div
3244 * are such that they allow at least one integer value of the div,
3245 * the we can eliminate the div using Fourier-Motzkin without
3246 * introducing any spurious solutions.
3248 static struct isl_basic_map *drop_more_redundant_divs(
3249 struct isl_basic_map *bmap, int *pairs, int n)
3251 struct isl_tab *tab = NULL;
3252 struct isl_vec *vec = NULL;
3253 unsigned dim;
3254 int remove = -1;
3255 isl_int g, fl, fu;
3257 isl_int_init(g);
3258 isl_int_init(fl);
3259 isl_int_init(fu);
3261 if (!bmap)
3262 goto error;
3264 dim = isl_space_dim(bmap->dim, isl_dim_all);
3265 vec = isl_vec_alloc(bmap->ctx, 1 + dim + bmap->n_div);
3266 if (!vec)
3267 goto error;
3269 tab = isl_tab_from_basic_map(bmap, 0);
3271 while (n > 0) {
3272 int i, l, u;
3273 int best = -1;
3274 enum isl_lp_result res;
3276 for (i = 0; i < bmap->n_div; ++i) {
3277 if (!pairs[i])
3278 continue;
3279 if (best >= 0 && pairs[best] <= pairs[i])
3280 continue;
3281 best = i;
3284 i = best;
3285 for (l = 0; l < bmap->n_ineq; ++l) {
3286 if (!isl_int_is_pos(bmap->ineq[l][1 + dim + i]))
3287 continue;
3288 for (u = 0; u < bmap->n_ineq; ++u) {
3289 if (!isl_int_is_neg(bmap->ineq[u][1 + dim + i]))
3290 continue;
3291 construct_test_ineq(bmap, i, l, u,
3292 vec->el, g, fl, fu);
3293 res = isl_tab_min(tab, vec->el,
3294 bmap->ctx->one, &g, NULL, 0);
3295 if (res == isl_lp_error)
3296 goto error;
3297 if (res == isl_lp_empty) {
3298 bmap = isl_basic_map_set_to_empty(bmap);
3299 break;
3301 if (res != isl_lp_ok || isl_int_is_neg(g))
3302 break;
3304 if (u < bmap->n_ineq)
3305 break;
3307 if (l == bmap->n_ineq) {
3308 remove = i;
3309 break;
3311 pairs[i] = 0;
3312 --n;
3315 isl_tab_free(tab);
3316 isl_vec_free(vec);
3318 isl_int_clear(g);
3319 isl_int_clear(fl);
3320 isl_int_clear(fu);
3322 free(pairs);
3324 if (remove < 0)
3325 return bmap;
3327 bmap = isl_basic_map_remove_dims(bmap, isl_dim_div, remove, 1);
3328 return isl_basic_map_drop_redundant_divs(bmap);
3329 error:
3330 free(pairs);
3331 isl_basic_map_free(bmap);
3332 isl_tab_free(tab);
3333 isl_vec_free(vec);
3334 isl_int_clear(g);
3335 isl_int_clear(fl);
3336 isl_int_clear(fu);
3337 return NULL;
3340 /* Given a pair of divs div1 and div2 such that, expect for the lower bound l
3341 * and the upper bound u, div1 always occurs together with div2 in the form
3342 * (div1 + m div2), where m is the constant range on the variable div1
3343 * allowed by l and u, replace the pair div1 and div2 by a single
3344 * div that is equal to div1 + m div2.
3346 * The new div will appear in the location that contains div2.
3347 * We need to modify all constraints that contain
3348 * div2 = (div - div1) / m
3349 * (If a constraint does not contain div2, it will also not contain div1.)
3350 * If the constraint also contains div1, then we know they appear
3351 * as f (div1 + m div2) and we can simply replace (div1 + m div2) by div,
3352 * i.e., the coefficient of div is f.
3354 * Otherwise, we first need to introduce div1 into the constraint.
3355 * Let the l be
3357 * div1 + f >=0
3359 * and u
3361 * -div1 + f' >= 0
3363 * A lower bound on div2
3365 * n div2 + t >= 0
3367 * can be replaced by
3369 * (n * (m div 2 + div1) + m t + n f)/g >= 0
3371 * with g = gcd(m,n).
3372 * An upper bound
3374 * -n div2 + t >= 0
3376 * can be replaced by
3378 * (-n * (m div2 + div1) + m t + n f')/g >= 0
3380 * These constraint are those that we would obtain from eliminating
3381 * div1 using Fourier-Motzkin.
3383 * After all constraints have been modified, we drop the lower and upper
3384 * bound and then drop div1.
3386 static struct isl_basic_map *coalesce_divs(struct isl_basic_map *bmap,
3387 unsigned div1, unsigned div2, unsigned l, unsigned u)
3389 isl_int a;
3390 isl_int b;
3391 isl_int m;
3392 unsigned dim, total;
3393 int i;
3395 dim = isl_space_dim(bmap->dim, isl_dim_all);
3396 total = 1 + dim + bmap->n_div;
3398 isl_int_init(a);
3399 isl_int_init(b);
3400 isl_int_init(m);
3401 isl_int_add(m, bmap->ineq[l][0], bmap->ineq[u][0]);
3402 isl_int_add_ui(m, m, 1);
3404 for (i = 0; i < bmap->n_ineq; ++i) {
3405 if (i == l || i == u)
3406 continue;
3407 if (isl_int_is_zero(bmap->ineq[i][1 + dim + div2]))
3408 continue;
3409 if (isl_int_is_zero(bmap->ineq[i][1 + dim + div1])) {
3410 isl_int_gcd(b, m, bmap->ineq[i][1 + dim + div2]);
3411 isl_int_divexact(a, m, b);
3412 isl_int_divexact(b, bmap->ineq[i][1 + dim + div2], b);
3413 if (isl_int_is_pos(b)) {
3414 isl_seq_combine(bmap->ineq[i], a, bmap->ineq[i],
3415 b, bmap->ineq[l], total);
3416 } else {
3417 isl_int_neg(b, b);
3418 isl_seq_combine(bmap->ineq[i], a, bmap->ineq[i],
3419 b, bmap->ineq[u], total);
3422 isl_int_set(bmap->ineq[i][1 + dim + div2],
3423 bmap->ineq[i][1 + dim + div1]);
3424 isl_int_set_si(bmap->ineq[i][1 + dim + div1], 0);
3427 isl_int_clear(a);
3428 isl_int_clear(b);
3429 isl_int_clear(m);
3430 if (l > u) {
3431 isl_basic_map_drop_inequality(bmap, l);
3432 isl_basic_map_drop_inequality(bmap, u);
3433 } else {
3434 isl_basic_map_drop_inequality(bmap, u);
3435 isl_basic_map_drop_inequality(bmap, l);
3437 bmap = isl_basic_map_drop_div(bmap, div1);
3438 return bmap;
3441 /* First check if we can coalesce any pair of divs and
3442 * then continue with dropping more redundant divs.
3444 * We loop over all pairs of lower and upper bounds on a div
3445 * with coefficient 1 and -1, respectively, check if there
3446 * is any other div "c" with which we can coalesce the div
3447 * and if so, perform the coalescing.
3449 static struct isl_basic_map *coalesce_or_drop_more_redundant_divs(
3450 struct isl_basic_map *bmap, int *pairs, int n)
3452 int i, l, u;
3453 unsigned dim;
3455 dim = isl_space_dim(bmap->dim, isl_dim_all);
3457 for (i = 0; i < bmap->n_div; ++i) {
3458 if (!pairs[i])
3459 continue;
3460 for (l = 0; l < bmap->n_ineq; ++l) {
3461 if (!isl_int_is_one(bmap->ineq[l][1 + dim + i]))
3462 continue;
3463 for (u = 0; u < bmap->n_ineq; ++u) {
3464 int c;
3466 if (!isl_int_is_negone(bmap->ineq[u][1+dim+i]))
3467 continue;
3468 c = div_find_coalesce(bmap, pairs, i, l, u);
3469 if (c < 0)
3470 continue;
3471 free(pairs);
3472 bmap = coalesce_divs(bmap, i, c, l, u);
3473 return isl_basic_map_drop_redundant_divs(bmap);
3478 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
3479 return bmap;
3481 return drop_more_redundant_divs(bmap, pairs, n);
3484 /* Remove divs that are not strictly needed.
3485 * In particular, if a div only occurs positively (or negatively)
3486 * in constraints, then it can simply be dropped.
3487 * Also, if a div occurs in only two constraints and if moreover
3488 * those two constraints are opposite to each other, except for the constant
3489 * term and if the sum of the constant terms is such that for any value
3490 * of the other values, there is always at least one integer value of the
3491 * div, i.e., if one plus this sum is greater than or equal to
3492 * the (absolute value) of the coefficent of the div in the constraints,
3493 * then we can also simply drop the div.
3495 * We skip divs that appear in equalities or in the definition of other divs.
3496 * Divs that appear in the definition of other divs usually occur in at least
3497 * 4 constraints, but the constraints may have been simplified.
3499 * If any divs are left after these simple checks then we move on
3500 * to more complicated cases in drop_more_redundant_divs.
3502 struct isl_basic_map *isl_basic_map_drop_redundant_divs(
3503 struct isl_basic_map *bmap)
3505 int i, j;
3506 unsigned off;
3507 int *pairs = NULL;
3508 int n = 0;
3510 if (!bmap)
3511 goto error;
3512 if (bmap->n_div == 0)
3513 return bmap;
3515 off = isl_space_dim(bmap->dim, isl_dim_all);
3516 pairs = isl_calloc_array(bmap->ctx, int, bmap->n_div);
3517 if (!pairs)
3518 goto error;
3520 for (i = 0; i < bmap->n_div; ++i) {
3521 int pos, neg;
3522 int last_pos, last_neg;
3523 int redundant;
3524 int defined;
3526 defined = !isl_int_is_zero(bmap->div[i][0]);
3527 for (j = i; j < bmap->n_div; ++j)
3528 if (!isl_int_is_zero(bmap->div[j][1 + 1 + off + i]))
3529 break;
3530 if (j < bmap->n_div)
3531 continue;
3532 for (j = 0; j < bmap->n_eq; ++j)
3533 if (!isl_int_is_zero(bmap->eq[j][1 + off + i]))
3534 break;
3535 if (j < bmap->n_eq)
3536 continue;
3537 ++n;
3538 pos = neg = 0;
3539 for (j = 0; j < bmap->n_ineq; ++j) {
3540 if (isl_int_is_pos(bmap->ineq[j][1 + off + i])) {
3541 last_pos = j;
3542 ++pos;
3544 if (isl_int_is_neg(bmap->ineq[j][1 + off + i])) {
3545 last_neg = j;
3546 ++neg;
3549 pairs[i] = pos * neg;
3550 if (pairs[i] == 0) {
3551 for (j = bmap->n_ineq - 1; j >= 0; --j)
3552 if (!isl_int_is_zero(bmap->ineq[j][1+off+i]))
3553 isl_basic_map_drop_inequality(bmap, j);
3554 bmap = isl_basic_map_drop_div(bmap, i);
3555 free(pairs);
3556 return isl_basic_map_drop_redundant_divs(bmap);
3558 if (pairs[i] != 1)
3559 continue;
3560 if (!isl_seq_is_neg(bmap->ineq[last_pos] + 1,
3561 bmap->ineq[last_neg] + 1,
3562 off + bmap->n_div))
3563 continue;
3565 isl_int_add(bmap->ineq[last_pos][0],
3566 bmap->ineq[last_pos][0], bmap->ineq[last_neg][0]);
3567 isl_int_add_ui(bmap->ineq[last_pos][0],
3568 bmap->ineq[last_pos][0], 1);
3569 redundant = isl_int_ge(bmap->ineq[last_pos][0],
3570 bmap->ineq[last_pos][1+off+i]);
3571 isl_int_sub_ui(bmap->ineq[last_pos][0],
3572 bmap->ineq[last_pos][0], 1);
3573 isl_int_sub(bmap->ineq[last_pos][0],
3574 bmap->ineq[last_pos][0], bmap->ineq[last_neg][0]);
3575 if (!redundant) {
3576 if (defined ||
3577 !ok_to_set_div_from_bound(bmap, i, last_pos)) {
3578 pairs[i] = 0;
3579 --n;
3580 continue;
3582 bmap = set_div_from_lower_bound(bmap, i, last_pos);
3583 bmap = isl_basic_map_simplify(bmap);
3584 free(pairs);
3585 return isl_basic_map_drop_redundant_divs(bmap);
3587 if (last_pos > last_neg) {
3588 isl_basic_map_drop_inequality(bmap, last_pos);
3589 isl_basic_map_drop_inequality(bmap, last_neg);
3590 } else {
3591 isl_basic_map_drop_inequality(bmap, last_neg);
3592 isl_basic_map_drop_inequality(bmap, last_pos);
3594 bmap = isl_basic_map_drop_div(bmap, i);
3595 free(pairs);
3596 return isl_basic_map_drop_redundant_divs(bmap);
3599 if (n > 0)
3600 return coalesce_or_drop_more_redundant_divs(bmap, pairs, n);
3602 free(pairs);
3603 return bmap;
3604 error:
3605 free(pairs);
3606 isl_basic_map_free(bmap);
3607 return NULL;
3610 struct isl_basic_set *isl_basic_set_drop_redundant_divs(
3611 struct isl_basic_set *bset)
3613 return (struct isl_basic_set *)
3614 isl_basic_map_drop_redundant_divs((struct isl_basic_map *)bset);
3617 struct isl_map *isl_map_drop_redundant_divs(struct isl_map *map)
3619 int i;
3621 if (!map)
3622 return NULL;
3623 for (i = 0; i < map->n; ++i) {
3624 map->p[i] = isl_basic_map_drop_redundant_divs(map->p[i]);
3625 if (!map->p[i])
3626 goto error;
3628 ISL_F_CLR(map, ISL_MAP_NORMALIZED);
3629 return map;
3630 error:
3631 isl_map_free(map);
3632 return NULL;
3635 struct isl_set *isl_set_drop_redundant_divs(struct isl_set *set)
3637 return (struct isl_set *)
3638 isl_map_drop_redundant_divs((struct isl_map *)set);
3641 /* Does "bmap" satisfy any equality that involves more than 2 variables
3642 * and/or has coefficients different from -1 and 1?
3644 static int has_multiple_var_equality(__isl_keep isl_basic_map *bmap)
3646 int i;
3647 unsigned total;
3649 total = isl_basic_map_dim(bmap, isl_dim_all);
3651 for (i = 0; i < bmap->n_eq; ++i) {
3652 int j, k;
3654 j = isl_seq_first_non_zero(bmap->eq[i] + 1, total);
3655 if (j < 0)
3656 continue;
3657 if (!isl_int_is_one(bmap->eq[i][1 + j]) &&
3658 !isl_int_is_negone(bmap->eq[i][1 + j]))
3659 return 1;
3661 j += 1;
3662 k = isl_seq_first_non_zero(bmap->eq[i] + 1 + j, total - j);
3663 if (k < 0)
3664 continue;
3665 j += k;
3666 if (!isl_int_is_one(bmap->eq[i][1 + j]) &&
3667 !isl_int_is_negone(bmap->eq[i][1 + j]))
3668 return 1;
3670 j += 1;
3671 k = isl_seq_first_non_zero(bmap->eq[i] + 1 + j, total - j);
3672 if (k >= 0)
3673 return 1;
3676 return 0;
3679 /* Remove any common factor g from the constraint coefficients in "v".
3680 * The constant term is stored in the first position and is replaced
3681 * by floor(c/g). If any common factor is removed and if this results
3682 * in a tightening of the constraint, then set *tightened.
3684 static __isl_give isl_vec *normalize_constraint(__isl_take isl_vec *v,
3685 int *tightened)
3687 isl_ctx *ctx;
3689 if (!v)
3690 return NULL;
3691 ctx = isl_vec_get_ctx(v);
3692 isl_seq_gcd(v->el + 1, v->size - 1, &ctx->normalize_gcd);
3693 if (isl_int_is_zero(ctx->normalize_gcd))
3694 return v;
3695 if (isl_int_is_one(ctx->normalize_gcd))
3696 return v;
3697 v = isl_vec_cow(v);
3698 if (!v)
3699 return NULL;
3700 if (tightened && !isl_int_is_divisible_by(v->el[0], ctx->normalize_gcd))
3701 *tightened = 1;
3702 isl_int_fdiv_q(v->el[0], v->el[0], ctx->normalize_gcd);
3703 isl_seq_scale_down(v->el + 1, v->el + 1, ctx->normalize_gcd,
3704 v->size - 1);
3705 return v;
3708 /* If "bmap" is an integer set that satisfies any equality involving
3709 * more than 2 variables and/or has coefficients different from -1 and 1,
3710 * then use variable compression to reduce the coefficients by removing
3711 * any (hidden) common factor.
3712 * In particular, apply the variable compression to each constraint,
3713 * factor out any common factor in the non-constant coefficients and
3714 * then apply the inverse of the compression.
3715 * At the end, we mark the basic map as having reduced constants.
3716 * If this flag is still set on the next invocation of this function,
3717 * then we skip the computation.
3719 * Removing a common factor may result in a tightening of some of
3720 * the constraints. If this happens, then we may end up with two
3721 * opposite inequalities that can be replaced by an equality.
3722 * We therefore call isl_basic_map_detect_inequality_pairs,
3723 * which checks for such pairs of inequalities as well as eliminate_divs_eq
3724 * and isl_basic_map_gauss if such a pair was found.
3726 __isl_give isl_basic_map *isl_basic_map_reduce_coefficients(
3727 __isl_take isl_basic_map *bmap)
3729 unsigned total;
3730 isl_ctx *ctx;
3731 isl_vec *v;
3732 isl_mat *eq, *T, *T2;
3733 int i;
3734 int tightened;
3736 if (!bmap)
3737 return NULL;
3738 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_REDUCED_COEFFICIENTS))
3739 return bmap;
3740 if (isl_basic_map_is_rational(bmap))
3741 return bmap;
3742 if (bmap->n_eq == 0)
3743 return bmap;
3744 if (!has_multiple_var_equality(bmap))
3745 return bmap;
3747 total = isl_basic_map_dim(bmap, isl_dim_all);
3748 ctx = isl_basic_map_get_ctx(bmap);
3749 v = isl_vec_alloc(ctx, 1 + total);
3750 if (!v)
3751 return isl_basic_map_free(bmap);
3753 eq = isl_mat_sub_alloc6(ctx, bmap->eq, 0, bmap->n_eq, 0, 1 + total);
3754 T = isl_mat_variable_compression(eq, &T2);
3755 if (!T || !T2)
3756 goto error;
3757 if (T->n_col == 0) {
3758 isl_mat_free(T);
3759 isl_mat_free(T2);
3760 isl_vec_free(v);
3761 return isl_basic_map_set_to_empty(bmap);
3764 tightened = 0;
3765 for (i = 0; i < bmap->n_ineq; ++i) {
3766 isl_seq_cpy(v->el, bmap->ineq[i], 1 + total);
3767 v = isl_vec_mat_product(v, isl_mat_copy(T));
3768 v = normalize_constraint(v, &tightened);
3769 v = isl_vec_mat_product(v, isl_mat_copy(T2));
3770 if (!v)
3771 goto error;
3772 isl_seq_cpy(bmap->ineq[i], v->el, 1 + total);
3775 isl_mat_free(T);
3776 isl_mat_free(T2);
3777 isl_vec_free(v);
3779 ISL_F_SET(bmap, ISL_BASIC_MAP_REDUCED_COEFFICIENTS);
3781 if (tightened) {
3782 int progress = 0;
3784 bmap = isl_basic_map_detect_inequality_pairs(bmap, &progress);
3785 if (progress) {
3786 bmap = eliminate_divs_eq(bmap, &progress);
3787 bmap = isl_basic_map_gauss(bmap, NULL);
3791 return bmap;
3792 error:
3793 isl_mat_free(T);
3794 isl_mat_free(T2);
3795 isl_vec_free(v);
3796 return isl_basic_map_free(bmap);
3799 /* Shift the integer division at position "div" of "bmap"
3800 * by "shift" times the variable at position "pos".
3801 * "pos" is as determined by isl_basic_map_offset, i.e., pos == 0
3802 * corresponds to the constant term.
3804 * That is, if the integer division has the form
3806 * floor(f(x)/d)
3808 * then replace it by
3810 * floor((f(x) + shift * d * x_pos)/d) - shift * x_pos
3812 __isl_give isl_basic_map *isl_basic_map_shift_div(
3813 __isl_take isl_basic_map *bmap, int div, int pos, isl_int shift)
3815 int i;
3816 unsigned total;
3818 if (!bmap)
3819 return NULL;
3821 total = isl_basic_map_dim(bmap, isl_dim_all);
3822 total -= isl_basic_map_dim(bmap, isl_dim_div);
3824 isl_int_addmul(bmap->div[div][1 + pos], shift, bmap->div[div][0]);
3826 for (i = 0; i < bmap->n_eq; ++i) {
3827 if (isl_int_is_zero(bmap->eq[i][1 + total + div]))
3828 continue;
3829 isl_int_submul(bmap->eq[i][pos],
3830 shift, bmap->eq[i][1 + total + div]);
3832 for (i = 0; i < bmap->n_ineq; ++i) {
3833 if (isl_int_is_zero(bmap->ineq[i][1 + total + div]))
3834 continue;
3835 isl_int_submul(bmap->ineq[i][pos],
3836 shift, bmap->ineq[i][1 + total + div]);
3838 for (i = 0; i < bmap->n_div; ++i) {
3839 if (isl_int_is_zero(bmap->div[i][0]))
3840 continue;
3841 if (isl_int_is_zero(bmap->div[i][1 + 1 + total + div]))
3842 continue;
3843 isl_int_submul(bmap->div[i][1 + pos],
3844 shift, bmap->div[i][1 + 1 + total + div]);
3847 return bmap;