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