?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
add isl_pw_qpolynomial_fold_dim
[isl.git] / isl_map_simplify.c
blobd4b7af9e6fe08b14c0e173d08e29014d9560afef
1 /*
2 * Copyright 2008-2009 Katholieke Universiteit Leuven
3 *
4 * Use of this software is governed by the GNU LGPLv2.1 license
5 *
6 * Written by Sven Verdoolaege, K.U.Leuven, Departement
7 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
8 */
9
10 #include "isl_equalities.h"
11 #include "isl_map.h"
12 #include "isl_map_private.h"
13 #include "isl_seq.h"
14 #include "isl_tab.h"
15
16 static void swap_equality(struct isl_basic_map *bmap, int a, int b)
17 {
18 isl_int *t = bmap->eq[a];
19 bmap->eq[a] = bmap->eq[b];
20 bmap->eq[b] = t;
21 }
22
23 static void swap_inequality(struct isl_basic_map *bmap, int a, int b)
24 {
25 if (a != b) {
26 isl_int *t = bmap->ineq[a];
27 bmap->ineq[a] = bmap->ineq[b];
28 bmap->ineq[b] = t;
29 }
30 }
31
32 static void set_swap_inequality(struct isl_basic_set *bset, int a, int b)
33 {
34 swap_inequality((struct isl_basic_map *)bset, a, b);
35 }
36
37 static void constraint_drop_vars(isl_int *c, unsigned n, unsigned rem)
38 {
39 isl_seq_cpy(c, c + n, rem);
40 isl_seq_clr(c + rem, n);
41 }
42
43 /* Drop n dimensions starting at first.
44 *
45 * In principle, this frees up some extra variables as the number
46 * of columns remains constant, but we would have to extend
47 * the div array too as the number of rows in this array is assumed
48 * to be equal to extra.
49 */
50 struct isl_basic_set *isl_basic_set_drop_dims(
51 struct isl_basic_set *bset, unsigned first, unsigned n)
52 {
53 int i;
54
55 if (!bset)
56 goto error;
57
58 isl_assert(bset->ctx, first + n <= bset->dim->n_out, goto error);
59
60 if (n == 0)
61 return bset;
62
63 bset = isl_basic_set_cow(bset);
64 if (!bset)
65 return NULL;
66
67 for (i = 0; i < bset->n_eq; ++i)
68 constraint_drop_vars(bset->eq[i]+1+bset->dim->nparam+first, n,
69 (bset->dim->n_out-first-n)+bset->extra);
70
71 for (i = 0; i < bset->n_ineq; ++i)
72 constraint_drop_vars(bset->ineq[i]+1+bset->dim->nparam+first, n,
73 (bset->dim->n_out-first-n)+bset->extra);
74
75 for (i = 0; i < bset->n_div; ++i)
76 constraint_drop_vars(bset->div[i]+1+1+bset->dim->nparam+first, n,
77 (bset->dim->n_out-first-n)+bset->extra);
78
79 bset->dim = isl_dim_drop_outputs(bset->dim, first, n);
80 if (!bset->dim)
81 goto error;
82
83 ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED);
84 bset = isl_basic_set_simplify(bset);
85 return isl_basic_set_finalize(bset);
86 error:
87 isl_basic_set_free(bset);
88 return NULL;
89 }
90
91 struct isl_set *isl_set_drop_dims(
92 struct isl_set *set, unsigned first, unsigned n)
93 {
94 int i;
95
96 if (!set)
97 goto error;
98
99 isl_assert(set->ctx, first + n <= set->dim->n_out, goto error);
100
101 if (n == 0)
102 return set;
103 set = isl_set_cow(set);
104 if (!set)
105 goto error;
106 set->dim = isl_dim_drop_outputs(set->dim, first, n);
107 if (!set->dim)
108 goto error;
109
110 for (i = 0; i < set->n; ++i) {
111 set->p[i] = isl_basic_set_drop_dims(set->p[i], first, n);
112 if (!set->p[i])
113 goto error;
114 }
115
116 ISL_F_CLR(set, ISL_SET_NORMALIZED);
117 return set;
118 error:
119 isl_set_free(set);
120 return NULL;
121 }
122
123 /* Move "n" divs starting at "first" to the end of the list of divs.
124 */
125 static struct isl_basic_map *move_divs_last(struct isl_basic_map *bmap,
126 unsigned first, unsigned n)
127 {
128 isl_int **div;
129 int i;
130
131 if (first + n == bmap->n_div)
132 return bmap;
133
134 div = isl_alloc_array(bmap->ctx, isl_int *, n);
135 if (!div)
136 goto error;
137 for (i = 0; i < n; ++i)
138 div[i] = bmap->div[first + i];
139 for (i = 0; i < bmap->n_div - first - n; ++i)
140 bmap->div[first + i] = bmap->div[first + n + i];
141 for (i = 0; i < n; ++i)
142 bmap->div[bmap->n_div - n + i] = div[i];
143 free(div);
144 return bmap;
145 error:
146 isl_basic_map_free(bmap);
147 return NULL;
148 }
149
150 /* Drop "n" dimensions of type "type" starting at "first".
151 *
152 * In principle, this frees up some extra variables as the number
153 * of columns remains constant, but we would have to extend
154 * the div array too as the number of rows in this array is assumed
155 * to be equal to extra.
156 */
157 struct isl_basic_map *isl_basic_map_drop(struct isl_basic_map *bmap,
158 enum isl_dim_type type, unsigned first, unsigned n)
159 {
160 int i;
161 unsigned dim;
162 unsigned offset;
163 unsigned left;
164
165 if (!bmap)
166 goto error;
167
168 dim = isl_basic_map_dim(bmap, type);
169 isl_assert(bmap->ctx, first + n <= dim, goto error);
170
171 if (n == 0)
172 return bmap;
173
174 bmap = isl_basic_map_cow(bmap);
175 if (!bmap)
176 return NULL;
177
178 offset = isl_basic_map_offset(bmap, type) + first;
179 left = isl_basic_map_total_dim(bmap) - (offset - 1) - n;
180 for (i = 0; i < bmap->n_eq; ++i)
181 constraint_drop_vars(bmap->eq[i]+offset, n, left);
182
183 for (i = 0; i < bmap->n_ineq; ++i)
184 constraint_drop_vars(bmap->ineq[i]+offset, n, left);
185
186 for (i = 0; i < bmap->n_div; ++i)
187 constraint_drop_vars(bmap->div[i]+1+offset, n, left);
188
189 if (type == isl_dim_div) {
190 bmap = move_divs_last(bmap, first, n);
191 if (!bmap)
192 goto error;
193 isl_basic_map_free_div(bmap, n);
194 } else
195 bmap->dim = isl_dim_drop(bmap->dim, type, first, n);
196 if (!bmap->dim)
197 goto error;
198
199 ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
200 bmap = isl_basic_map_simplify(bmap);
201 return isl_basic_map_finalize(bmap);
202 error:
203 isl_basic_map_free(bmap);
204 return NULL;
205 }
206
207 struct isl_basic_map *isl_basic_map_drop_inputs(
208 struct isl_basic_map *bmap, unsigned first, unsigned n)
209 {
210 return isl_basic_map_drop(bmap, isl_dim_in, first, n);
211 }
212
213 struct isl_map *isl_map_drop(struct isl_map *map,
214 enum isl_dim_type type, unsigned first, unsigned n)
215 {
216 int i;
217
218 if (!map)
219 goto error;
220
221 isl_assert(map->ctx, first + n <= isl_map_dim(map, type), goto error);
222
223 if (n == 0)
224 return map;
225 map = isl_map_cow(map);
226 if (!map)
227 goto error;
228 map->dim = isl_dim_drop(map->dim, type, first, n);
229 if (!map->dim)
230 goto error;
231
232 for (i = 0; i < map->n; ++i) {
233 map->p[i] = isl_basic_map_drop(map->p[i], type, first, n);
234 if (!map->p[i])
235 goto error;
236 }
237 ISL_F_CLR(map, ISL_MAP_NORMALIZED);
238
239 return map;
240 error:
241 isl_map_free(map);
242 return NULL;
243 }
244
245 struct isl_set *isl_set_drop(struct isl_set *set,
246 enum isl_dim_type type, unsigned first, unsigned n)
247 {
248 return (isl_set *)isl_map_drop((isl_map *)set, type, first, n);
249 }
250
251 struct isl_map *isl_map_drop_inputs(
252 struct isl_map *map, unsigned first, unsigned n)
253 {
254 return isl_map_drop(map, isl_dim_in, first, n);
255 }
256
257 /*
258 * We don't cow, as the div is assumed to be redundant.
259 */
260 static struct isl_basic_map *isl_basic_map_drop_div(
261 struct isl_basic_map *bmap, unsigned div)
262 {
263 int i;
264 unsigned pos;
265
266 if (!bmap)
267 goto error;
268
269 pos = 1 + isl_dim_total(bmap->dim) + div;
270
271 isl_assert(bmap->ctx, div < bmap->n_div, goto error);
272
273 for (i = 0; i < bmap->n_eq; ++i)
274 constraint_drop_vars(bmap->eq[i]+pos, 1, bmap->extra-div-1);
275
276 for (i = 0; i < bmap->n_ineq; ++i) {
277 if (!isl_int_is_zero(bmap->ineq[i][pos])) {
278 isl_basic_map_drop_inequality(bmap, i);
279 --i;
280 continue;
281 }
282 constraint_drop_vars(bmap->ineq[i]+pos, 1, bmap->extra-div-1);
283 }
284
285 for (i = 0; i < bmap->n_div; ++i)
286 constraint_drop_vars(bmap->div[i]+1+pos, 1, bmap->extra-div-1);
287
288 if (div != bmap->n_div - 1) {
289 int j;
290 isl_int *t = bmap->div[div];
291
292 for (j = div; j < bmap->n_div - 1; ++j)
293 bmap->div[j] = bmap->div[j+1];
294
295 bmap->div[bmap->n_div - 1] = t;
296 }
297 ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
298 isl_basic_map_free_div(bmap, 1);
299
300 return bmap;
301 error:
302 isl_basic_map_free(bmap);
303 return NULL;
304 }
305
306 struct isl_basic_map *isl_basic_map_normalize_constraints(
307 struct isl_basic_map *bmap)
308 {
309 int i;
310 isl_int gcd;
311 unsigned total = isl_basic_map_total_dim(bmap);
312
313 isl_int_init(gcd);
314 for (i = bmap->n_eq - 1; i >= 0; --i) {
315 isl_seq_gcd(bmap->eq[i]+1, total, &gcd);
316 if (isl_int_is_zero(gcd)) {
317 if (!isl_int_is_zero(bmap->eq[i][0])) {
318 bmap = isl_basic_map_set_to_empty(bmap);
319 break;
320 }
321 isl_basic_map_drop_equality(bmap, i);
322 continue;
323 }
324 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
325 isl_int_gcd(gcd, gcd, bmap->eq[i][0]);
326 if (isl_int_is_one(gcd))
327 continue;
328 if (!isl_int_is_divisible_by(bmap->eq[i][0], gcd)) {
329 bmap = isl_basic_map_set_to_empty(bmap);
330 break;
331 }
332 isl_seq_scale_down(bmap->eq[i], bmap->eq[i], gcd, 1+total);
333 }
334
335 for (i = bmap->n_ineq - 1; i >= 0; --i) {
336 isl_seq_gcd(bmap->ineq[i]+1, total, &gcd);
337 if (isl_int_is_zero(gcd)) {
338 if (isl_int_is_neg(bmap->ineq[i][0])) {
339 bmap = isl_basic_map_set_to_empty(bmap);
340 break;
341 }
342 isl_basic_map_drop_inequality(bmap, i);
343 continue;
344 }
345 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
346 isl_int_gcd(gcd, gcd, bmap->ineq[i][0]);
347 if (isl_int_is_one(gcd))
348 continue;
349 isl_int_fdiv_q(bmap->ineq[i][0], bmap->ineq[i][0], gcd);
350 isl_seq_scale_down(bmap->ineq[i]+1, bmap->ineq[i]+1, gcd, total);
351 }
352 isl_int_clear(gcd);
353
354 return bmap;
355 }
356
357 struct isl_basic_set *isl_basic_set_normalize_constraints(
358 struct isl_basic_set *bset)
359 {
360 return (struct isl_basic_set *)isl_basic_map_normalize_constraints(
361 (struct isl_basic_map *)bset);
362 }
363
364 /* Assumes divs have been ordered if keep_divs is set.
365 */
366 static void eliminate_var_using_equality(struct isl_basic_map *bmap,
367 unsigned pos, isl_int *eq, int keep_divs, int *progress)
368 {
369 unsigned total;
370 int k;
371 int last_div;
372
373 total = isl_basic_map_total_dim(bmap);
374 last_div = isl_seq_last_non_zero(eq + 1 + isl_dim_total(bmap->dim),
375 bmap->n_div);
376 for (k = 0; k < bmap->n_eq; ++k) {
377 if (bmap->eq[k] == eq)
378 continue;
379 if (isl_int_is_zero(bmap->eq[k][1+pos]))
380 continue;
381 if (progress)
382 *progress = 1;
383 isl_seq_elim(bmap->eq[k], eq, 1+pos, 1+total, NULL);
384 }
385
386 for (k = 0; k < bmap->n_ineq; ++k) {
387 if (isl_int_is_zero(bmap->ineq[k][1+pos]))
388 continue;
389 if (progress)
390 *progress = 1;
391 isl_seq_elim(bmap->ineq[k], eq, 1+pos, 1+total, NULL);
392 ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
393 }
394
395 for (k = 0; k < bmap->n_div; ++k) {
396 if (isl_int_is_zero(bmap->div[k][0]))
397 continue;
398 if (isl_int_is_zero(bmap->div[k][1+1+pos]))
399 continue;
400 if (progress)
401 *progress = 1;
402 /* We need to be careful about circular definitions,
403 * so for now we just remove the definition of div k
404 * if the equality contains any divs.
405 * If keep_divs is set, then the divs have been ordered
406 * and we can keep the definition as long as the result
407 * is still ordered.
408 */
409 if (last_div == -1 || (keep_divs && last_div < k))
410 isl_seq_elim(bmap->div[k]+1, eq,
411 1+pos, 1+total, &bmap->div[k][0]);
412 else
413 isl_seq_clr(bmap->div[k], 1 + total);
414 ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
415 }
416 }
417
418 /* Assumes divs have been ordered if keep_divs is set.
419 */
420 static void eliminate_div(struct isl_basic_map *bmap, isl_int *eq,
421 unsigned div, int keep_divs)
422 {
423 unsigned pos = isl_dim_total(bmap->dim) + div;
424
425 eliminate_var_using_equality(bmap, pos, eq, keep_divs, NULL);
426
427 isl_basic_map_drop_div(bmap, div);
428 }
429
430 /* Elimininate divs based on equalities
431 */
432 static struct isl_basic_map *eliminate_divs_eq(
433 struct isl_basic_map *bmap, int *progress)
434 {
435 int d;
436 int i;
437 int modified = 0;
438 unsigned off;
439
440 bmap = isl_basic_map_order_divs(bmap);
441
442 if (!bmap)
443 return NULL;
444
445 off = 1 + isl_dim_total(bmap->dim);
446
447 for (d = bmap->n_div - 1; d >= 0 ; --d) {
448 for (i = 0; i < bmap->n_eq; ++i) {
449 if (!isl_int_is_one(bmap->eq[i][off + d]) &&
450 !isl_int_is_negone(bmap->eq[i][off + d]))
451 continue;
452 modified = 1;
453 *progress = 1;
454 eliminate_div(bmap, bmap->eq[i], d, 1);
455 isl_basic_map_drop_equality(bmap, i);
456 break;
457 }
458 }
459 if (modified)
460 return eliminate_divs_eq(bmap, progress);
461 return bmap;
462 }
463
464 /* Elimininate divs based on inequalities
465 */
466 static struct isl_basic_map *eliminate_divs_ineq(
467 struct isl_basic_map *bmap, int *progress)
468 {
469 int d;
470 int i;
471 unsigned off;
472 struct isl_ctx *ctx;
473
474 if (!bmap)
475 return NULL;
476
477 ctx = bmap->ctx;
478 off = 1 + isl_dim_total(bmap->dim);
479
480 for (d = bmap->n_div - 1; d >= 0 ; --d) {
481 for (i = 0; i < bmap->n_eq; ++i)
482 if (!isl_int_is_zero(bmap->eq[i][off + d]))
483 break;
484 if (i < bmap->n_eq)
485 continue;
486 for (i = 0; i < bmap->n_ineq; ++i)
487 if (isl_int_abs_gt(bmap->ineq[i][off + d], ctx->one))
488 break;
489 if (i < bmap->n_ineq)
490 continue;
491 *progress = 1;
492 bmap = isl_basic_map_eliminate_vars(bmap, (off-1)+d, 1);
493 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
494 break;
495 bmap = isl_basic_map_drop_div(bmap, d);
496 if (!bmap)
497 break;
498 }
499 return bmap;
500 }
501
502 struct isl_basic_map *isl_basic_map_gauss(
503 struct isl_basic_map *bmap, int *progress)
504 {
505 int k;
506 int done;
507 int last_var;
508 unsigned total_var;
509 unsigned total;
510
511 bmap = isl_basic_map_order_divs(bmap);
512
513 if (!bmap)
514 return NULL;
515
516 total = isl_basic_map_total_dim(bmap);
517 total_var = total - bmap->n_div;
518
519 last_var = total - 1;
520 for (done = 0; done < bmap->n_eq; ++done) {
521 for (; last_var >= 0; --last_var) {
522 for (k = done; k < bmap->n_eq; ++k)
523 if (!isl_int_is_zero(bmap->eq[k][1+last_var]))
524 break;
525 if (k < bmap->n_eq)
526 break;
527 }
528 if (last_var < 0)
529 break;
530 if (k != done)
531 swap_equality(bmap, k, done);
532 if (isl_int_is_neg(bmap->eq[done][1+last_var]))
533 isl_seq_neg(bmap->eq[done], bmap->eq[done], 1+total);
534
535 eliminate_var_using_equality(bmap, last_var, bmap->eq[done], 1,
536 progress);
537
538 if (last_var >= total_var &&
539 isl_int_is_zero(bmap->div[last_var - total_var][0])) {
540 unsigned div = last_var - total_var;
541 isl_seq_neg(bmap->div[div]+1, bmap->eq[done], 1+total);
542 isl_int_set_si(bmap->div[div][1+1+last_var], 0);
543 isl_int_set(bmap->div[div][0],
544 bmap->eq[done][1+last_var]);
545 ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
546 }
547 }
548 if (done == bmap->n_eq)
549 return bmap;
550 for (k = done; k < bmap->n_eq; ++k) {
551 if (isl_int_is_zero(bmap->eq[k][0]))
552 continue;
553 return isl_basic_map_set_to_empty(bmap);
554 }
555 isl_basic_map_free_equality(bmap, bmap->n_eq-done);
556 return bmap;
557 }
558
559 struct isl_basic_set *isl_basic_set_gauss(
560 struct isl_basic_set *bset, int *progress)
561 {
562 return (struct isl_basic_set*)isl_basic_map_gauss(
563 (struct isl_basic_map *)bset, progress);
564 }
565
566
567 static unsigned int round_up(unsigned int v)
568 {
569 int old_v = v;
570
571 while (v) {
572 old_v = v;
573 v ^= v & -v;
574 }
575 return old_v << 1;
576 }
577
578 static int hash_index(isl_int ***index, unsigned int size, int bits,
579 struct isl_basic_map *bmap, int k)
580 {
581 int h;
582 unsigned total = isl_basic_map_total_dim(bmap);
583 uint32_t hash = isl_seq_get_hash_bits(bmap->ineq[k]+1, total, bits);
584 for (h = hash; index[h]; h = (h+1) % size)
585 if (&bmap->ineq[k] != index[h] &&
586 isl_seq_eq(bmap->ineq[k]+1, index[h][0]+1, total))
587 break;
588 return h;
589 }
590
591 static int set_hash_index(isl_int ***index, unsigned int size, int bits,
592 struct isl_basic_set *bset, int k)
593 {
594 return hash_index(index, size, bits, (struct isl_basic_map *)bset, k);
595 }
596
597 /* If we can eliminate more than one div, then we need to make
598 * sure we do it from last div to first div, in order not to
599 * change the position of the other divs that still need to
600 * be removed.
601 */
602 static struct isl_basic_map *remove_duplicate_divs(
603 struct isl_basic_map *bmap, int *progress)
604 {
605 unsigned int size;
606 int *index;
607 int *elim_for;
608 int k, l, h;
609 int bits;
610 struct isl_blk eq;
611 unsigned total_var = isl_dim_total(bmap->dim);
612 unsigned total = total_var + bmap->n_div;
613 struct isl_ctx *ctx;
614
615 if (bmap->n_div <= 1)
616 return bmap;
617
618 ctx = bmap->ctx;
619 for (k = bmap->n_div - 1; k >= 0; --k)
620 if (!isl_int_is_zero(bmap->div[k][0]))
621 break;
622 if (k <= 0)
623 return bmap;
624
625 elim_for = isl_calloc_array(ctx, int, bmap->n_div);
626 size = round_up(4 * bmap->n_div / 3 - 1);
627 bits = ffs(size) - 1;
628 index = isl_calloc_array(ctx, int, size);
629 if (!index)
630 return bmap;
631 eq = isl_blk_alloc(ctx, 1+total);
632 if (isl_blk_is_error(eq))
633 goto out;
634
635 isl_seq_clr(eq.data, 1+total);
636 index[isl_seq_get_hash_bits(bmap->div[k], 2+total, bits)] = k + 1;
637 for (--k; k >= 0; --k) {
638 uint32_t hash;
639
640 if (isl_int_is_zero(bmap->div[k][0]))
641 continue;
642
643 hash = isl_seq_get_hash_bits(bmap->div[k], 2+total, bits);
644 for (h = hash; index[h]; h = (h+1) % size)
645 if (isl_seq_eq(bmap->div[k],
646 bmap->div[index[h]-1], 2+total))
647 break;
648 if (index[h]) {
649 *progress = 1;
650 l = index[h] - 1;
651 elim_for[l] = k + 1;
652 }
653 index[h] = k+1;
654 }
655 for (l = bmap->n_div - 1; l >= 0; --l) {
656 if (!elim_for[l])
657 continue;
658 k = elim_for[l] - 1;
659 isl_int_set_si(eq.data[1+total_var+k], -1);
660 isl_int_set_si(eq.data[1+total_var+l], 1);
661 eliminate_div(bmap, eq.data, l, 0);
662 isl_int_set_si(eq.data[1+total_var+k], 0);
663 isl_int_set_si(eq.data[1+total_var+l], 0);
664 }
665
666 isl_blk_free(ctx, eq);
667 out:
668 free(index);
669 free(elim_for);
670 return bmap;
671 }
672
673 static int n_pure_div_eq(struct isl_basic_map *bmap)
674 {
675 int i, j;
676 unsigned total;
677
678 total = isl_dim_total(bmap->dim);
679 for (i = 0, j = bmap->n_div-1; i < bmap->n_eq; ++i) {
680 while (j >= 0 && isl_int_is_zero(bmap->eq[i][1 + total + j]))
681 --j;
682 if (j < 0)
683 break;
684 if (isl_seq_first_non_zero(bmap->eq[i] + 1 + total, j) != -1)
685 return 0;
686 }
687 return i;
688 }
689
690 /* Normalize divs that appear in equalities.
691 *
692 * In particular, we assume that bmap contains some equalities
693 * of the form
694 *
695 * a x = m * e_i
696 *
697 * and we want to replace the set of e_i by a minimal set and
698 * such that the new e_i have a canonical representation in terms
699 * of the vector x.
700 * If any of the equalities involves more than one divs, then
701 * we currently simply bail out.
702 *
703 * Let us first additionally assume that all equalities involve
704 * a div. The equalities then express modulo constraints on the
705 * remaining variables and we can use "parameter compression"
706 * to find a minimal set of constraints. The result is a transformation
707 *
708 * x = T(x') = x_0 + G x'
709 *
710 * with G a lower-triangular matrix with all elements below the diagonal
711 * non-negative and smaller than the diagonal element on the same row.
712 * We first normalize x_0 by making the same property hold in the affine
713 * T matrix.
714 * The rows i of G with a 1 on the diagonal do not impose any modulo
715 * constraint and simply express x_i = x'_i.
716 * For each of the remaining rows i, we introduce a div and a corresponding
717 * equality. In particular
718 *
719 * g_ii e_j = x_i - g_i(x')
720 *
721 * where each x'_k is replaced either by x_k (if g_kk = 1) or the
722 * corresponding div (if g_kk != 1).
723 *
724 * If there are any equalities not involving any div, then we
725 * first apply a variable compression on the variables x:
726 *
727 * x = C x'' x'' = C_2 x
728 *
729 * and perform the above parameter compression on A C instead of on A.
730 * The resulting compression is then of the form
731 *
732 * x'' = T(x') = x_0 + G x'
733 *
734 * and in constructing the new divs and the corresponding equalities,
735 * we have to replace each x'', i.e., the x'_k with (g_kk = 1),
736 * by the corresponding row from C_2.
737 */
738 static struct isl_basic_map *normalize_divs(
739 struct isl_basic_map *bmap, int *progress)
740 {
741 int i, j, k;
742 int total;
743 int div_eq;
744 struct isl_mat *B;
745 struct isl_vec *d;
746 struct isl_mat *T = NULL;
747 struct isl_mat *C = NULL;
748 struct isl_mat *C2 = NULL;
749 isl_int v;
750 int *pos;
751 int dropped, needed;
752
753 if (!bmap)
754 return NULL;
755
756 if (bmap->n_div == 0)
757 return bmap;
758
759 if (bmap->n_eq == 0)
760 return bmap;
761
762 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_NORMALIZED_DIVS))
763 return bmap;
764
765 total = isl_dim_total(bmap->dim);
766 div_eq = n_pure_div_eq(bmap);
767 if (div_eq == 0)
768 return bmap;
769
770 if (div_eq < bmap->n_eq) {
771 B = isl_mat_sub_alloc(bmap->ctx, bmap->eq, div_eq,
772 bmap->n_eq - div_eq, 0, 1 + total);
773 C = isl_mat_variable_compression(B, &C2);
774 if (!C || !C2)
775 goto error;
776 if (C->n_col == 0) {
777 bmap = isl_basic_map_set_to_empty(bmap);
778 isl_mat_free(C);
779 isl_mat_free(C2);
780 goto done;
781 }
782 }
783
784 d = isl_vec_alloc(bmap->ctx, div_eq);
785 if (!d)
786 goto error;
787 for (i = 0, j = bmap->n_div-1; i < div_eq; ++i) {
788 while (j >= 0 && isl_int_is_zero(bmap->eq[i][1 + total + j]))
789 --j;
790 isl_int_set(d->block.data[i], bmap->eq[i][1 + total + j]);
791 }
792 B = isl_mat_sub_alloc(bmap->ctx, bmap->eq, 0, div_eq, 0, 1 + total);
793
794 if (C) {
795 B = isl_mat_product(B, C);
796 C = NULL;
797 }
798
799 T = isl_mat_parameter_compression(B, d);
800 if (!T)
801 goto error;
802 if (T->n_col == 0) {
803 bmap = isl_basic_map_set_to_empty(bmap);
804 isl_mat_free(C2);
805 isl_mat_free(T);
806 goto done;
807 }
808 isl_int_init(v);
809 for (i = 0; i < T->n_row - 1; ++i) {
810 isl_int_fdiv_q(v, T->row[1 + i][0], T->row[1 + i][1 + i]);
811 if (isl_int_is_zero(v))
812 continue;
813 isl_mat_col_submul(T, 0, v, 1 + i);
814 }
815 isl_int_clear(v);
816 pos = isl_alloc_array(bmap->ctx, int, T->n_row);
817 /* We have to be careful because dropping equalities may reorder them */
818 dropped = 0;
819 for (j = bmap->n_div - 1; j >= 0; --j) {
820 for (i = 0; i < bmap->n_eq; ++i)
821 if (!isl_int_is_zero(bmap->eq[i][1 + total + j]))
822 break;
823 if (i < bmap->n_eq) {
824 bmap = isl_basic_map_drop_div(bmap, j);
825 isl_basic_map_drop_equality(bmap, i);
826 ++dropped;
827 }
828 }
829 pos[0] = 0;
830 needed = 0;
831 for (i = 1; i < T->n_row; ++i) {
832 if (isl_int_is_one(T->row[i][i]))
833 pos[i] = i;
834 else
835 needed++;
836 }
837 if (needed > dropped) {
838 bmap = isl_basic_map_extend_dim(bmap, isl_dim_copy(bmap->dim),
839 needed, needed, 0);
840 if (!bmap)
841 goto error;
842 }
843 for (i = 1; i < T->n_row; ++i) {
844 if (isl_int_is_one(T->row[i][i]))
845 continue;
846 k = isl_basic_map_alloc_div(bmap);
847 pos[i] = 1 + total + k;
848 isl_seq_clr(bmap->div[k] + 1, 1 + total + bmap->n_div);
849 isl_int_set(bmap->div[k][0], T->row[i][i]);
850 if (C2)
851 isl_seq_cpy(bmap->div[k] + 1, C2->row[i], 1 + total);
852 else
853 isl_int_set_si(bmap->div[k][1 + i], 1);
854 for (j = 0; j < i; ++j) {
855 if (isl_int_is_zero(T->row[i][j]))
856 continue;
857 if (pos[j] < T->n_row && C2)
858 isl_seq_submul(bmap->div[k] + 1, T->row[i][j],
859 C2->row[pos[j]], 1 + total);
860 else
861 isl_int_neg(bmap->div[k][1 + pos[j]],
862 T->row[i][j]);
863 }
864 j = isl_basic_map_alloc_equality(bmap);
865 isl_seq_neg(bmap->eq[j], bmap->div[k]+1, 1+total+bmap->n_div);
866 isl_int_set(bmap->eq[j][pos[i]], bmap->div[k][0]);
867 }
868 free(pos);
869 isl_mat_free(C2);
870 isl_mat_free(T);
871
872 if (progress)
873 *progress = 1;
874 done:
875 ISL_F_SET(bmap, ISL_BASIC_MAP_NORMALIZED_DIVS);
876
877 return bmap;
878 error:
879 isl_mat_free(C);
880 isl_mat_free(C2);
881 isl_mat_free(T);
882 return bmap;
883 }
884
885 static struct isl_basic_map *set_div_from_lower_bound(
886 struct isl_basic_map *bmap, int div, int ineq)
887 {
888 unsigned total = 1 + isl_dim_total(bmap->dim);
889
890 isl_seq_neg(bmap->div[div] + 1, bmap->ineq[ineq], total + bmap->n_div);
891 isl_int_set(bmap->div[div][0], bmap->ineq[ineq][total + div]);
892 isl_int_add(bmap->div[div][1], bmap->div[div][1], bmap->div[div][0]);
893 isl_int_sub_ui(bmap->div[div][1], bmap->div[div][1], 1);
894 isl_int_set_si(bmap->div[div][1 + total + div], 0);
895
896 return bmap;
897 }
898
899 /* Check whether it is ok to define a div based on an inequality.
900 * To avoid the introduction of circular definitions of divs, we
901 * do not allow such a definition if the resulting expression would refer to
902 * any other undefined divs or if any known div is defined in
903 * terms of the unknown div.
904 */
905 static int ok_to_set_div_from_bound(struct isl_basic_map *bmap,
906 int div, int ineq)
907 {
908 int j;
909 unsigned total = 1 + isl_dim_total(bmap->dim);
910
911 /* Not defined in terms of unknown divs */
912 for (j = 0; j < bmap->n_div; ++j) {
913 if (div == j)
914 continue;
915 if (isl_int_is_zero(bmap->ineq[ineq][total + j]))
916 continue;
917 if (isl_int_is_zero(bmap->div[j][0]))
918 return 0;
919 }
920
921 /* No other div defined in terms of this one => avoid loops */
922 for (j = 0; j < bmap->n_div; ++j) {
923 if (div == j)
924 continue;
925 if (isl_int_is_zero(bmap->div[j][0]))
926 continue;
927 if (!isl_int_is_zero(bmap->div[j][1 + total + div]))
928 return 0;
929 }
930
931 return 1;
932 }
933
934 /* Given two constraints "k" and "l" that are opposite to each other,
935 * except for the constant term, check if we can use them
936 * to obtain an expression for one of the hitherto unknown divs.
937 * "sum" is the sum of the constant terms of the constraints.
938 * If this sum is strictly smaller than the coefficient of one
939 * of the divs, then this pair can be used define the div.
940 * To avoid the introduction of circular definitions of divs, we
941 * do not use the pair if the resulting expression would refer to
942 * any other undefined divs or if any known div is defined in
943 * terms of the unknown div.
944 */
945 static struct isl_basic_map *check_for_div_constraints(
946 struct isl_basic_map *bmap, int k, int l, isl_int sum, int *progress)
947 {
948 int i;
949 unsigned total = 1 + isl_dim_total(bmap->dim);
950
951 for (i = 0; i < bmap->n_div; ++i) {
952 if (!isl_int_is_zero(bmap->div[i][0]))
953 continue;
954 if (isl_int_is_zero(bmap->ineq[k][total + i]))
955 continue;
956 if (isl_int_abs_ge(sum, bmap->ineq[k][total + i]))
957 continue;
958 if (!ok_to_set_div_from_bound(bmap, i, k))
959 break;
960 if (isl_int_is_pos(bmap->ineq[k][total + i]))
961 bmap = set_div_from_lower_bound(bmap, i, k);
962 else
963 bmap = set_div_from_lower_bound(bmap, i, l);
964 if (progress)
965 *progress = 1;
966 break;
967 }
968 return bmap;
969 }
970
971 static struct isl_basic_map *remove_duplicate_constraints(
972 struct isl_basic_map *bmap, int *progress)
973 {
974 unsigned int size;
975 isl_int ***index;
976 int k, l, h;
977 int bits;
978 unsigned total = isl_basic_map_total_dim(bmap);
979 isl_int sum;
980
981 if (bmap->n_ineq <= 1)
982 return bmap;
983
984 size = round_up(4 * (bmap->n_ineq+1) / 3 - 1);
985 bits = ffs(size) - 1;
986 index = isl_calloc_array(ctx, isl_int **, size);
987 if (!index)
988 return bmap;
989
990 index[isl_seq_get_hash_bits(bmap->ineq[0]+1, total, bits)] = &bmap->ineq[0];
991 for (k = 1; k < bmap->n_ineq; ++k) {
992 h = hash_index(index, size, bits, bmap, k);
993 if (!index[h]) {
994 index[h] = &bmap->ineq[k];
995 continue;
996 }
997 if (progress)
998 *progress = 1;
999 l = index[h] - &bmap->ineq[0];
1000 if (isl_int_lt(bmap->ineq[k][0], bmap->ineq[l][0]))
1001 swap_inequality(bmap, k, l);
1002 isl_basic_map_drop_inequality(bmap, k);
1003 --k;
1004 }
1005 isl_int_init(sum);
1006 for (k = 0; k < bmap->n_ineq-1; ++k) {
1007 isl_seq_neg(bmap->ineq[k]+1, bmap->ineq[k]+1, total);
1008 h = hash_index(index, size, bits, bmap, k);
1009 isl_seq_neg(bmap->ineq[k]+1, bmap->ineq[k]+1, total);
1010 if (!index[h])
1011 continue;
1012 l = index[h] - &bmap->ineq[0];
1013 isl_int_add(sum, bmap->ineq[k][0], bmap->ineq[l][0]);
1014 if (isl_int_is_pos(sum)) {
1015 bmap = check_for_div_constraints(bmap, k, l, sum,
1016 progress);
1017 continue;
1018 }
1019 if (isl_int_is_zero(sum)) {
1020 /* We need to break out of the loop after these
1021 * changes since the contents of the hash
1022 * will no longer be valid.
1023 * Plus, we probably we want to regauss first.
1024 */
1025 if (progress)
1026 *progress = 1;
1027 isl_basic_map_drop_inequality(bmap, l);
1028 isl_basic_map_inequality_to_equality(bmap, k);
1029 } else
1030 bmap = isl_basic_map_set_to_empty(bmap);
1031 break;
1032 }
1033 isl_int_clear(sum);
1034
1035 free(index);
1036 return bmap;
1037 }
1038
1039
1040 struct isl_basic_map *isl_basic_map_simplify(struct isl_basic_map *bmap)
1041 {
1042 int progress = 1;
1043 if (!bmap)
1044 return NULL;
1045 while (progress) {
1046 progress = 0;
1047 bmap = isl_basic_map_normalize_constraints(bmap);
1048 bmap = remove_duplicate_divs(bmap, &progress);
1049 bmap = eliminate_divs_eq(bmap, &progress);
1050 bmap = eliminate_divs_ineq(bmap, &progress);
1051 bmap = isl_basic_map_gauss(bmap, &progress);
1052 /* requires equalities in normal form */
1053 bmap = normalize_divs(bmap, &progress);
1054 bmap = remove_duplicate_constraints(bmap, &progress);
1055 }
1056 return bmap;
1057 }
1058
1059 struct isl_basic_set *isl_basic_set_simplify(struct isl_basic_set *bset)
1060 {
1061 return (struct isl_basic_set *)
1062 isl_basic_map_simplify((struct isl_basic_map *)bset);
1063 }
1064
1065
1066 /* If the only constraints a div d=floor(f/m)
1067 * appears in are its two defining constraints
1068 *
1069 * f - m d >=0
1070 * -(f - (m - 1)) + m d >= 0
1071 *
1072 * then it can safely be removed.
1073 */
1074 static int div_is_redundant(struct isl_basic_map *bmap, int div)
1075 {
1076 int i;
1077 unsigned pos = 1 + isl_dim_total(bmap->dim) + div;
1078
1079 for (i = 0; i < bmap->n_eq; ++i)
1080 if (!isl_int_is_zero(bmap->eq[i][pos]))
1081 return 0;
1082
1083 for (i = 0; i < bmap->n_ineq; ++i) {
1084 if (isl_int_is_zero(bmap->ineq[i][pos]))
1085 continue;
1086 if (isl_int_eq(bmap->ineq[i][pos], bmap->div[div][0])) {
1087 int neg;
1088 isl_int_sub(bmap->div[div][1],
1089 bmap->div[div][1], bmap->div[div][0]);
1090 isl_int_add_ui(bmap->div[div][1], bmap->div[div][1], 1);
1091 neg = isl_seq_is_neg(bmap->ineq[i], bmap->div[div]+1, pos);
1092 isl_int_sub_ui(bmap->div[div][1], bmap->div[div][1], 1);
1093 isl_int_add(bmap->div[div][1],
1094 bmap->div[div][1], bmap->div[div][0]);
1095 if (!neg)
1096 return 0;
1097 if (isl_seq_first_non_zero(bmap->ineq[i]+pos+1,
1098 bmap->n_div-div-1) != -1)
1099 return 0;
1100 } else if (isl_int_abs_eq(bmap->ineq[i][pos], bmap->div[div][0])) {
1101 if (!isl_seq_eq(bmap->ineq[i], bmap->div[div]+1, pos))
1102 return 0;
1103 if (isl_seq_first_non_zero(bmap->ineq[i]+pos+1,
1104 bmap->n_div-div-1) != -1)
1105 return 0;
1106 } else
1107 return 0;
1108 }
1109
1110 for (i = 0; i < bmap->n_div; ++i)
1111 if (!isl_int_is_zero(bmap->div[i][1+pos]))
1112 return 0;
1113
1114 return 1;
1115 }
1116
1117 /*
1118 * Remove divs that don't occur in any of the constraints or other divs.
1119 * These can arise when dropping some of the variables in a quast
1120 * returned by piplib.
1121 */
1122 static struct isl_basic_map *remove_redundant_divs(struct isl_basic_map *bmap)
1123 {
1124 int i;
1125
1126 if (!bmap)
1127 return NULL;
1128
1129 for (i = bmap->n_div-1; i >= 0; --i) {
1130 if (!div_is_redundant(bmap, i))
1131 continue;
1132 bmap = isl_basic_map_drop_div(bmap, i);
1133 }
1134 return bmap;
1135 }
1136
1137 struct isl_basic_map *isl_basic_map_finalize(struct isl_basic_map *bmap)
1138 {
1139 bmap = remove_redundant_divs(bmap);
1140 if (!bmap)
1141 return NULL;
1142 ISL_F_SET(bmap, ISL_BASIC_SET_FINAL);
1143 return bmap;
1144 }
1145
1146 struct isl_basic_set *isl_basic_set_finalize(struct isl_basic_set *bset)
1147 {
1148 return (struct isl_basic_set *)
1149 isl_basic_map_finalize((struct isl_basic_map *)bset);
1150 }
1151
1152 struct isl_set *isl_set_finalize(struct isl_set *set)
1153 {
1154 int i;
1155
1156 if (!set)
1157 return NULL;
1158 for (i = 0; i < set->n; ++i) {
1159 set->p[i] = isl_basic_set_finalize(set->p[i]);
1160 if (!set->p[i])
1161 goto error;
1162 }
1163 return set;
1164 error:
1165 isl_set_free(set);
1166 return NULL;
1167 }
1168
1169 struct isl_map *isl_map_finalize(struct isl_map *map)
1170 {
1171 int i;
1172
1173 if (!map)
1174 return NULL;
1175 for (i = 0; i < map->n; ++i) {
1176 map->p[i] = isl_basic_map_finalize(map->p[i]);
1177 if (!map->p[i])
1178 goto error;
1179 }
1180 ISL_F_CLR(map, ISL_MAP_NORMALIZED);
1181 return map;
1182 error:
1183 isl_map_free(map);
1184 return NULL;
1185 }
1186
1187
1188 /* Remove definition of any div that is defined in terms of the given variable.
1189 * The div itself is not removed. Functions such as
1190 * eliminate_divs_ineq depend on the other divs remaining in place.
1191 */
1192 static struct isl_basic_map *remove_dependent_vars(struct isl_basic_map *bmap,
1193 int pos)
1194 {
1195 int i;
1196
1197 for (i = 0; i < bmap->n_div; ++i) {
1198 if (isl_int_is_zero(bmap->div[i][0]))
1199 continue;
1200 if (isl_int_is_zero(bmap->div[i][1+1+pos]))
1201 continue;
1202 isl_int_set_si(bmap->div[i][0], 0);
1203 }
1204 return bmap;
1205 }
1206
1207 /* Eliminate the specified variables from the constraints using
1208 * Fourier-Motzkin. The variables themselves are not removed.
1209 */
1210 struct isl_basic_map *isl_basic_map_eliminate_vars(
1211 struct isl_basic_map *bmap, unsigned pos, unsigned n)
1212 {
1213 int d;
1214 int i, j, k;
1215 unsigned total;
1216
1217 if (n == 0)
1218 return bmap;
1219 if (!bmap)
1220 return NULL;
1221 total = isl_basic_map_total_dim(bmap);
1222
1223 bmap = isl_basic_map_cow(bmap);
1224 for (d = pos + n - 1; d >= 0 && d >= pos; --d)
1225 bmap = remove_dependent_vars(bmap, d);
1226
1227 for (d = pos + n - 1;
1228 d >= 0 && d >= total - bmap->n_div && d >= pos; --d)
1229 isl_seq_clr(bmap->div[d-(total-bmap->n_div)], 2+total);
1230 for (d = pos + n - 1; d >= 0 && d >= pos; --d) {
1231 int n_lower, n_upper;
1232 if (!bmap)
1233 return NULL;
1234 for (i = 0; i < bmap->n_eq; ++i) {
1235 if (isl_int_is_zero(bmap->eq[i][1+d]))
1236 continue;
1237 eliminate_var_using_equality(bmap, d, bmap->eq[i], 0, NULL);
1238 isl_basic_map_drop_equality(bmap, i);
1239 break;
1240 }
1241 if (i < bmap->n_eq)
1242 continue;
1243 n_lower = 0;
1244 n_upper = 0;
1245 for (i = 0; i < bmap->n_ineq; ++i) {
1246 if (isl_int_is_pos(bmap->ineq[i][1+d]))
1247 n_lower++;
1248 else if (isl_int_is_neg(bmap->ineq[i][1+d]))
1249 n_upper++;
1250 }
1251 bmap = isl_basic_map_extend_constraints(bmap,
1252 0, n_lower * n_upper);
1253 for (i = bmap->n_ineq - 1; i >= 0; --i) {
1254 int last;
1255 if (isl_int_is_zero(bmap->ineq[i][1+d]))
1256 continue;
1257 last = -1;
1258 for (j = 0; j < i; ++j) {
1259 if (isl_int_is_zero(bmap->ineq[j][1+d]))
1260 continue;
1261 last = j;
1262 if (isl_int_sgn(bmap->ineq[i][1+d]) ==
1263 isl_int_sgn(bmap->ineq[j][1+d]))
1264 continue;
1265 k = isl_basic_map_alloc_inequality(bmap);
1266 if (k < 0)
1267 goto error;
1268 isl_seq_cpy(bmap->ineq[k], bmap->ineq[i],
1269 1+total);
1270 isl_seq_elim(bmap->ineq[k], bmap->ineq[j],
1271 1+d, 1+total, NULL);
1272 }
1273 isl_basic_map_drop_inequality(bmap, i);
1274 i = last + 1;
1275 }
1276 if (n_lower > 0 && n_upper > 0) {
1277 bmap = isl_basic_map_normalize_constraints(bmap);
1278 bmap = remove_duplicate_constraints(bmap, NULL);
1279 bmap = isl_basic_map_gauss(bmap, NULL);
1280 bmap = isl_basic_map_convex_hull(bmap);
1281 if (!bmap)
1282 goto error;
1283 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
1284 break;
1285 }
1286 }
1287 ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
1288 return bmap;
1289 error:
1290 isl_basic_map_free(bmap);
1291 return NULL;
1292 }
1293
1294 struct isl_basic_set *isl_basic_set_eliminate_vars(
1295 struct isl_basic_set *bset, unsigned pos, unsigned n)
1296 {
1297 return (struct isl_basic_set *)isl_basic_map_eliminate_vars(
1298 (struct isl_basic_map *)bset, pos, n);
1299 }
1300
1301 /* Don't assume equalities are in order, because align_divs
1302 * may have changed the order of the divs.
1303 */
1304 static void compute_elimination_index(struct isl_basic_map *bmap, int *elim)
1305 {
1306 int d, i;
1307 unsigned total;
1308
1309 total = isl_dim_total(bmap->dim);
1310 for (d = 0; d < total; ++d)
1311 elim[d] = -1;
1312 for (i = 0; i < bmap->n_eq; ++i) {
1313 for (d = total - 1; d >= 0; --d) {
1314 if (isl_int_is_zero(bmap->eq[i][1+d]))
1315 continue;
1316 elim[d] = i;
1317 break;
1318 }
1319 }
1320 }
1321
1322 static void set_compute_elimination_index(struct isl_basic_set *bset, int *elim)
1323 {
1324 compute_elimination_index((struct isl_basic_map *)bset, elim);
1325 }
1326
1327 static int reduced_using_equalities(isl_int *dst, isl_int *src,
1328 struct isl_basic_map *bmap, int *elim)
1329 {
1330 int d;
1331 int copied = 0;
1332 unsigned total;
1333
1334 total = isl_dim_total(bmap->dim);
1335 for (d = total - 1; d >= 0; --d) {
1336 if (isl_int_is_zero(src[1+d]))
1337 continue;
1338 if (elim[d] == -1)
1339 continue;
1340 if (!copied) {
1341 isl_seq_cpy(dst, src, 1 + total);
1342 copied = 1;
1343 }
1344 isl_seq_elim(dst, bmap->eq[elim[d]], 1 + d, 1 + total, NULL);
1345 }
1346 return copied;
1347 }
1348
1349 static int set_reduced_using_equalities(isl_int *dst, isl_int *src,
1350 struct isl_basic_set *bset, int *elim)
1351 {
1352 return reduced_using_equalities(dst, src,
1353 (struct isl_basic_map *)bset, elim);
1354 }
1355
1356 static struct isl_basic_set *isl_basic_set_reduce_using_equalities(
1357 struct isl_basic_set *bset, struct isl_basic_set *context)
1358 {
1359 int i;
1360 int *elim;
1361
1362 if (!bset || !context)
1363 goto error;
1364
1365 bset = isl_basic_set_cow(bset);
1366 if (!bset)
1367 goto error;
1368
1369 elim = isl_alloc_array(ctx, int, isl_basic_set_n_dim(bset));
1370 if (!elim)
1371 goto error;
1372 set_compute_elimination_index(context, elim);
1373 for (i = 0; i < bset->n_eq; ++i)
1374 set_reduced_using_equalities(bset->eq[i], bset->eq[i],
1375 context, elim);
1376 for (i = 0; i < bset->n_ineq; ++i)
1377 set_reduced_using_equalities(bset->ineq[i], bset->ineq[i],
1378 context, elim);
1379 isl_basic_set_free(context);
1380 free(elim);
1381 bset = isl_basic_set_simplify(bset);
1382 bset = isl_basic_set_finalize(bset);
1383 return bset;
1384 error:
1385 isl_basic_set_free(bset);
1386 isl_basic_set_free(context);
1387 return NULL;
1388 }
1389
1390 static struct isl_basic_set *remove_shifted_constraints(
1391 struct isl_basic_set *bset, struct isl_basic_set *context)
1392 {
1393 unsigned int size;
1394 isl_int ***index;
1395 int bits;
1396 int k, h, l;
1397
1398 if (!bset)
1399 return NULL;
1400
1401 size = round_up(4 * (context->n_ineq+1) / 3 - 1);
1402 bits = ffs(size) - 1;
1403 index = isl_calloc_array(ctx, isl_int **, size);
1404 if (!index)
1405 return bset;
1406
1407 for (k = 0; k < context->n_ineq; ++k) {
1408 h = set_hash_index(index, size, bits, context, k);
1409 index[h] = &context->ineq[k];
1410 }
1411 for (k = 0; k < bset->n_ineq; ++k) {
1412 h = set_hash_index(index, size, bits, bset, k);
1413 if (!index[h])
1414 continue;
1415 l = index[h] - &context->ineq[0];
1416 if (isl_int_lt(bset->ineq[k][0], context->ineq[l][0]))
1417 continue;
1418 bset = isl_basic_set_cow(bset);
1419 if (!bset)
1420 goto error;
1421 isl_basic_set_drop_inequality(bset, k);
1422 --k;
1423 }
1424 free(index);
1425 return bset;
1426 error:
1427 free(index);
1428 return bset;
1429 }
1430
1431 /* Tighten (decrease) the constant terms of the inequalities based
1432 * on the equalities, without removing any integer points.
1433 * For example, if there is an equality
1434 *
1435 * i = 3 * j
1436 *
1437 * and an inequality
1438 *
1439 * i >= 1
1440 *
1441 * then we want to replace the inequality by
1442 *
1443 * i >= 3
1444 *
1445 * We do this by computing a variable compression and translating
1446 * the constraints to the compressed space.
1447 * If any constraint has coefficients (except the contant term)
1448 * with a common factor "f", then we can replace the constant term "c"
1449 * by
1450 *
1451 * f * floor(c/f)
1452 *
1453 * That is, we add
1454 *
1455 * f * floor(c/f) - c = -fract(c/f)
1456 *
1457 * and we can add the same value to the original constraint.
1458 *
1459 * In the example, the compressed space only contains "j",
1460 * and the inequality translates to
1461 *
1462 * 3 * j - 1 >= 0
1463 *
1464 * We add -fract(-1/3) = -2 to the original constraint to obtain
1465 *
1466 * i - 3 >= 0
1467 */
1468 static struct isl_basic_set *normalize_constraints_in_compressed_space(
1469 struct isl_basic_set *bset)
1470 {
1471 int i;
1472 unsigned total;
1473 struct isl_mat *B, *C;
1474 isl_int gcd;
1475
1476 if (!bset)
1477 return NULL;
1478
1479 if (ISL_F_ISSET(bset, ISL_BASIC_SET_RATIONAL))
1480 return bset;
1481
1482 if (!bset->n_ineq)
1483 return bset;
1484
1485 bset = isl_basic_set_cow(bset);
1486 if (!bset)
1487 return NULL;
1488
1489 total = isl_basic_set_total_dim(bset);
1490 B = isl_mat_sub_alloc(bset->ctx, bset->eq, 0, bset->n_eq, 0, 1 + total);
1491 C = isl_mat_variable_compression(B, NULL);
1492 if (!C)
1493 return bset;
1494 if (C->n_col == 0) {
1495 isl_mat_free(C);
1496 return isl_basic_set_set_to_empty(bset);
1497 }
1498 B = isl_mat_sub_alloc(bset->ctx, bset->ineq,
1499 0, bset->n_ineq, 0, 1 + total);
1500 C = isl_mat_product(B, C);
1501 if (!C)
1502 return bset;
1503
1504 isl_int_init(gcd);
1505 for (i = 0; i < bset->n_ineq; ++i) {
1506 isl_seq_gcd(C->row[i] + 1, C->n_col - 1, &gcd);
1507 if (isl_int_is_one(gcd))
1508 continue;
1509 isl_int_fdiv_r(C->row[i][0], C->row[i][0], gcd);
1510 isl_int_sub(bset->ineq[i][0], bset->ineq[i][0], C->row[i][0]);
1511 }
1512 isl_int_clear(gcd);
1513
1514 isl_mat_free(C);
1515
1516 return bset;
1517 }
1518
1519 /* Remove all information from bset that is redundant in the context
1520 * of context. In particular, equalities that are linear combinations
1521 * of those in context are removed. Then the inequalities that are
1522 * redundant in the context of the equalities and inequalities of
1523 * context are removed.
1524 *
1525 * We first simplify the constraints of "bset" in the context of the
1526 * equalities of "context".
1527 * Then we simplify the inequalities of the context in the context
1528 * of the equalities of bset and remove the inequalities from "bset"
1529 * that are obviously redundant with respect to some inequality in "context".
1530 *
1531 * If there are any inequalities left, we construct a tableau for
1532 * the context and then add the inequalities of "bset".
1533 * Before adding these equalities, we freeze all constraints such that
1534 * they won't be considered redundant in terms of the constraints of "bset".
1535 * Then we detect all equalities and redundant constraints (among the
1536 * constraints that weren't frozen) and update bset according to the results.
1537 * We have to be careful here because we don't want any of the context
1538 * constraints to remain and because we haven't added the equalities of "bset"
1539 * to the tableau so we temporarily have to pretend that there were no
1540 * equalities.
1541 */
1542 static struct isl_basic_set *uset_gist(struct isl_basic_set *bset,
1543 struct isl_basic_set *context)
1544 {
1545 int i;
1546 struct isl_tab *tab;
1547 unsigned context_ineq;
1548 struct isl_basic_set *combined = NULL;
1549
1550 if (!context || !bset)
1551 goto error;
1552
1553 if (context->n_eq > 0)
1554 bset = isl_basic_set_reduce_using_equalities(bset,
1555 isl_basic_set_copy(context));
1556 if (!bset)
1557 goto error;
1558 if (isl_basic_set_fast_is_empty(bset))
1559 goto done;
1560 if (!bset->n_ineq)
1561 goto done;
1562
1563 if (bset->n_eq > 0) {
1564 struct isl_basic_set *affine_hull;
1565 affine_hull = isl_basic_set_copy(bset);
1566 affine_hull = isl_basic_set_cow(affine_hull);
1567 if (!affine_hull)
1568 goto error;
1569 isl_basic_set_free_inequality(affine_hull, affine_hull->n_ineq);
1570 context = isl_basic_set_intersect(context, affine_hull);
1571 context = isl_basic_set_gauss(context, NULL);
1572 context = normalize_constraints_in_compressed_space(context);
1573 }
1574 if (!context)
1575 goto error;
1576 if (ISL_F_ISSET(context, ISL_BASIC_SET_EMPTY)) {
1577 isl_basic_set_free(bset);
1578 return context;
1579 }
1580 if (!context->n_ineq)
1581 goto done;
1582 bset = remove_shifted_constraints(bset, context);
1583 if (!bset->n_ineq)
1584 goto done;
1585 context_ineq = context->n_ineq;
1586 combined = isl_basic_set_cow(isl_basic_set_copy(context));
1587 if (isl_basic_set_free_equality(combined, context->n_eq) < 0)
1588 goto error;
1589 combined = isl_basic_set_extend_constraints(combined,
1590 bset->n_eq, bset->n_ineq);
1591 tab = isl_tab_from_basic_set(combined);
1592 if (!tab)
1593 goto error;
1594 for (i = 0; i < context_ineq; ++i)
1595 if (isl_tab_freeze_constraint(tab, i) < 0)
1596 goto error;
1597 tab = isl_tab_extend(tab, bset->n_ineq);
1598 if (!tab)
1599 goto error;
1600 for (i = 0; i < bset->n_ineq; ++i)
1601 if (isl_tab_add_ineq(tab, bset->ineq[i]) < 0)
1602 goto error;
1603 bset = isl_basic_set_add_constraints(combined, bset, 0);
1604 tab = isl_tab_detect_implicit_equalities(tab);
1605 if (isl_tab_detect_redundant(tab) < 0) {
1606 isl_tab_free(tab);
1607 goto error2;
1608 }
1609 for (i = 0; i < context_ineq; ++i) {
1610 tab->con[i].is_zero = 0;
1611 tab->con[i].is_redundant = 1;
1612 }
1613 bset = isl_basic_set_update_from_tab(bset, tab);
1614 isl_tab_free(tab);
1615 ISL_F_SET(bset, ISL_BASIC_SET_NO_IMPLICIT);
1616 ISL_F_SET(bset, ISL_BASIC_SET_NO_REDUNDANT);
1617 done:
1618 bset = isl_basic_set_simplify(bset);
1619 bset = isl_basic_set_finalize(bset);
1620 isl_basic_set_free(context);
1621 return bset;
1622 error:
1623 isl_basic_set_free(combined);
1624 error2:
1625 isl_basic_set_free(bset);
1626 isl_basic_set_free(context);
1627 return NULL;
1628 }
1629
1630 /* Normalize the divs in "bmap" in the context of the equalities in "context".
1631 * We simply add the equalities in context to bmap and then do a regular
1632 * div normalizations. Better results can be obtained by normalizing
1633 * only the divs in bmap than do not also appear in context.
1634 * We need to be careful to reduce the divs using the equalities
1635 * so that later calls to isl_basic_map_overlying_set wouldn't introduce
1636 * spurious constraints.
1637 */
1638 static struct isl_basic_map *normalize_divs_in_context(
1639 struct isl_basic_map *bmap, struct isl_basic_map *context)
1640 {
1641 int i;
1642 unsigned total_context;
1643 int div_eq;
1644
1645 div_eq = n_pure_div_eq(bmap);
1646 if (div_eq == 0)
1647 return bmap;
1648
1649 if (context->n_div > 0)
1650 bmap = isl_basic_map_align_divs(bmap, context);
1651
1652 total_context = isl_basic_map_total_dim(context);
1653 bmap = isl_basic_map_extend_constraints(bmap, context->n_eq, 0);
1654 for (i = 0; i < context->n_eq; ++i) {
1655 int k;
1656 k = isl_basic_map_alloc_equality(bmap);
1657 isl_seq_cpy(bmap->eq[k], context->eq[i], 1 + total_context);
1658 isl_seq_clr(bmap->eq[k] + 1 + total_context,
1659 isl_basic_map_total_dim(bmap) - total_context);
1660 }
1661 bmap = isl_basic_map_gauss(bmap, NULL);
1662 bmap = normalize_divs(bmap, NULL);
1663 bmap = isl_basic_map_gauss(bmap, NULL);
1664 return bmap;
1665 }
1666
1667 struct isl_basic_map *isl_basic_map_gist(struct isl_basic_map *bmap,
1668 struct isl_basic_map *context)
1669 {
1670 struct isl_basic_set *bset;
1671
1672 if (!bmap || !context)
1673 goto error;
1674
1675 if (isl_basic_map_is_universe(context)) {
1676 isl_basic_map_free(context);
1677 return bmap;
1678 }
1679 if (isl_basic_map_is_universe(bmap)) {
1680 isl_basic_map_free(context);
1681 return bmap;
1682 }
1683 if (isl_basic_map_fast_is_empty(context)) {
1684 struct isl_dim *dim = isl_dim_copy(bmap->dim);
1685 isl_basic_map_free(context);
1686 isl_basic_map_free(bmap);
1687 return isl_basic_map_universe(dim);
1688 }
1689 if (isl_basic_map_fast_is_empty(bmap)) {
1690 isl_basic_map_free(context);
1691 return bmap;
1692 }
1693
1694 bmap = isl_basic_map_convex_hull(bmap);
1695 context = isl_basic_map_convex_hull(context);
1696
1697 if (context->n_eq)
1698 bmap = normalize_divs_in_context(bmap, context);
1699
1700 context = isl_basic_map_align_divs(context, bmap);
1701 bmap = isl_basic_map_align_divs(bmap, context);
1702
1703 bset = uset_gist(isl_basic_map_underlying_set(isl_basic_map_copy(bmap)),
1704 isl_basic_map_underlying_set(context));
1705
1706 return isl_basic_map_overlying_set(bset, bmap);
1707 error:
1708 isl_basic_map_free(bmap);
1709 isl_basic_map_free(context);
1710 return NULL;
1711 }
1712
1713 /*
1714 * Assumes context has no implicit divs.
1715 */
1716 struct isl_map *isl_map_gist(struct isl_map *map, struct isl_basic_map *context)
1717 {
1718 int i;
1719
1720 if (!map || !context)
1721 goto error;;
1722
1723 if (isl_basic_map_is_universe(context)) {
1724 isl_basic_map_free(context);
1725 return map;
1726 }
1727 if (isl_basic_map_fast_is_empty(context)) {
1728 struct isl_dim *dim = isl_dim_copy(map->dim);
1729 isl_basic_map_free(context);
1730 isl_map_free(map);
1731 return isl_map_universe(dim);
1732 }
1733
1734 context = isl_basic_map_convex_hull(context);
1735 map = isl_map_cow(map);
1736 if (!map || !context)
1737 goto error;;
1738 isl_assert(map->ctx, isl_dim_equal(map->dim, context->dim), goto error);
1739 map = isl_map_compute_divs(map);
1740 for (i = 0; i < map->n; ++i)
1741 context = isl_basic_map_align_divs(context, map->p[i]);
1742 for (i = 0; i < map->n; ++i) {
1743 map->p[i] = isl_basic_map_gist(map->p[i],
1744 isl_basic_map_copy(context));
1745 if (!map->p[i])
1746 goto error;
1747 }
1748 isl_basic_map_free(context);
1749 ISL_F_CLR(map, ISL_MAP_NORMALIZED);
1750 return map;
1751 error:
1752 isl_map_free(map);
1753 isl_basic_map_free(context);
1754 return NULL;
1755 }
1756
1757 struct isl_basic_set *isl_basic_set_gist(struct isl_basic_set *bset,
1758 struct isl_basic_set *context)
1759 {
1760 return (struct isl_basic_set *)isl_basic_map_gist(
1761 (struct isl_basic_map *)bset, (struct isl_basic_map *)context);
1762 }
1763
1764 struct isl_set *isl_set_gist(struct isl_set *set, struct isl_basic_set *context)
1765 {
1766 return (struct isl_set *)isl_map_gist((struct isl_map *)set,
1767 (struct isl_basic_map *)context);
1768 }
1769
1770 /* Quick check to see if two basic maps are disjoint.
1771 * In particular, we reduce the equalities and inequalities of
1772 * one basic map in the context of the equalities of the other
1773 * basic map and check if we get a contradiction.
1774 */
1775 int isl_basic_map_fast_is_disjoint(struct isl_basic_map *bmap1,
1776 struct isl_basic_map *bmap2)
1777 {
1778 struct isl_vec *v = NULL;
1779 int *elim = NULL;
1780 unsigned total;
1781 int i;
1782
1783 if (!bmap1 || !bmap2)
1784 return -1;
1785 isl_assert(bmap1->ctx, isl_dim_equal(bmap1->dim, bmap2->dim),
1786 return -1);
1787 if (bmap1->n_div || bmap2->n_div)
1788 return 0;
1789 if (!bmap1->n_eq && !bmap2->n_eq)
1790 return 0;
1791
1792 total = isl_dim_total(bmap1->dim);
1793 if (total == 0)
1794 return 0;
1795 v = isl_vec_alloc(bmap1->ctx, 1 + total);
1796 if (!v)
1797 goto error;
1798 elim = isl_alloc_array(bmap1->ctx, int, total);
1799 if (!elim)
1800 goto error;
1801 compute_elimination_index(bmap1, elim);
1802 for (i = 0; i < bmap2->n_eq; ++i) {
1803 int reduced;
1804 reduced = reduced_using_equalities(v->block.data, bmap2->eq[i],
1805 bmap1, elim);
1806 if (reduced && !isl_int_is_zero(v->block.data[0]) &&
1807 isl_seq_first_non_zero(v->block.data + 1, total) == -1)
1808 goto disjoint;
1809 }
1810 for (i = 0; i < bmap2->n_ineq; ++i) {
1811 int reduced;
1812 reduced = reduced_using_equalities(v->block.data,
1813 bmap2->ineq[i], bmap1, elim);
1814 if (reduced && isl_int_is_neg(v->block.data[0]) &&
1815 isl_seq_first_non_zero(v->block.data + 1, total) == -1)
1816 goto disjoint;
1817 }
1818 compute_elimination_index(bmap2, elim);
1819 for (i = 0; i < bmap1->n_ineq; ++i) {
1820 int reduced;
1821 reduced = reduced_using_equalities(v->block.data,
1822 bmap1->ineq[i], bmap2, elim);
1823 if (reduced && isl_int_is_neg(v->block.data[0]) &&
1824 isl_seq_first_non_zero(v->block.data + 1, total) == -1)
1825 goto disjoint;
1826 }
1827 isl_vec_free(v);
1828 free(elim);
1829 return 0;
1830 disjoint:
1831 isl_vec_free(v);
1832 free(elim);
1833 return 1;
1834 error:
1835 isl_vec_free(v);
1836 free(elim);
1837 return -1;
1838 }
1839
1840 int isl_basic_set_fast_is_disjoint(struct isl_basic_set *bset1,
1841 struct isl_basic_set *bset2)
1842 {
1843 return isl_basic_map_fast_is_disjoint((struct isl_basic_map *)bset1,
1844 (struct isl_basic_map *)bset2);
1845 }
1846
1847 int isl_map_fast_is_disjoint(struct isl_map *map1, struct isl_map *map2)
1848 {
1849 int i, j;
1850
1851 if (!map1 || !map2)
1852 return -1;
1853
1854 if (isl_map_fast_is_equal(map1, map2))
1855 return 0;
1856
1857 for (i = 0; i < map1->n; ++i) {
1858 for (j = 0; j < map2->n; ++j) {
1859 int d = isl_basic_map_fast_is_disjoint(map1->p[i],
1860 map2->p[j]);
1861 if (d != 1)
1862 return d;
1863 }
1864 }
1865 return 1;
1866 }
1867
1868 int isl_set_fast_is_disjoint(struct isl_set *set1, struct isl_set *set2)
1869 {
1870 return isl_map_fast_is_disjoint((struct isl_map *)set1,
1871 (struct isl_map *)set2);
1872 }
1873
1874 /* Check if we can combine a given div with lower bound l and upper
1875 * bound u with some other div and if so return that other div.
1876 * Otherwise return -1.
1877 *
1878 * We first check that
1879 * - the bounds are opposites of each other (except for the constant
1880 * term)
1881 * - the bounds do not reference any other div
1882 * - no div is defined in terms of this div
1883 *
1884 * Let m be the size of the range allowed on the div by the bounds.
1885 * That is, the bounds are of the form
1886 *
1887 * e <= a <= e + m - 1
1888 *
1889 * with e some expression in the other variables.
1890 * We look for another div b such that no third div is defined in terms
1891 * of this second div b and such that in any constraint that contains
1892 * a (except for the given lower and upper bound), also contains b
1893 * with a coefficient that is m times that of b.
1894 * That is, all constraints (execpt for the lower and upper bound)
1895 * are of the form
1896 *
1897 * e + f (a + m b) >= 0
1898 *
1899 * If so, we return b so that "a + m b" can be replaced by
1900 * a single div "c = a + m b".
1901 */
1902 static int div_find_coalesce(struct isl_basic_map *bmap, int *pairs,
1903 unsigned div, unsigned l, unsigned u)
1904 {
1905 int i, j;
1906 unsigned dim;
1907 int coalesce = -1;
1908
1909 if (bmap->n_div <= 1)
1910 return -1;
1911 dim = isl_dim_total(bmap->dim);
1912 if (isl_seq_first_non_zero(bmap->ineq[l] + 1 + dim, div) != -1)
1913 return -1;
1914 if (isl_seq_first_non_zero(bmap->ineq[l] + 1 + dim + div + 1,
1915 bmap->n_div - div - 1) != -1)
1916 return -1;
1917 if (!isl_seq_is_neg(bmap->ineq[l] + 1, bmap->ineq[u] + 1,
1918 dim + bmap->n_div))
1919 return -1;
1920
1921 for (i = 0; i < bmap->n_div; ++i) {
1922 if (isl_int_is_zero(bmap->div[i][0]))
1923 continue;
1924 if (!isl_int_is_zero(bmap->div[i][1 + 1 + dim + div]))
1925 return -1;
1926 }
1927
1928 isl_int_add(bmap->ineq[l][0], bmap->ineq[l][0], bmap->ineq[u][0]);
1929 if (isl_int_is_neg(bmap->ineq[l][0])) {
1930 isl_int_sub(bmap->ineq[l][0],
1931 bmap->ineq[l][0], bmap->ineq[u][0]);
1932 bmap = isl_basic_map_copy(bmap);
1933 bmap = isl_basic_map_set_to_empty(bmap);
1934 isl_basic_map_free(bmap);
1935 return -1;
1936 }
1937 isl_int_add_ui(bmap->ineq[l][0], bmap->ineq[l][0], 1);
1938 for (i = 0; i < bmap->n_div; ++i) {
1939 if (i == div)
1940 continue;
1941 if (!pairs[i])
1942 continue;
1943 for (j = 0; j < bmap->n_div; ++j) {
1944 if (isl_int_is_zero(bmap->div[j][0]))
1945 continue;
1946 if (!isl_int_is_zero(bmap->div[j][1 + 1 + dim + i]))
1947 break;
1948 }
1949 if (j < bmap->n_div)
1950 continue;
1951 for (j = 0; j < bmap->n_ineq; ++j) {
1952 int valid;
1953 if (j == l || j == u)
1954 continue;
1955 if (isl_int_is_zero(bmap->ineq[j][1 + dim + div]))
1956 continue;
1957 if (isl_int_is_zero(bmap->ineq[j][1 + dim + i]))
1958 break;
1959 isl_int_mul(bmap->ineq[j][1 + dim + div],
1960 bmap->ineq[j][1 + dim + div],
1961 bmap->ineq[l][0]);
1962 valid = isl_int_eq(bmap->ineq[j][1 + dim + div],
1963 bmap->ineq[j][1 + dim + i]);
1964 isl_int_divexact(bmap->ineq[j][1 + dim + div],
1965 bmap->ineq[j][1 + dim + div],
1966 bmap->ineq[l][0]);
1967 if (!valid)
1968 break;
1969 }
1970 if (j < bmap->n_ineq)
1971 continue;
1972 coalesce = i;
1973 break;
1974 }
1975 isl_int_sub_ui(bmap->ineq[l][0], bmap->ineq[l][0], 1);
1976 isl_int_sub(bmap->ineq[l][0], bmap->ineq[l][0], bmap->ineq[u][0]);
1977 return coalesce;
1978 }
1979
1980 /* Given a lower and an upper bound on div i, construct an inequality
1981 * that when nonnegative ensures that this pair of bounds always allows
1982 * for an integer value of the given div.
1983 * The lower bound is inequality l, while the upper bound is inequality u.
1984 * The constructed inequality is stored in ineq.
1985 * g, fl, fu are temporary scalars.
1986 *
1987 * Let the upper bound be
1988 *
1989 * -n_u a + e_u >= 0
1990 *
1991 * and the lower bound
1992 *
1993 * n_l a + e_l >= 0
1994 *
1995 * Let n_u = f_u g and n_l = f_l g, with g = gcd(n_u, n_l).
1996 * We have
1997 *
1998 * - f_u e_l <= f_u f_l g a <= f_l e_u
1999 *
2000 * Since all variables are integer valued, this is equivalent to
2001 *
2002 * - f_u e_l - (f_u - 1) <= f_u f_l g a <= f_l e_u + (f_l - 1)
2003 *
2004 * If this interval is at least f_u f_l g, then it contains at least
2005 * one integer value for a.
2006 * That is, the test constraint is
2007 *
2008 * f_l e_u + f_u e_l + f_l - 1 + f_u - 1 + 1 >= f_u f_l g
2009 */
2010 static void construct_test_ineq(struct isl_basic_map *bmap, int i,
2011 int l, int u, isl_int *ineq, isl_int g, isl_int fl, isl_int fu)
2012 {
2013 unsigned dim;
2014 dim = isl_dim_total(bmap->dim);
2015
2016 isl_int_gcd(g, bmap->ineq[l][1 + dim + i], bmap->ineq[u][1 + dim + i]);
2017 isl_int_divexact(fl, bmap->ineq[l][1 + dim + i], g);
2018 isl_int_divexact(fu, bmap->ineq[u][1 + dim + i], g);
2019 isl_int_neg(fu, fu);
2020 isl_seq_combine(ineq, fl, bmap->ineq[u], fu, bmap->ineq[l],
2021 1 + dim + bmap->n_div);
2022 isl_int_add(ineq[0], ineq[0], fl);
2023 isl_int_add(ineq[0], ineq[0], fu);
2024 isl_int_sub_ui(ineq[0], ineq[0], 1);
2025 isl_int_mul(g, g, fl);
2026 isl_int_mul(g, g, fu);
2027 isl_int_sub(ineq[0], ineq[0], g);
2028 }
2029
2030 /* Remove more kinds of divs that are not strictly needed.
2031 * In particular, if all pairs of lower and upper bounds on a div
2032 * are such that they allow at least one integer value of the div,
2033 * the we can eliminate the div using Fourier-Motzkin without
2034 * introducing any spurious solutions.
2035 */
2036 static struct isl_basic_map *drop_more_redundant_divs(
2037 struct isl_basic_map *bmap, int *pairs, int n)
2038 {
2039 struct isl_tab *tab = NULL;
2040 struct isl_vec *vec = NULL;
2041 unsigned dim;
2042 int remove = -1;
2043 isl_int g, fl, fu;
2044
2045 isl_int_init(g);
2046 isl_int_init(fl);
2047 isl_int_init(fu);
2048
2049 if (!bmap)
2050 goto error;
2051
2052 dim = isl_dim_total(bmap->dim);
2053 vec = isl_vec_alloc(bmap->ctx, 1 + dim + bmap->n_div);
2054 if (!vec)
2055 goto error;
2056
2057 tab = isl_tab_from_basic_map(bmap);
2058
2059 while (n > 0) {
2060 int i, l, u;
2061 int best = -1;
2062 enum isl_lp_result res;
2063
2064 for (i = 0; i < bmap->n_div; ++i) {
2065 if (!pairs[i])
2066 continue;
2067 if (best >= 0 && pairs[best] <= pairs[i])
2068 continue;
2069 best = i;
2070 }
2071
2072 i = best;
2073 for (l = 0; l < bmap->n_ineq; ++l) {
2074 if (!isl_int_is_pos(bmap->ineq[l][1 + dim + i]))
2075 continue;
2076 for (u = 0; u < bmap->n_ineq; ++u) {
2077 if (!isl_int_is_neg(bmap->ineq[u][1 + dim + i]))
2078 continue;
2079 construct_test_ineq(bmap, i, l, u,
2080 vec->el, g, fl, fu);
2081 res = isl_tab_min(tab, vec->el,
2082 bmap->ctx->one, &g, NULL, 0);
2083 if (res == isl_lp_error)
2084 goto error;
2085 if (res == isl_lp_empty) {
2086 bmap = isl_basic_map_set_to_empty(bmap);
2087 break;
2088 }
2089 if (res != isl_lp_ok || isl_int_is_neg(g))
2090 break;
2091 }
2092 if (u < bmap->n_ineq)
2093 break;
2094 }
2095 if (l == bmap->n_ineq) {
2096 remove = i;
2097 break;
2098 }
2099 pairs[i] = 0;
2100 --n;
2101 }
2102
2103 isl_tab_free(tab);
2104 isl_vec_free(vec);
2105
2106 isl_int_clear(g);
2107 isl_int_clear(fl);
2108 isl_int_clear(fu);
2109
2110 free(pairs);
2111
2112 if (remove < 0)
2113 return bmap;
2114
2115 bmap = isl_basic_map_remove(bmap, isl_dim_div, remove, 1);
2116 return isl_basic_map_drop_redundant_divs(bmap);
2117 error:
2118 free(pairs);
2119 isl_basic_map_free(bmap);
2120 isl_tab_free(tab);
2121 isl_vec_free(vec);
2122 isl_int_clear(g);
2123 isl_int_clear(fl);
2124 isl_int_clear(fu);
2125 return NULL;
2126 }
2127
2128 /* Given a pair of divs div1 and div2 such that, expect for the lower bound l
2129 * and the upper bound u, div1 always occurs together with div2 in the form
2130 * (div1 + m div2), where m is the constant range on the variable div1
2131 * allowed by l and u, replace the pair div1 and div2 by a single
2132 * div that is equal to div1 + m div2.
2133 *
2134 * The new div will appear in the location that contains div2.
2135 * We need to modify all constraints that contain
2136 * div2 = (div - div1) / m
2137 * (If a constraint does not contain div2, it will also not contain div1.)
2138 * If the constraint also contains div1, then we know they appear
2139 * as f (div1 + m div2) and we can simply replace (div1 + m div2) by div,
2140 * i.e., the coefficient of div is f.
2141 *
2142 * Otherwise, we first need to introduce div1 into the constraint.
2143 * Let the l be
2144 *
2145 * div1 + f >=0
2146 *
2147 * and u
2148 *
2149 * -div1 + f' >= 0
2150 *
2151 * A lower bound on div2
2152 *
2153 * n div2 + t >= 0
2154 *
2155 * can be replaced by
2156 *
2157 * (n * (m div 2 + div1) + m t + n f)/g >= 0
2158 *
2159 * with g = gcd(m,n).
2160 * An upper bound
2161 *
2162 * -n div2 + t >= 0
2163 *
2164 * can be replaced by
2165 *
2166 * (-n * (m div2 + div1) + m t + n f')/g >= 0
2167 *
2168 * These constraint are those that we would obtain from eliminating
2169 * div1 using Fourier-Motzkin.
2170 *
2171 * After all constraints have been modified, we drop the lower and upper
2172 * bound and then drop div1.
2173 */
2174 static struct isl_basic_map *coalesce_divs(struct isl_basic_map *bmap,
2175 unsigned div1, unsigned div2, unsigned l, unsigned u)
2176 {
2177 isl_int a;
2178 isl_int b;
2179 isl_int m;
2180 unsigned dim, total;
2181 int i;
2182
2183 dim = isl_dim_total(bmap->dim);
2184 total = 1 + dim + bmap->n_div;
2185
2186 isl_int_init(a);
2187 isl_int_init(b);
2188 isl_int_init(m);
2189 isl_int_add(m, bmap->ineq[l][0], bmap->ineq[u][0]);
2190 isl_int_add_ui(m, m, 1);
2191
2192 for (i = 0; i < bmap->n_ineq; ++i) {
2193 if (i == l || i == u)
2194 continue;
2195 if (isl_int_is_zero(bmap->ineq[i][1 + dim + div2]))
2196 continue;
2197 if (isl_int_is_zero(bmap->ineq[i][1 + dim + div1])) {
2198 isl_int_gcd(b, m, bmap->ineq[i][1 + dim + div2]);
2199 isl_int_divexact(a, m, b);
2200 isl_int_divexact(b, bmap->ineq[i][1 + dim + div2], b);
2201 if (isl_int_is_pos(b)) {
2202 isl_seq_combine(bmap->ineq[i], a, bmap->ineq[i],
2203 b, bmap->ineq[l], total);
2204 } else {
2205 isl_int_neg(b, b);
2206 isl_seq_combine(bmap->ineq[i], a, bmap->ineq[i],
2207 b, bmap->ineq[u], total);
2208 }
2209 }
2210 isl_int_set(bmap->ineq[i][1 + dim + div2],
2211 bmap->ineq[i][1 + dim + div1]);
2212 isl_int_set_si(bmap->ineq[i][1 + dim + div1], 0);
2213 }
2214
2215 isl_int_clear(a);
2216 isl_int_clear(b);
2217 isl_int_clear(m);
2218 if (l > u) {
2219 isl_basic_map_drop_inequality(bmap, l);
2220 isl_basic_map_drop_inequality(bmap, u);
2221 } else {
2222 isl_basic_map_drop_inequality(bmap, u);
2223 isl_basic_map_drop_inequality(bmap, l);
2224 }
2225 bmap = isl_basic_map_drop_div(bmap, div1);
2226 return bmap;
2227 }
2228
2229 /* First check if we can coalesce any pair of divs and
2230 * then continue with dropping more redundant divs.
2231 *
2232 * We loop over all pairs of lower and upper bounds on a div
2233 * with coefficient 1 and -1, respectively, check if there
2234 * is any other div "c" with which we can coalesce the div
2235 * and if so, perform the coalescing.
2236 */
2237 static struct isl_basic_map *coalesce_or_drop_more_redundant_divs(
2238 struct isl_basic_map *bmap, int *pairs, int n)
2239 {
2240 int i, l, u;
2241 unsigned dim;
2242
2243 dim = isl_dim_total(bmap->dim);
2244
2245 for (i = 0; i < bmap->n_div; ++i) {
2246 if (!pairs[i])
2247 continue;
2248 for (l = 0; l < bmap->n_ineq; ++l) {
2249 if (!isl_int_is_one(bmap->ineq[l][1 + dim + i]))
2250 continue;
2251 for (u = 0; u < bmap->n_ineq; ++u) {
2252 int c;
2253
2254 if (!isl_int_is_negone(bmap->ineq[u][1+dim+i]))
2255 continue;
2256 c = div_find_coalesce(bmap, pairs, i, l, u);
2257 if (c < 0)
2258 continue;
2259 free(pairs);
2260 bmap = coalesce_divs(bmap, i, c, l, u);
2261 return isl_basic_map_drop_redundant_divs(bmap);
2262 }
2263 }
2264 }
2265
2266 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
2267 return bmap;
2268
2269 return drop_more_redundant_divs(bmap, pairs, n);
2270 }
2271
2272 /* Remove divs that are not strictly needed.
2273 * In particular, if a div only occurs positively (or negatively)
2274 * in constraints, then it can simply be dropped.
2275 * Also, if a div occurs only occurs in two constraints and if moreover
2276 * those two constraints are opposite to each other, except for the constant
2277 * term and if the sum of the constant terms is such that for any value
2278 * of the other values, there is always at least one integer value of the
2279 * div, i.e., if one plus this sum is greater than or equal to
2280 * the (absolute value) of the coefficent of the div in the constraints,
2281 * then we can also simply drop the div.
2282 *
2283 * If any divs are left after these simple checks then we move on
2284 * to more complicated cases in drop_more_redundant_divs.
2285 */
2286 struct isl_basic_map *isl_basic_map_drop_redundant_divs(
2287 struct isl_basic_map *bmap)
2288 {
2289 int i, j;
2290 unsigned off;
2291 int *pairs = NULL;
2292 int n = 0;
2293
2294 if (!bmap)
2295 goto error;
2296
2297 off = isl_dim_total(bmap->dim);
2298 pairs = isl_calloc_array(bmap->ctx, int, bmap->n_div);
2299 if (!pairs)
2300 goto error;
2301
2302 for (i = 0; i < bmap->n_div; ++i) {
2303 int pos, neg;
2304 int last_pos, last_neg;
2305 int redundant;
2306 int defined;
2307
2308 defined = !isl_int_is_zero(bmap->div[i][0]);
2309 for (j = 0; j < bmap->n_eq; ++j)
2310 if (!isl_int_is_zero(bmap->eq[j][1 + off + i]))
2311 break;
2312 if (j < bmap->n_eq)
2313 continue;
2314 ++n;
2315 pos = neg = 0;
2316 for (j = 0; j < bmap->n_ineq; ++j) {
2317 if (isl_int_is_pos(bmap->ineq[j][1 + off + i])) {
2318 last_pos = j;
2319 ++pos;
2320 }
2321 if (isl_int_is_neg(bmap->ineq[j][1 + off + i])) {
2322 last_neg = j;
2323 ++neg;
2324 }
2325 }
2326 pairs[i] = pos * neg;
2327 if (pairs[i] == 0) {
2328 for (j = bmap->n_ineq - 1; j >= 0; --j)
2329 if (!isl_int_is_zero(bmap->ineq[j][1+off+i]))
2330 isl_basic_map_drop_inequality(bmap, j);
2331 bmap = isl_basic_map_drop_div(bmap, i);
2332 free(pairs);
2333 return isl_basic_map_drop_redundant_divs(bmap);
2334 }
2335 if (pairs[i] != 1)
2336 continue;
2337 if (!isl_seq_is_neg(bmap->ineq[last_pos] + 1,
2338 bmap->ineq[last_neg] + 1,
2339 off + bmap->n_div))
2340 continue;
2341
2342 isl_int_add(bmap->ineq[last_pos][0],
2343 bmap->ineq[last_pos][0], bmap->ineq[last_neg][0]);
2344 isl_int_add_ui(bmap->ineq[last_pos][0],
2345 bmap->ineq[last_pos][0], 1);
2346 redundant = isl_int_ge(bmap->ineq[last_pos][0],
2347 bmap->ineq[last_pos][1+off+i]);
2348 isl_int_sub_ui(bmap->ineq[last_pos][0],
2349 bmap->ineq[last_pos][0], 1);
2350 isl_int_sub(bmap->ineq[last_pos][0],
2351 bmap->ineq[last_pos][0], bmap->ineq[last_neg][0]);
2352 if (!redundant) {
2353 if (defined ||
2354 !ok_to_set_div_from_bound(bmap, i, last_pos)) {
2355 pairs[i] = 0;
2356 --n;
2357 continue;
2358 }
2359 bmap = set_div_from_lower_bound(bmap, i, last_pos);
2360 bmap = isl_basic_map_simplify(bmap);
2361 free(pairs);
2362 return isl_basic_map_drop_redundant_divs(bmap);
2363 }
2364 if (last_pos > last_neg) {
2365 isl_basic_map_drop_inequality(bmap, last_pos);
2366 isl_basic_map_drop_inequality(bmap, last_neg);
2367 } else {
2368 isl_basic_map_drop_inequality(bmap, last_neg);
2369 isl_basic_map_drop_inequality(bmap, last_pos);
2370 }
2371 bmap = isl_basic_map_drop_div(bmap, i);
2372 free(pairs);
2373 return isl_basic_map_drop_redundant_divs(bmap);
2374 }
2375
2376 if (n > 0)
2377 return coalesce_or_drop_more_redundant_divs(bmap, pairs, n);
2378
2379 free(pairs);
2380 return bmap;
2381 error:
2382 free(pairs);
2383 isl_basic_map_free(bmap);
2384 return NULL;
2385 }
2386
2387 struct isl_basic_set *isl_basic_set_drop_redundant_divs(
2388 struct isl_basic_set *bset)
2389 {
2390 return (struct isl_basic_set *)
2391 isl_basic_map_drop_redundant_divs((struct isl_basic_map *)bset);
2392 }
2393
2394 struct isl_map *isl_map_drop_redundant_divs(struct isl_map *map)
2395 {
2396 int i;
2397
2398 if (!map)
2399 return NULL;
2400 for (i = 0; i < map->n; ++i) {
2401 map->p[i] = isl_basic_map_drop_redundant_divs(map->p[i]);
2402 if (!map->p[i])
2403 goto error;
2404 }
2405 ISL_F_CLR(map, ISL_MAP_NORMALIZED);
2406 return map;
2407 error:
2408 isl_map_free(map);
2409 return NULL;
2410 }
2411
2412 struct isl_set *isl_set_drop_redundant_divs(struct isl_set *set)
2413 {
2414 return (struct isl_set *)
2415 isl_map_drop_redundant_divs((struct isl_map *)set);
2416 }
Integer set library