Line data Source code
1 : /*
2 : * Copyright 2008-2009 Katholieke Universiteit Leuven
3 : * Copyright 2012-2013 Ecole Normale Superieure
4 : * Copyright 2014-2015 INRIA Rocquencourt
5 : * Copyright 2016 Sven Verdoolaege
6 : *
7 : * Use of this software is governed by the MIT license
8 : *
9 : * Written by Sven Verdoolaege, K.U.Leuven, Departement
10 : * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
11 : * and Ecole Normale Superieure, 45 rue d’Ulm, 75230 Paris, France
12 : * and Inria Paris - Rocquencourt, Domaine de Voluceau - Rocquencourt,
13 : * B.P. 105 - 78153 Le Chesnay, France
14 : */
15 :
16 : #include <isl_ctx_private.h>
17 : #include <isl_map_private.h>
18 : #include "isl_equalities.h"
19 : #include <isl/map.h>
20 : #include <isl_seq.h>
21 : #include "isl_tab.h"
22 : #include <isl_space_private.h>
23 : #include <isl_mat_private.h>
24 : #include <isl_vec_private.h>
25 :
26 : #include <bset_to_bmap.c>
27 : #include <bset_from_bmap.c>
28 : #include <set_to_map.c>
29 : #include <set_from_map.c>
30 :
31 10848825273 : static void swap_equality(struct isl_basic_map *bmap, int a, int b)
32 : {
33 10848825273 : isl_int *t = bmap->eq[a];
34 10848825273 : bmap->eq[a] = bmap->eq[b];
35 10848825273 : bmap->eq[b] = t;
36 10848825273 : }
37 :
38 30940916 : static void swap_inequality(struct isl_basic_map *bmap, int a, int b)
39 : {
40 30940916 : if (a != b) {
41 30940916 : isl_int *t = bmap->ineq[a];
42 30940916 : bmap->ineq[a] = bmap->ineq[b];
43 30940916 : bmap->ineq[b] = t;
44 : }
45 30940916 : }
46 :
47 19194043916 : __isl_give isl_basic_map *isl_basic_map_normalize_constraints(
48 : __isl_take isl_basic_map *bmap)
49 : {
50 : int i;
51 : isl_int gcd;
52 19194043916 : unsigned total = isl_basic_map_total_dim(bmap);
53 :
54 19194043916 : if (!bmap)
55 0 : return NULL;
56 :
57 19194043916 : isl_int_init(gcd);
58 54241642462 : for (i = bmap->n_eq - 1; i >= 0; --i) {
59 35523001973 : isl_seq_gcd(bmap->eq[i]+1, total, &gcd);
60 35523001973 : if (isl_int_is_zero(gcd)) {
61 1576581836 : if (!isl_int_is_zero(bmap->eq[i][0])) {
62 474390804 : bmap = isl_basic_map_set_to_empty(bmap);
63 474390804 : break;
64 : }
65 1102191032 : isl_basic_map_drop_equality(bmap, i);
66 1102191032 : continue;
67 : }
68 33946420137 : if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
69 33543505272 : isl_int_gcd(gcd, gcd, bmap->eq[i][0]);
70 33946420137 : if (isl_int_is_one(gcd))
71 30540955301 : continue;
72 3405464836 : if (!isl_int_is_divisible_by(bmap->eq[i][0], gcd)) {
73 1012623 : bmap = isl_basic_map_set_to_empty(bmap);
74 1012623 : break;
75 : }
76 3404452213 : isl_seq_scale_down(bmap->eq[i], bmap->eq[i], gcd, 1+total);
77 : }
78 :
79 33137036911 : for (i = bmap->n_ineq - 1; i >= 0; --i) {
80 13962593570 : isl_seq_gcd(bmap->ineq[i]+1, total, &gcd);
81 13962593570 : if (isl_int_is_zero(gcd)) {
82 1309153033 : if (isl_int_is_neg(bmap->ineq[i][0])) {
83 19600575 : bmap = isl_basic_map_set_to_empty(bmap);
84 19600575 : break;
85 : }
86 1289552458 : isl_basic_map_drop_inequality(bmap, i);
87 1289552458 : continue;
88 : }
89 12653440537 : if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
90 9497125210 : isl_int_gcd(gcd, gcd, bmap->ineq[i][0]);
91 12653440537 : if (isl_int_is_one(gcd))
92 12582698644 : continue;
93 70741893 : isl_int_fdiv_q(bmap->ineq[i][0], bmap->ineq[i][0], gcd);
94 70741893 : isl_seq_scale_down(bmap->ineq[i]+1, bmap->ineq[i]+1, gcd, total);
95 : }
96 19194043916 : isl_int_clear(gcd);
97 :
98 19194043916 : return bmap;
99 : }
100 :
101 2979789573 : __isl_give isl_basic_set *isl_basic_set_normalize_constraints(
102 : __isl_take isl_basic_set *bset)
103 : {
104 2979789573 : isl_basic_map *bmap = bset_to_bmap(bset);
105 2979789573 : return bset_from_bmap(isl_basic_map_normalize_constraints(bmap));
106 : }
107 :
108 : /* Reduce the coefficient of the variable at position "pos"
109 : * in integer division "div", such that it lies in the half-open
110 : * interval (1/2,1/2], extracting any excess value from this integer division.
111 : * "pos" is as determined by isl_basic_map_offset, i.e., pos == 0
112 : * corresponds to the constant term.
113 : *
114 : * That is, the integer division is of the form
115 : *
116 : * floor((... + (c * d + r) * x_pos + ...)/d)
117 : *
118 : * with -d < 2 * r <= d.
119 : * Replace it by
120 : *
121 : * floor((... + r * x_pos + ...)/d) + c * x_pos
122 : *
123 : * If 2 * ((c * d + r) % d) <= d, then c = floor((c * d + r)/d).
124 : * Otherwise, c = floor((c * d + r)/d) + 1.
125 : *
126 : * This is the same normalization that is performed by isl_aff_floor.
127 : */
128 0 : static __isl_give isl_basic_map *reduce_coefficient_in_div(
129 : __isl_take isl_basic_map *bmap, int div, int pos)
130 : {
131 : isl_int shift;
132 : int add_one;
133 :
134 0 : isl_int_init(shift);
135 0 : isl_int_fdiv_r(shift, bmap->div[div][1 + pos], bmap->div[div][0]);
136 0 : isl_int_mul_ui(shift, shift, 2);
137 0 : add_one = isl_int_gt(shift, bmap->div[div][0]);
138 0 : isl_int_fdiv_q(shift, bmap->div[div][1 + pos], bmap->div[div][0]);
139 0 : if (add_one)
140 0 : isl_int_add_ui(shift, shift, 1);
141 0 : isl_int_neg(shift, shift);
142 0 : bmap = isl_basic_map_shift_div(bmap, div, pos, shift);
143 0 : isl_int_clear(shift);
144 :
145 0 : return bmap;
146 : }
147 :
148 : /* Does the coefficient of the variable at position "pos"
149 : * in integer division "div" need to be reduced?
150 : * That is, does it lie outside the half-open interval (1/2,1/2]?
151 : * The coefficient c/d lies outside this interval if abs(2 * c) >= d and
152 : * 2 * c != d.
153 : */
154 0 : static isl_bool needs_reduction(__isl_keep isl_basic_map *bmap, int div,
155 : int pos)
156 : {
157 : isl_bool r;
158 :
159 0 : if (isl_int_is_zero(bmap->div[div][1 + pos]))
160 0 : return isl_bool_false;
161 :
162 0 : isl_int_mul_ui(bmap->div[div][1 + pos], bmap->div[div][1 + pos], 2);
163 0 : r = isl_int_abs_ge(bmap->div[div][1 + pos], bmap->div[div][0]) &&
164 0 : !isl_int_eq(bmap->div[div][1 + pos], bmap->div[div][0]);
165 0 : isl_int_divexact_ui(bmap->div[div][1 + pos],
166 : bmap->div[div][1 + pos], 2);
167 :
168 0 : return r;
169 : }
170 :
171 : /* Reduce the coefficients (including the constant term) of
172 : * integer division "div", if needed.
173 : * In particular, make sure all coefficients lie in
174 : * the half-open interval (1/2,1/2].
175 : */
176 0 : static __isl_give isl_basic_map *reduce_div_coefficients_of_div(
177 : __isl_take isl_basic_map *bmap, int div)
178 : {
179 : int i;
180 0 : unsigned total = 1 + isl_basic_map_total_dim(bmap);
181 :
182 0 : for (i = 0; i < total; ++i) {
183 : isl_bool reduce;
184 :
185 0 : reduce = needs_reduction(bmap, div, i);
186 0 : if (reduce < 0)
187 0 : return isl_basic_map_free(bmap);
188 0 : if (!reduce)
189 0 : continue;
190 0 : bmap = reduce_coefficient_in_div(bmap, div, i);
191 0 : if (!bmap)
192 0 : break;
193 : }
194 :
195 0 : return bmap;
196 : }
197 :
198 : /* Reduce the coefficients (including the constant term) of
199 : * the known integer divisions, if needed
200 : * In particular, make sure all coefficients lie in
201 : * the half-open interval (1/2,1/2].
202 : */
203 16214254343 : static __isl_give isl_basic_map *reduce_div_coefficients(
204 : __isl_take isl_basic_map *bmap)
205 : {
206 : int i;
207 :
208 16214254343 : if (!bmap)
209 0 : return NULL;
210 16214254343 : if (bmap->n_div == 0)
211 16214254343 : return bmap;
212 :
213 0 : for (i = 0; i < bmap->n_div; ++i) {
214 0 : if (isl_int_is_zero(bmap->div[i][0]))
215 0 : continue;
216 0 : bmap = reduce_div_coefficients_of_div(bmap, i);
217 0 : if (!bmap)
218 0 : break;
219 : }
220 :
221 0 : return bmap;
222 : }
223 :
224 : /* Remove any common factor in numerator and denominator of the div expression,
225 : * not taking into account the constant term.
226 : * That is, if the div is of the form
227 : *
228 : * floor((a + m f(x))/(m d))
229 : *
230 : * then replace it by
231 : *
232 : * floor((floor(a/m) + f(x))/d)
233 : *
234 : * The difference {a/m}/d in the argument satisfies 0 <= {a/m}/d < 1/d
235 : * and can therefore not influence the result of the floor.
236 : */
237 0 : static void normalize_div_expression(__isl_keep isl_basic_map *bmap, int div)
238 : {
239 0 : unsigned total = isl_basic_map_total_dim(bmap);
240 0 : isl_ctx *ctx = bmap->ctx;
241 :
242 0 : if (isl_int_is_zero(bmap->div[div][0]))
243 0 : return;
244 0 : isl_seq_gcd(bmap->div[div] + 2, total, &ctx->normalize_gcd);
245 0 : isl_int_gcd(ctx->normalize_gcd, ctx->normalize_gcd, bmap->div[div][0]);
246 0 : if (isl_int_is_one(ctx->normalize_gcd))
247 0 : return;
248 0 : isl_int_fdiv_q(bmap->div[div][1], bmap->div[div][1],
249 : ctx->normalize_gcd);
250 0 : isl_int_divexact(bmap->div[div][0], bmap->div[div][0],
251 : ctx->normalize_gcd);
252 0 : isl_seq_scale_down(bmap->div[div] + 2, bmap->div[div] + 2,
253 0 : ctx->normalize_gcd, total);
254 : }
255 :
256 : /* Remove any common factor in numerator and denominator of a div expression,
257 : * not taking into account the constant term.
258 : * That is, look for any div of the form
259 : *
260 : * floor((a + m f(x))/(m d))
261 : *
262 : * and replace it by
263 : *
264 : * floor((floor(a/m) + f(x))/d)
265 : *
266 : * The difference {a/m}/d in the argument satisfies 0 <= {a/m}/d < 1/d
267 : * and can therefore not influence the result of the floor.
268 : */
269 16214254343 : static __isl_give isl_basic_map *normalize_div_expressions(
270 : __isl_take isl_basic_map *bmap)
271 : {
272 : int i;
273 :
274 16214254343 : if (!bmap)
275 0 : return NULL;
276 16214254343 : if (bmap->n_div == 0)
277 16214254343 : return bmap;
278 :
279 0 : for (i = 0; i < bmap->n_div; ++i)
280 0 : normalize_div_expression(bmap, i);
281 :
282 0 : return bmap;
283 : }
284 :
285 : /* Assumes divs have been ordered if keep_divs is set.
286 : */
287 62338595215 : static void eliminate_var_using_equality(struct isl_basic_map *bmap,
288 : unsigned pos, isl_int *eq, int keep_divs, int *progress)
289 : {
290 : unsigned total;
291 : unsigned space_total;
292 : int k;
293 : int last_div;
294 :
295 62338595215 : total = isl_basic_map_total_dim(bmap);
296 62338595215 : space_total = isl_space_dim(bmap->dim, isl_dim_all);
297 62338595215 : last_div = isl_seq_last_non_zero(eq + 1 + space_total, bmap->n_div);
298 >34603*10^7 : for (k = 0; k < bmap->n_eq; ++k) {
299 >28369*10^7 : if (bmap->eq[k] == eq)
300 62338595215 : continue;
301 >22135*10^7 : if (isl_int_is_zero(bmap->eq[k][1+pos]))
302 >19484*10^7 : continue;
303 26511451875 : if (progress)
304 10108038752 : *progress = 1;
305 26511451875 : isl_seq_elim(bmap->eq[k], eq, 1+pos, 1+total, NULL);
306 26511451875 : isl_seq_normalize(bmap->ctx, bmap->eq[k], 1 + total);
307 : }
308 :
309 >11834*10^7 : for (k = 0; k < bmap->n_ineq; ++k) {
310 56006708240 : if (isl_int_is_zero(bmap->ineq[k][1+pos]))
311 48459366240 : continue;
312 7547342000 : if (progress)
313 489853963 : *progress = 1;
314 7547342000 : isl_seq_elim(bmap->ineq[k], eq, 1+pos, 1+total, NULL);
315 7547342000 : isl_seq_normalize(bmap->ctx, bmap->ineq[k], 1 + total);
316 7547342000 : ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
317 : }
318 :
319 62338595215 : for (k = 0; k < bmap->n_div; ++k) {
320 0 : if (isl_int_is_zero(bmap->div[k][0]))
321 0 : continue;
322 0 : if (isl_int_is_zero(bmap->div[k][1+1+pos]))
323 0 : continue;
324 0 : if (progress)
325 0 : *progress = 1;
326 : /* We need to be careful about circular definitions,
327 : * so for now we just remove the definition of div k
328 : * if the equality contains any divs.
329 : * If keep_divs is set, then the divs have been ordered
330 : * and we can keep the definition as long as the result
331 : * is still ordered.
332 : */
333 0 : if (last_div == -1 || (keep_divs && last_div < k)) {
334 0 : isl_seq_elim(bmap->div[k]+1, eq,
335 0 : 1+pos, 1+total, &bmap->div[k][0]);
336 0 : normalize_div_expression(bmap, k);
337 : } else
338 0 : isl_seq_clr(bmap->div[k], 1 + total);
339 0 : ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
340 : }
341 62338595215 : }
342 :
343 : /* Assumes divs have been ordered if keep_divs is set.
344 : */
345 0 : static __isl_give isl_basic_map *eliminate_div(__isl_take isl_basic_map *bmap,
346 : isl_int *eq, unsigned div, int keep_divs)
347 : {
348 0 : unsigned pos = isl_space_dim(bmap->dim, isl_dim_all) + div;
349 :
350 0 : eliminate_var_using_equality(bmap, pos, eq, keep_divs, NULL);
351 :
352 0 : bmap = isl_basic_map_drop_div(bmap, div);
353 :
354 0 : return bmap;
355 : }
356 :
357 : /* Check if elimination of div "div" using equality "eq" would not
358 : * result in a div depending on a later div.
359 : */
360 0 : static isl_bool ok_to_eliminate_div(__isl_keep isl_basic_map *bmap, isl_int *eq,
361 : unsigned div)
362 : {
363 : int k;
364 : int last_div;
365 0 : unsigned space_total = isl_space_dim(bmap->dim, isl_dim_all);
366 0 : unsigned pos = space_total + div;
367 :
368 0 : last_div = isl_seq_last_non_zero(eq + 1 + space_total, bmap->n_div);
369 0 : if (last_div < 0 || last_div <= div)
370 0 : return isl_bool_true;
371 :
372 0 : for (k = 0; k <= last_div; ++k) {
373 0 : if (isl_int_is_zero(bmap->div[k][0]))
374 0 : continue;
375 0 : if (!isl_int_is_zero(bmap->div[k][1 + 1 + pos]))
376 0 : return isl_bool_false;
377 : }
378 :
379 0 : return isl_bool_true;
380 : }
381 :
382 : /* Eliminate divs based on equalities
383 : */
384 16214254561 : static __isl_give isl_basic_map *eliminate_divs_eq(
385 : __isl_take isl_basic_map *bmap, int *progress)
386 : {
387 : int d;
388 : int i;
389 16214254561 : int modified = 0;
390 : unsigned off;
391 :
392 16214254561 : bmap = isl_basic_map_order_divs(bmap);
393 :
394 16214254561 : if (!bmap)
395 0 : return NULL;
396 :
397 16214254561 : off = 1 + isl_space_dim(bmap->dim, isl_dim_all);
398 :
399 16214254561 : for (d = bmap->n_div - 1; d >= 0 ; --d) {
400 0 : for (i = 0; i < bmap->n_eq; ++i) {
401 : isl_bool ok;
402 :
403 0 : if (!isl_int_is_one(bmap->eq[i][off + d]) &&
404 0 : !isl_int_is_negone(bmap->eq[i][off + d]))
405 0 : continue;
406 0 : ok = ok_to_eliminate_div(bmap, bmap->eq[i], d);
407 0 : if (ok < 0)
408 0 : return isl_basic_map_free(bmap);
409 0 : if (!ok)
410 0 : continue;
411 0 : modified = 1;
412 0 : *progress = 1;
413 0 : bmap = eliminate_div(bmap, bmap->eq[i], d, 1);
414 0 : if (isl_basic_map_drop_equality(bmap, i) < 0)
415 0 : return isl_basic_map_free(bmap);
416 0 : break;
417 : }
418 : }
419 16214254561 : if (modified)
420 0 : return eliminate_divs_eq(bmap, progress);
421 16214254561 : return bmap;
422 : }
423 :
424 : /* Eliminate divs based on inequalities
425 : */
426 16214254343 : static __isl_give isl_basic_map *eliminate_divs_ineq(
427 : __isl_take isl_basic_map *bmap, int *progress)
428 : {
429 : int d;
430 : int i;
431 : unsigned off;
432 : struct isl_ctx *ctx;
433 :
434 16214254343 : if (!bmap)
435 0 : return NULL;
436 :
437 16214254343 : ctx = bmap->ctx;
438 16214254343 : off = 1 + isl_space_dim(bmap->dim, isl_dim_all);
439 :
440 16214254343 : for (d = bmap->n_div - 1; d >= 0 ; --d) {
441 0 : for (i = 0; i < bmap->n_eq; ++i)
442 0 : if (!isl_int_is_zero(bmap->eq[i][off + d]))
443 : break;
444 0 : if (i < bmap->n_eq)
445 0 : continue;
446 0 : for (i = 0; i < bmap->n_ineq; ++i)
447 0 : if (isl_int_abs_gt(bmap->ineq[i][off + d], ctx->one))
448 0 : break;
449 0 : if (i < bmap->n_ineq)
450 0 : continue;
451 0 : *progress = 1;
452 0 : bmap = isl_basic_map_eliminate_vars(bmap, (off-1)+d, 1);
453 0 : if (!bmap || ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
454 : break;
455 0 : bmap = isl_basic_map_drop_div(bmap, d);
456 0 : if (!bmap)
457 0 : break;
458 : }
459 16214254343 : return bmap;
460 : }
461 :
462 : /* Does the equality constraint at position "eq" in "bmap" involve
463 : * any local variables in the range [first, first + n)
464 : * that are not marked as having an explicit representation?
465 : */
466 0 : static isl_bool bmap_eq_involves_unknown_divs(__isl_keep isl_basic_map *bmap,
467 : int eq, unsigned first, unsigned n)
468 : {
469 : unsigned o_div;
470 : int i;
471 :
472 0 : if (!bmap)
473 0 : return isl_bool_error;
474 :
475 0 : o_div = isl_basic_map_offset(bmap, isl_dim_div);
476 0 : for (i = 0; i < n; ++i) {
477 : isl_bool unknown;
478 :
479 0 : if (isl_int_is_zero(bmap->eq[eq][o_div + first + i]))
480 0 : continue;
481 0 : unknown = isl_basic_map_div_is_marked_unknown(bmap, first + i);
482 0 : if (unknown < 0)
483 0 : return isl_bool_error;
484 0 : if (unknown)
485 0 : return isl_bool_true;
486 : }
487 :
488 0 : return isl_bool_false;
489 : }
490 :
491 : /* The last local variable involved in the equality constraint
492 : * at position "eq" in "bmap" is the local variable at position "div".
493 : * It can therefore be used to extract an explicit representation
494 : * for that variable.
495 : * Do so unless the local variable already has an explicit representation or
496 : * the explicit representation would involve any other local variables
497 : * that in turn do not have an explicit representation.
498 : * An equality constraint involving local variables without an explicit
499 : * representation can be used in isl_basic_map_drop_redundant_divs
500 : * to separate out an independent local variable. Introducing
501 : * an explicit representation here would block this transformation,
502 : * while the partial explicit representation in itself is not very useful.
503 : * Set *progress if anything is changed.
504 : *
505 : * The equality constraint is of the form
506 : *
507 : * f(x) + n e >= 0
508 : *
509 : * with n a positive number. The explicit representation derived from
510 : * this constraint is
511 : *
512 : * floor((-f(x))/n)
513 : */
514 0 : static __isl_give isl_basic_map *set_div_from_eq(__isl_take isl_basic_map *bmap,
515 : int div, int eq, int *progress)
516 : {
517 : unsigned total, o_div;
518 : isl_bool involves;
519 :
520 0 : if (!bmap)
521 0 : return NULL;
522 :
523 0 : if (!isl_int_is_zero(bmap->div[div][0]))
524 0 : return bmap;
525 :
526 0 : involves = bmap_eq_involves_unknown_divs(bmap, eq, 0, div);
527 0 : if (involves < 0)
528 0 : return isl_basic_map_free(bmap);
529 0 : if (involves)
530 0 : return bmap;
531 :
532 0 : total = isl_basic_map_dim(bmap, isl_dim_all);
533 0 : o_div = isl_basic_map_offset(bmap, isl_dim_div);
534 0 : isl_seq_neg(bmap->div[div] + 1, bmap->eq[eq], 1 + total);
535 0 : isl_int_set_si(bmap->div[div][1 + o_div + div], 0);
536 0 : isl_int_set(bmap->div[div][0], bmap->eq[eq][o_div + div]);
537 0 : if (progress)
538 0 : *progress = 1;
539 0 : ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
540 :
541 0 : return bmap;
542 : }
543 :
544 : /* Perform fangcheng (Gaussian elimination) on the equality
545 : * constraints of "bmap".
546 : * That is, put them into row-echelon form, starting from the last column
547 : * backward and use them to eliminate the corresponding coefficients
548 : * from all constraints.
549 : *
550 : * If "progress" is not NULL, then it gets set if the elimination
551 : * result in any changes.
552 : * The elimination process may result in some equality constraints
553 : * getting interchanged or removed.
554 : * If "swap" or "drop" are not NULL, then they get called when
555 : * two equality constraints get interchanged or
556 : * when a number of final equality constraints get removed.
557 : * As a special case, if the input turns out to be empty,
558 : * then drop gets called with the number of removed equality
559 : * constraints set to the total number of equality constraints.
560 : * If "swap" or "drop" are not NULL, then the local variables (if any)
561 : * are assumed to be in a valid order.
562 : */
563 30448952040 : __isl_give isl_basic_map *isl_basic_map_gauss5(__isl_take isl_basic_map *bmap,
564 : int *progress,
565 : isl_stat (*swap)(unsigned a, unsigned b, void *user),
566 : isl_stat (*drop)(unsigned n, void *user), void *user)
567 : {
568 : int k;
569 : int done;
570 : int last_var;
571 : unsigned total_var;
572 : unsigned total;
573 : unsigned n_drop;
574 :
575 30448952040 : if (!swap && !drop)
576 30448899974 : bmap = isl_basic_map_order_divs(bmap);
577 :
578 30448952040 : if (!bmap)
579 0 : return NULL;
580 :
581 30448952040 : total = isl_basic_map_total_dim(bmap);
582 30448952040 : total_var = total - bmap->n_div;
583 :
584 30448952040 : last_var = total - 1;
585 92786271998 : for (done = 0; done < bmap->n_eq; ++done) {
586 >11071*10^7 : for (; last_var >= 0; --last_var) {
587 >28195*10^7 : for (k = done; k < bmap->n_eq; ++k)
588 >23890*10^7 : if (!isl_int_is_zero(bmap->eq[k][1+last_var]))
589 : break;
590 >10538*10^7 : if (k < bmap->n_eq)
591 62337319958 : break;
592 : }
593 67668298380 : if (last_var < 0)
594 5330978422 : break;
595 62337319958 : if (k != done) {
596 10848825273 : swap_equality(bmap, k, done);
597 10848825273 : if (swap && swap(k, done, user) < 0)
598 0 : return isl_basic_map_free(bmap);
599 : }
600 62337319958 : if (isl_int_is_neg(bmap->eq[done][1+last_var]))
601 5745305914 : isl_seq_neg(bmap->eq[done], bmap->eq[done], 1+total);
602 :
603 62337319958 : eliminate_var_using_equality(bmap, last_var, bmap->eq[done], 1,
604 : progress);
605 :
606 62337319958 : if (last_var >= total_var)
607 0 : bmap = set_div_from_eq(bmap, last_var - total_var,
608 : done, progress);
609 62337319958 : if (!bmap)
610 0 : return NULL;
611 : }
612 30448952040 : if (done == bmap->n_eq)
613 25117973618 : return bmap;
614 17476516574 : for (k = done; k < bmap->n_eq; ++k) {
615 5683201980 : if (isl_int_is_zero(bmap->eq[k][0]))
616 3407279865 : continue;
617 2275922115 : if (drop && drop(bmap->n_eq, user) < 0)
618 0 : return isl_basic_map_free(bmap);
619 2275922115 : return isl_basic_map_set_to_empty(bmap);
620 : }
621 3055056307 : n_drop = bmap->n_eq - done;
622 3055056307 : isl_basic_map_free_equality(bmap, n_drop);
623 3055056307 : if (drop && drop(n_drop, user) < 0)
624 0 : return isl_basic_map_free(bmap);
625 3055056307 : return bmap;
626 : }
627 :
628 30448899974 : __isl_give isl_basic_map *isl_basic_map_gauss(__isl_take isl_basic_map *bmap,
629 : int *progress)
630 : {
631 30448899974 : return isl_basic_map_gauss5(bmap, progress, NULL, NULL, NULL);
632 : }
633 :
634 2458845818 : __isl_give isl_basic_set *isl_basic_set_gauss(
635 : __isl_take isl_basic_set *bset, int *progress)
636 : {
637 2458845818 : return bset_from_bmap(isl_basic_map_gauss(bset_to_bmap(bset),
638 : progress));
639 : }
640 :
641 :
642 1992419621 : static unsigned int round_up(unsigned int v)
643 : {
644 1992419621 : int old_v = v;
645 :
646 7532214623 : while (v) {
647 3547375381 : old_v = v;
648 3547375381 : v ^= v & -v;
649 : }
650 1992419621 : return old_v << 1;
651 : }
652 :
653 : /* Hash table of inequalities in a basic map.
654 : * "index" is an array of addresses of inequalities in the basic map, some
655 : * of which are NULL. The inequalities are hashed on the coefficients
656 : * except the constant term.
657 : * "size" is the number of elements in the array and is always a power of two
658 : * "bits" is the number of bits need to represent an index into the array.
659 : * "total" is the total dimension of the basic map.
660 : */
661 : struct isl_constraint_index {
662 : unsigned int size;
663 : int bits;
664 : isl_int ***index;
665 : unsigned total;
666 : };
667 :
668 : /* Fill in the "ci" data structure for holding the inequalities of "bmap".
669 : */
670 1992419621 : static isl_stat create_constraint_index(struct isl_constraint_index *ci,
671 : __isl_keep isl_basic_map *bmap)
672 : {
673 : isl_ctx *ctx;
674 :
675 1992419621 : ci->index = NULL;
676 1992419621 : if (!bmap)
677 0 : return isl_stat_error;
678 1992419621 : ci->total = isl_basic_set_total_dim(bmap);
679 1992419621 : if (bmap->n_ineq == 0)
680 0 : return isl_stat_ok;
681 1992419621 : ci->size = round_up(4 * (bmap->n_ineq + 1) / 3 - 1);
682 1992419621 : ci->bits = ffs(ci->size) - 1;
683 1992419621 : ctx = isl_basic_map_get_ctx(bmap);
684 1992419621 : ci->index = isl_calloc_array(ctx, isl_int **, ci->size);
685 1992419621 : if (!ci->index)
686 0 : return isl_stat_error;
687 :
688 1992419621 : return isl_stat_ok;
689 : }
690 :
691 : /* Free the memory allocated by create_constraint_index.
692 : */
693 1992419621 : static void constraint_index_free(struct isl_constraint_index *ci)
694 : {
695 1992419621 : free(ci->index);
696 1992419621 : }
697 :
698 : /* Return the position in ci->index that contains the address of
699 : * an inequality that is equal to *ineq up to the constant term,
700 : * provided this address is not identical to "ineq".
701 : * If there is no such inequality, then return the position where
702 : * such an inequality should be inserted.
703 : */
704 10103746287 : static int hash_index_ineq(struct isl_constraint_index *ci, isl_int **ineq)
705 : {
706 : int h;
707 10103746287 : uint32_t hash = isl_seq_get_hash_bits((*ineq) + 1, ci->total, ci->bits);
708 15149639054 : for (h = hash; ci->index[h]; h = (h+1) % ci->size)
709 12170524815 : if (ineq != ci->index[h] &&
710 5825045912 : isl_seq_eq((*ineq) + 1, ci->index[h][0]+1, ci->total))
711 1299586136 : break;
712 10103746287 : return h;
713 : }
714 :
715 : /* Return the position in ci->index that contains the address of
716 : * an inequality that is equal to the k'th inequality of "bmap"
717 : * up to the constant term, provided it does not point to the very
718 : * same inequality.
719 : * If there is no such inequality, then return the position where
720 : * such an inequality should be inserted.
721 : */
722 10103746287 : static int hash_index(struct isl_constraint_index *ci,
723 : __isl_keep isl_basic_map *bmap, int k)
724 : {
725 10103746287 : return hash_index_ineq(ci, &bmap->ineq[k]);
726 : }
727 :
728 0 : static int set_hash_index(struct isl_constraint_index *ci,
729 : __isl_keep isl_basic_set *bset, int k)
730 : {
731 0 : return hash_index(ci, bset, k);
732 : }
733 :
734 : /* Fill in the "ci" data structure with the inequalities of "bset".
735 : */
736 0 : static isl_stat setup_constraint_index(struct isl_constraint_index *ci,
737 : __isl_keep isl_basic_set *bset)
738 : {
739 : int k, h;
740 :
741 0 : if (create_constraint_index(ci, bset) < 0)
742 0 : return isl_stat_error;
743 :
744 0 : for (k = 0; k < bset->n_ineq; ++k) {
745 0 : h = set_hash_index(ci, bset, k);
746 0 : ci->index[h] = &bset->ineq[k];
747 : }
748 :
749 0 : return isl_stat_ok;
750 : }
751 :
752 : /* Is the inequality ineq (obviously) redundant with respect
753 : * to the constraints in "ci"?
754 : *
755 : * Look for an inequality in "ci" with the same coefficients and then
756 : * check if the contant term of "ineq" is greater than or equal
757 : * to the constant term of that inequality. If so, "ineq" is clearly
758 : * redundant.
759 : *
760 : * Note that hash_index_ineq ignores a stored constraint if it has
761 : * the same address as the passed inequality. It is ok to pass
762 : * the address of a local variable here since it will never be
763 : * the same as the address of a constraint in "ci".
764 : */
765 0 : static isl_bool constraint_index_is_redundant(struct isl_constraint_index *ci,
766 : isl_int *ineq)
767 : {
768 : int h;
769 :
770 0 : h = hash_index_ineq(ci, &ineq);
771 0 : if (!ci->index[h])
772 0 : return isl_bool_false;
773 0 : return isl_int_ge(ineq[0], (*ci->index[h])[0]);
774 : }
775 :
776 : /* If we can eliminate more than one div, then we need to make
777 : * sure we do it from last div to first div, in order not to
778 : * change the position of the other divs that still need to
779 : * be removed.
780 : */
781 16214254343 : static __isl_give isl_basic_map *remove_duplicate_divs(
782 : __isl_take isl_basic_map *bmap, int *progress)
783 : {
784 : unsigned int size;
785 : int *index;
786 : int *elim_for;
787 : int k, l, h;
788 : int bits;
789 : struct isl_blk eq;
790 : unsigned total_var;
791 : unsigned total;
792 : struct isl_ctx *ctx;
793 :
794 16214254343 : bmap = isl_basic_map_order_divs(bmap);
795 16214254343 : if (!bmap || bmap->n_div <= 1)
796 16214254343 : return bmap;
797 :
798 0 : total_var = isl_space_dim(bmap->dim, isl_dim_all);
799 0 : total = total_var + bmap->n_div;
800 :
801 0 : ctx = bmap->ctx;
802 0 : for (k = bmap->n_div - 1; k >= 0; --k)
803 0 : if (!isl_int_is_zero(bmap->div[k][0]))
804 : break;
805 0 : if (k <= 0)
806 0 : return bmap;
807 :
808 0 : size = round_up(4 * bmap->n_div / 3 - 1);
809 0 : if (size == 0)
810 0 : return bmap;
811 0 : elim_for = isl_calloc_array(ctx, int, bmap->n_div);
812 0 : bits = ffs(size) - 1;
813 0 : index = isl_calloc_array(ctx, int, size);
814 0 : if (!elim_for || !index)
815 : goto out;
816 0 : eq = isl_blk_alloc(ctx, 1+total);
817 0 : if (isl_blk_is_error(eq))
818 0 : goto out;
819 :
820 0 : isl_seq_clr(eq.data, 1+total);
821 0 : index[isl_seq_get_hash_bits(bmap->div[k], 2+total, bits)] = k + 1;
822 0 : for (--k; k >= 0; --k) {
823 : uint32_t hash;
824 :
825 0 : if (isl_int_is_zero(bmap->div[k][0]))
826 0 : continue;
827 :
828 0 : hash = isl_seq_get_hash_bits(bmap->div[k], 2+total, bits);
829 0 : for (h = hash; index[h]; h = (h+1) % size)
830 0 : if (isl_seq_eq(bmap->div[k],
831 0 : bmap->div[index[h]-1], 2+total))
832 0 : break;
833 0 : if (index[h]) {
834 0 : *progress = 1;
835 0 : l = index[h] - 1;
836 0 : elim_for[l] = k + 1;
837 : }
838 0 : index[h] = k+1;
839 : }
840 0 : for (l = bmap->n_div - 1; l >= 0; --l) {
841 0 : if (!elim_for[l])
842 0 : continue;
843 0 : k = elim_for[l] - 1;
844 0 : isl_int_set_si(eq.data[1+total_var+k], -1);
845 0 : isl_int_set_si(eq.data[1+total_var+l], 1);
846 0 : bmap = eliminate_div(bmap, eq.data, l, 1);
847 0 : if (!bmap)
848 0 : break;
849 0 : isl_int_set_si(eq.data[1+total_var+k], 0);
850 0 : isl_int_set_si(eq.data[1+total_var+l], 0);
851 : }
852 :
853 0 : isl_blk_free(ctx, eq);
854 : out:
855 0 : free(index);
856 0 : free(elim_for);
857 0 : return bmap;
858 : }
859 :
860 0 : static int n_pure_div_eq(struct isl_basic_map *bmap)
861 : {
862 : int i, j;
863 : unsigned total;
864 :
865 0 : total = isl_space_dim(bmap->dim, isl_dim_all);
866 0 : for (i = 0, j = bmap->n_div-1; i < bmap->n_eq; ++i) {
867 0 : while (j >= 0 && isl_int_is_zero(bmap->eq[i][1 + total + j]))
868 0 : --j;
869 0 : if (j < 0)
870 0 : break;
871 0 : if (isl_seq_first_non_zero(bmap->eq[i] + 1 + total, j) != -1)
872 0 : return 0;
873 : }
874 0 : return i;
875 : }
876 :
877 : /* Normalize divs that appear in equalities.
878 : *
879 : * In particular, we assume that bmap contains some equalities
880 : * of the form
881 : *
882 : * a x = m * e_i
883 : *
884 : * and we want to replace the set of e_i by a minimal set and
885 : * such that the new e_i have a canonical representation in terms
886 : * of the vector x.
887 : * If any of the equalities involves more than one divs, then
888 : * we currently simply bail out.
889 : *
890 : * Let us first additionally assume that all equalities involve
891 : * a div. The equalities then express modulo constraints on the
892 : * remaining variables and we can use "parameter compression"
893 : * to find a minimal set of constraints. The result is a transformation
894 : *
895 : * x = T(x') = x_0 + G x'
896 : *
897 : * with G a lower-triangular matrix with all elements below the diagonal
898 : * non-negative and smaller than the diagonal element on the same row.
899 : * We first normalize x_0 by making the same property hold in the affine
900 : * T matrix.
901 : * The rows i of G with a 1 on the diagonal do not impose any modulo
902 : * constraint and simply express x_i = x'_i.
903 : * For each of the remaining rows i, we introduce a div and a corresponding
904 : * equality. In particular
905 : *
906 : * g_ii e_j = x_i - g_i(x')
907 : *
908 : * where each x'_k is replaced either by x_k (if g_kk = 1) or the
909 : * corresponding div (if g_kk != 1).
910 : *
911 : * If there are any equalities not involving any div, then we
912 : * first apply a variable compression on the variables x:
913 : *
914 : * x = C x'' x'' = C_2 x
915 : *
916 : * and perform the above parameter compression on A C instead of on A.
917 : * The resulting compression is then of the form
918 : *
919 : * x'' = T(x') = x_0 + G x'
920 : *
921 : * and in constructing the new divs and the corresponding equalities,
922 : * we have to replace each x'', i.e., the x'_k with (g_kk = 1),
923 : * by the corresponding row from C_2.
924 : */
925 16214254343 : static __isl_give isl_basic_map *normalize_divs(__isl_take isl_basic_map *bmap,
926 : int *progress)
927 : {
928 : int i, j, k;
929 : int total;
930 : int div_eq;
931 : struct isl_mat *B;
932 : struct isl_vec *d;
933 16214254343 : struct isl_mat *T = NULL;
934 16214254343 : struct isl_mat *C = NULL;
935 16214254343 : struct isl_mat *C2 = NULL;
936 : isl_int v;
937 16214254343 : int *pos = NULL;
938 : int dropped, needed;
939 :
940 16214254343 : if (!bmap)
941 0 : return NULL;
942 :
943 16214254343 : if (bmap->n_div == 0)
944 16214254343 : return bmap;
945 :
946 0 : if (bmap->n_eq == 0)
947 0 : return bmap;
948 :
949 0 : if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_NORMALIZED_DIVS))
950 0 : return bmap;
951 :
952 0 : total = isl_space_dim(bmap->dim, isl_dim_all);
953 0 : div_eq = n_pure_div_eq(bmap);
954 0 : if (div_eq == 0)
955 0 : return bmap;
956 :
957 0 : if (div_eq < bmap->n_eq) {
958 0 : B = isl_mat_sub_alloc6(bmap->ctx, bmap->eq, div_eq,
959 0 : bmap->n_eq - div_eq, 0, 1 + total);
960 0 : C = isl_mat_variable_compression(B, &C2);
961 0 : if (!C || !C2)
962 : goto error;
963 0 : if (C->n_col == 0) {
964 0 : bmap = isl_basic_map_set_to_empty(bmap);
965 0 : isl_mat_free(C);
966 0 : isl_mat_free(C2);
967 0 : goto done;
968 : }
969 : }
970 :
971 0 : d = isl_vec_alloc(bmap->ctx, div_eq);
972 0 : if (!d)
973 0 : goto error;
974 0 : for (i = 0, j = bmap->n_div-1; i < div_eq; ++i) {
975 0 : while (j >= 0 && isl_int_is_zero(bmap->eq[i][1 + total + j]))
976 0 : --j;
977 0 : isl_int_set(d->block.data[i], bmap->eq[i][1 + total + j]);
978 : }
979 0 : B = isl_mat_sub_alloc6(bmap->ctx, bmap->eq, 0, div_eq, 0, 1 + total);
980 :
981 0 : if (C) {
982 0 : B = isl_mat_product(B, C);
983 0 : C = NULL;
984 : }
985 :
986 0 : T = isl_mat_parameter_compression(B, d);
987 0 : if (!T)
988 0 : goto error;
989 0 : if (T->n_col == 0) {
990 0 : bmap = isl_basic_map_set_to_empty(bmap);
991 0 : isl_mat_free(C2);
992 0 : isl_mat_free(T);
993 0 : goto done;
994 : }
995 0 : isl_int_init(v);
996 0 : for (i = 0; i < T->n_row - 1; ++i) {
997 0 : isl_int_fdiv_q(v, T->row[1 + i][0], T->row[1 + i][1 + i]);
998 0 : if (isl_int_is_zero(v))
999 0 : continue;
1000 0 : isl_mat_col_submul(T, 0, v, 1 + i);
1001 : }
1002 0 : isl_int_clear(v);
1003 0 : pos = isl_alloc_array(bmap->ctx, int, T->n_row);
1004 0 : if (!pos)
1005 0 : goto error;
1006 : /* We have to be careful because dropping equalities may reorder them */
1007 0 : dropped = 0;
1008 0 : for (j = bmap->n_div - 1; j >= 0; --j) {
1009 0 : for (i = 0; i < bmap->n_eq; ++i)
1010 0 : if (!isl_int_is_zero(bmap->eq[i][1 + total + j]))
1011 : break;
1012 0 : if (i < bmap->n_eq) {
1013 0 : bmap = isl_basic_map_drop_div(bmap, j);
1014 0 : isl_basic_map_drop_equality(bmap, i);
1015 0 : ++dropped;
1016 : }
1017 : }
1018 0 : pos[0] = 0;
1019 0 : needed = 0;
1020 0 : for (i = 1; i < T->n_row; ++i) {
1021 0 : if (isl_int_is_one(T->row[i][i]))
1022 0 : pos[i] = i;
1023 : else
1024 0 : needed++;
1025 : }
1026 0 : if (needed > dropped) {
1027 0 : bmap = isl_basic_map_extend_space(bmap, isl_space_copy(bmap->dim),
1028 : needed, needed, 0);
1029 0 : if (!bmap)
1030 0 : goto error;
1031 : }
1032 0 : for (i = 1; i < T->n_row; ++i) {
1033 0 : if (isl_int_is_one(T->row[i][i]))
1034 0 : continue;
1035 0 : k = isl_basic_map_alloc_div(bmap);
1036 0 : pos[i] = 1 + total + k;
1037 0 : isl_seq_clr(bmap->div[k] + 1, 1 + total + bmap->n_div);
1038 0 : isl_int_set(bmap->div[k][0], T->row[i][i]);
1039 0 : if (C2)
1040 0 : isl_seq_cpy(bmap->div[k] + 1, C2->row[i], 1 + total);
1041 : else
1042 0 : isl_int_set_si(bmap->div[k][1 + i], 1);
1043 0 : for (j = 0; j < i; ++j) {
1044 0 : if (isl_int_is_zero(T->row[i][j]))
1045 0 : continue;
1046 0 : if (pos[j] < T->n_row && C2)
1047 0 : isl_seq_submul(bmap->div[k] + 1, T->row[i][j],
1048 0 : C2->row[pos[j]], 1 + total);
1049 : else
1050 0 : isl_int_neg(bmap->div[k][1 + pos[j]],
1051 : T->row[i][j]);
1052 : }
1053 0 : j = isl_basic_map_alloc_equality(bmap);
1054 0 : isl_seq_neg(bmap->eq[j], bmap->div[k]+1, 1+total+bmap->n_div);
1055 0 : isl_int_set(bmap->eq[j][pos[i]], bmap->div[k][0]);
1056 : }
1057 0 : free(pos);
1058 0 : isl_mat_free(C2);
1059 0 : isl_mat_free(T);
1060 :
1061 0 : if (progress)
1062 0 : *progress = 1;
1063 : done:
1064 0 : ISL_F_SET(bmap, ISL_BASIC_MAP_NORMALIZED_DIVS);
1065 :
1066 0 : return bmap;
1067 : error:
1068 0 : free(pos);
1069 0 : isl_mat_free(C);
1070 0 : isl_mat_free(C2);
1071 0 : isl_mat_free(T);
1072 0 : return bmap;
1073 : }
1074 :
1075 0 : static __isl_give isl_basic_map *set_div_from_lower_bound(
1076 : __isl_take isl_basic_map *bmap, int div, int ineq)
1077 : {
1078 0 : unsigned total = 1 + isl_space_dim(bmap->dim, isl_dim_all);
1079 :
1080 0 : isl_seq_neg(bmap->div[div] + 1, bmap->ineq[ineq], total + bmap->n_div);
1081 0 : isl_int_set(bmap->div[div][0], bmap->ineq[ineq][total + div]);
1082 0 : isl_int_add(bmap->div[div][1], bmap->div[div][1], bmap->div[div][0]);
1083 0 : isl_int_sub_ui(bmap->div[div][1], bmap->div[div][1], 1);
1084 0 : isl_int_set_si(bmap->div[div][1 + total + div], 0);
1085 :
1086 0 : return bmap;
1087 : }
1088 :
1089 : /* Check whether it is ok to define a div based on an inequality.
1090 : * To avoid the introduction of circular definitions of divs, we
1091 : * do not allow such a definition if the resulting expression would refer to
1092 : * any other undefined divs or if any known div is defined in
1093 : * terms of the unknown div.
1094 : */
1095 0 : static isl_bool ok_to_set_div_from_bound(__isl_keep isl_basic_map *bmap,
1096 : int div, int ineq)
1097 : {
1098 : int j;
1099 0 : unsigned total = 1 + isl_space_dim(bmap->dim, isl_dim_all);
1100 :
1101 : /* Not defined in terms of unknown divs */
1102 0 : for (j = 0; j < bmap->n_div; ++j) {
1103 0 : if (div == j)
1104 0 : continue;
1105 0 : if (isl_int_is_zero(bmap->ineq[ineq][total + j]))
1106 0 : continue;
1107 0 : if (isl_int_is_zero(bmap->div[j][0]))
1108 0 : return isl_bool_false;
1109 : }
1110 :
1111 : /* No other div defined in terms of this one => avoid loops */
1112 0 : for (j = 0; j < bmap->n_div; ++j) {
1113 0 : if (div == j)
1114 0 : continue;
1115 0 : if (isl_int_is_zero(bmap->div[j][0]))
1116 0 : continue;
1117 0 : if (!isl_int_is_zero(bmap->div[j][1 + total + div]))
1118 0 : return isl_bool_false;
1119 : }
1120 :
1121 0 : return isl_bool_true;
1122 : }
1123 :
1124 : /* Would an expression for div "div" based on inequality "ineq" of "bmap"
1125 : * be a better expression than the current one?
1126 : *
1127 : * If we do not have any expression yet, then any expression would be better.
1128 : * Otherwise we check if the last variable involved in the inequality
1129 : * (disregarding the div that it would define) is in an earlier position
1130 : * than the last variable involved in the current div expression.
1131 : */
1132 0 : static isl_bool better_div_constraint(__isl_keep isl_basic_map *bmap,
1133 : int div, int ineq)
1134 : {
1135 0 : unsigned total = 1 + isl_space_dim(bmap->dim, isl_dim_all);
1136 : int last_div;
1137 : int last_ineq;
1138 :
1139 0 : if (isl_int_is_zero(bmap->div[div][0]))
1140 0 : return isl_bool_true;
1141 :
1142 0 : if (isl_seq_last_non_zero(bmap->ineq[ineq] + total + div + 1,
1143 0 : bmap->n_div - (div + 1)) >= 0)
1144 0 : return isl_bool_false;
1145 :
1146 0 : last_ineq = isl_seq_last_non_zero(bmap->ineq[ineq], total + div);
1147 0 : last_div = isl_seq_last_non_zero(bmap->div[div] + 1,
1148 0 : total + bmap->n_div);
1149 :
1150 0 : return last_ineq < last_div;
1151 : }
1152 :
1153 : /* Given two constraints "k" and "l" that are opposite to each other,
1154 : * except for the constant term, check if we can use them
1155 : * to obtain an expression for one of the hitherto unknown divs or
1156 : * a "better" expression for a div for which we already have an expression.
1157 : * "sum" is the sum of the constant terms of the constraints.
1158 : * If this sum is strictly smaller than the coefficient of one
1159 : * of the divs, then this pair can be used define the div.
1160 : * To avoid the introduction of circular definitions of divs, we
1161 : * do not use the pair if the resulting expression would refer to
1162 : * any other undefined divs or if any known div is defined in
1163 : * terms of the unknown div.
1164 : */
1165 1012353148 : static __isl_give isl_basic_map *check_for_div_constraints(
1166 : __isl_take isl_basic_map *bmap, int k, int l, isl_int sum,
1167 : int *progress)
1168 : {
1169 : int i;
1170 1012353148 : unsigned total = 1 + isl_space_dim(bmap->dim, isl_dim_all);
1171 :
1172 1012353148 : for (i = 0; i < bmap->n_div; ++i) {
1173 : isl_bool set_div;
1174 :
1175 0 : if (isl_int_is_zero(bmap->ineq[k][total + i]))
1176 0 : continue;
1177 0 : if (isl_int_abs_ge(sum, bmap->ineq[k][total + i]))
1178 0 : continue;
1179 0 : set_div = better_div_constraint(bmap, i, k);
1180 0 : if (set_div >= 0 && set_div)
1181 0 : set_div = ok_to_set_div_from_bound(bmap, i, k);
1182 0 : if (set_div < 0)
1183 0 : return isl_basic_map_free(bmap);
1184 0 : if (!set_div)
1185 0 : break;
1186 0 : if (isl_int_is_pos(bmap->ineq[k][total + i]))
1187 0 : bmap = set_div_from_lower_bound(bmap, i, k);
1188 : else
1189 0 : bmap = set_div_from_lower_bound(bmap, i, l);
1190 0 : if (progress)
1191 0 : *progress = 1;
1192 0 : break;
1193 : }
1194 1012353148 : return bmap;
1195 : }
1196 :
1197 16214282001 : __isl_give isl_basic_map *isl_basic_map_remove_duplicate_constraints(
1198 : __isl_take isl_basic_map *bmap, int *progress, int detect_divs)
1199 : {
1200 : struct isl_constraint_index ci;
1201 : int k, l, h;
1202 16214282001 : unsigned total = isl_basic_map_total_dim(bmap);
1203 : isl_int sum;
1204 :
1205 16214282001 : if (!bmap || bmap->n_ineq <= 1)
1206 14221862380 : return bmap;
1207 :
1208 1992419621 : if (create_constraint_index(&ci, bmap) < 0)
1209 0 : return bmap;
1210 :
1211 1992419621 : h = isl_seq_get_hash_bits(bmap->ineq[0] + 1, total, ci.bits);
1212 1992419621 : ci.index[h] = &bmap->ineq[0];
1213 7456164042 : for (k = 1; k < bmap->n_ineq; ++k) {
1214 5463744421 : h = hash_index(&ci, bmap, k);
1215 5463744421 : if (!ci.index[h]) {
1216 5258891874 : ci.index[h] = &bmap->ineq[k];
1217 5258891874 : continue;
1218 : }
1219 204852547 : if (progress)
1220 204852547 : *progress = 1;
1221 204852547 : l = ci.index[h] - &bmap->ineq[0];
1222 204852547 : if (isl_int_lt(bmap->ineq[k][0], bmap->ineq[l][0]))
1223 30940916 : swap_inequality(bmap, k, l);
1224 204852547 : isl_basic_map_drop_inequality(bmap, k);
1225 204852547 : --k;
1226 : }
1227 1992419621 : isl_int_init(sum);
1228 6550137802 : for (k = 0; k < bmap->n_ineq-1; ++k) {
1229 4640001866 : isl_seq_neg(bmap->ineq[k]+1, bmap->ineq[k]+1, total);
1230 4640001866 : h = hash_index(&ci, bmap, k);
1231 4640001866 : isl_seq_neg(bmap->ineq[k]+1, bmap->ineq[k]+1, total);
1232 4640001866 : if (!ci.index[h])
1233 3545268277 : continue;
1234 1094733589 : l = ci.index[h] - &bmap->ineq[0];
1235 1094733589 : isl_int_add(sum, bmap->ineq[k][0], bmap->ineq[l][0]);
1236 1094733589 : if (isl_int_is_pos(sum)) {
1237 1012449904 : if (detect_divs)
1238 1012353148 : bmap = check_for_div_constraints(bmap, k, l,
1239 : sum, progress);
1240 1012449904 : continue;
1241 : }
1242 82283685 : if (isl_int_is_zero(sum)) {
1243 : /* We need to break out of the loop after these
1244 : * changes since the contents of the hash
1245 : * will no longer be valid.
1246 : * Plus, we probably we want to regauss first.
1247 : */
1248 37971387 : if (progress)
1249 37971387 : *progress = 1;
1250 37971387 : isl_basic_map_drop_inequality(bmap, l);
1251 37971387 : isl_basic_map_inequality_to_equality(bmap, k);
1252 : } else
1253 44312298 : bmap = isl_basic_map_set_to_empty(bmap);
1254 82283685 : break;
1255 : }
1256 1992419621 : isl_int_clear(sum);
1257 :
1258 1992419621 : constraint_index_free(&ci);
1259 1992419621 : return bmap;
1260 : }
1261 :
1262 : /* Detect all pairs of inequalities that form an equality.
1263 : *
1264 : * isl_basic_map_remove_duplicate_constraints detects at most one such pair.
1265 : * Call it repeatedly while it is making progress.
1266 : */
1267 27658 : __isl_give isl_basic_map *isl_basic_map_detect_inequality_pairs(
1268 : __isl_take isl_basic_map *bmap, int *progress)
1269 : {
1270 : int duplicate;
1271 :
1272 : do {
1273 27658 : duplicate = 0;
1274 27658 : bmap = isl_basic_map_remove_duplicate_constraints(bmap,
1275 : &duplicate, 0);
1276 27658 : if (progress && duplicate)
1277 221 : *progress = 1;
1278 27658 : } while (duplicate);
1279 :
1280 16142 : return bmap;
1281 : }
1282 :
1283 : /* Eliminate knowns divs from constraints where they appear with
1284 : * a (positive or negative) unit coefficient.
1285 : *
1286 : * That is, replace
1287 : *
1288 : * floor(e/m) + f >= 0
1289 : *
1290 : * by
1291 : *
1292 : * e + m f >= 0
1293 : *
1294 : * and
1295 : *
1296 : * -floor(e/m) + f >= 0
1297 : *
1298 : * by
1299 : *
1300 : * -e + m f + m - 1 >= 0
1301 : *
1302 : * The first conversion is valid because floor(e/m) >= -f is equivalent
1303 : * to e/m >= -f because -f is an integral expression.
1304 : * The second conversion follows from the fact that
1305 : *
1306 : * -floor(e/m) = ceil(-e/m) = floor((-e + m - 1)/m)
1307 : *
1308 : *
1309 : * Note that one of the div constraints may have been eliminated
1310 : * due to being redundant with respect to the constraint that is
1311 : * being modified by this function. The modified constraint may
1312 : * no longer imply this div constraint, so we add it back to make
1313 : * sure we do not lose any information.
1314 : *
1315 : * We skip integral divs, i.e., those with denominator 1, as we would
1316 : * risk eliminating the div from the div constraints. We do not need
1317 : * to handle those divs here anyway since the div constraints will turn
1318 : * out to form an equality and this equality can then be used to eliminate
1319 : * the div from all constraints.
1320 : */
1321 16214254343 : static __isl_give isl_basic_map *eliminate_unit_divs(
1322 : __isl_take isl_basic_map *bmap, int *progress)
1323 : {
1324 : int i, j;
1325 : isl_ctx *ctx;
1326 : unsigned total;
1327 :
1328 16214254343 : if (!bmap)
1329 0 : return NULL;
1330 :
1331 16214254343 : ctx = isl_basic_map_get_ctx(bmap);
1332 16214254343 : total = 1 + isl_space_dim(bmap->dim, isl_dim_all);
1333 :
1334 16214254343 : for (i = 0; i < bmap->n_div; ++i) {
1335 0 : if (isl_int_is_zero(bmap->div[i][0]))
1336 0 : continue;
1337 0 : if (isl_int_is_one(bmap->div[i][0]))
1338 0 : continue;
1339 0 : for (j = 0; j < bmap->n_ineq; ++j) {
1340 : int s;
1341 :
1342 0 : if (!isl_int_is_one(bmap->ineq[j][total + i]) &&
1343 0 : !isl_int_is_negone(bmap->ineq[j][total + i]))
1344 0 : continue;
1345 :
1346 0 : *progress = 1;
1347 :
1348 0 : s = isl_int_sgn(bmap->ineq[j][total + i]);
1349 0 : isl_int_set_si(bmap->ineq[j][total + i], 0);
1350 0 : if (s < 0)
1351 0 : isl_seq_combine(bmap->ineq[j],
1352 0 : ctx->negone, bmap->div[i] + 1,
1353 0 : bmap->div[i][0], bmap->ineq[j],
1354 0 : total + bmap->n_div);
1355 : else
1356 0 : isl_seq_combine(bmap->ineq[j],
1357 0 : ctx->one, bmap->div[i] + 1,
1358 0 : bmap->div[i][0], bmap->ineq[j],
1359 0 : total + bmap->n_div);
1360 0 : if (s < 0) {
1361 0 : isl_int_add(bmap->ineq[j][0],
1362 : bmap->ineq[j][0], bmap->div[i][0]);
1363 0 : isl_int_sub_ui(bmap->ineq[j][0],
1364 : bmap->ineq[j][0], 1);
1365 : }
1366 :
1367 0 : bmap = isl_basic_map_extend_constraints(bmap, 0, 1);
1368 0 : if (isl_basic_map_add_div_constraint(bmap, i, s) < 0)
1369 0 : return isl_basic_map_free(bmap);
1370 : }
1371 : }
1372 :
1373 16214254343 : return bmap;
1374 : }
1375 :
1376 12231209136 : __isl_give isl_basic_map *isl_basic_map_simplify(__isl_take isl_basic_map *bmap)
1377 : {
1378 12231209136 : int progress = 1;
1379 12231209136 : if (!bmap)
1380 0 : return NULL;
1381 40676672615 : while (progress) {
1382 : isl_bool empty;
1383 :
1384 17041287968 : progress = 0;
1385 17041287968 : empty = isl_basic_map_plain_is_empty(bmap);
1386 17041287968 : if (empty < 0)
1387 0 : return isl_basic_map_free(bmap);
1388 17041287968 : if (empty)
1389 827033625 : break;
1390 16214254343 : bmap = isl_basic_map_normalize_constraints(bmap);
1391 16214254343 : bmap = reduce_div_coefficients(bmap);
1392 16214254343 : bmap = normalize_div_expressions(bmap);
1393 16214254343 : bmap = remove_duplicate_divs(bmap, &progress);
1394 16214254343 : bmap = eliminate_unit_divs(bmap, &progress);
1395 16214254343 : bmap = eliminate_divs_eq(bmap, &progress);
1396 16214254343 : bmap = eliminate_divs_ineq(bmap, &progress);
1397 16214254343 : bmap = isl_basic_map_gauss(bmap, &progress);
1398 : /* requires equalities in normal form */
1399 16214254343 : bmap = normalize_divs(bmap, &progress);
1400 16214254343 : bmap = isl_basic_map_remove_duplicate_constraints(bmap,
1401 : &progress, 1);
1402 16214254343 : if (bmap && progress)
1403 4810078832 : ISL_F_CLR(bmap, ISL_BASIC_MAP_REDUCED_COEFFICIENTS);
1404 : }
1405 12231209136 : return bmap;
1406 : }
1407 :
1408 11421529127 : struct isl_basic_set *isl_basic_set_simplify(struct isl_basic_set *bset)
1409 : {
1410 11421529127 : return bset_from_bmap(isl_basic_map_simplify(bset_to_bmap(bset)));
1411 : }
1412 :
1413 :
1414 0 : isl_bool isl_basic_map_is_div_constraint(__isl_keep isl_basic_map *bmap,
1415 : isl_int *constraint, unsigned div)
1416 : {
1417 : unsigned pos;
1418 :
1419 0 : if (!bmap)
1420 0 : return isl_bool_error;
1421 :
1422 0 : pos = 1 + isl_space_dim(bmap->dim, isl_dim_all) + div;
1423 :
1424 0 : if (isl_int_eq(constraint[pos], bmap->div[div][0])) {
1425 : int neg;
1426 0 : isl_int_sub(bmap->div[div][1],
1427 : bmap->div[div][1], bmap->div[div][0]);
1428 0 : isl_int_add_ui(bmap->div[div][1], bmap->div[div][1], 1);
1429 0 : neg = isl_seq_is_neg(constraint, bmap->div[div]+1, pos);
1430 0 : isl_int_sub_ui(bmap->div[div][1], bmap->div[div][1], 1);
1431 0 : isl_int_add(bmap->div[div][1],
1432 : bmap->div[div][1], bmap->div[div][0]);
1433 0 : if (!neg)
1434 0 : return isl_bool_false;
1435 0 : if (isl_seq_first_non_zero(constraint+pos+1,
1436 0 : bmap->n_div-div-1) != -1)
1437 0 : return isl_bool_false;
1438 0 : } else if (isl_int_abs_eq(constraint[pos], bmap->div[div][0])) {
1439 0 : if (!isl_seq_eq(constraint, bmap->div[div]+1, pos))
1440 0 : return isl_bool_false;
1441 0 : if (isl_seq_first_non_zero(constraint+pos+1,
1442 0 : bmap->n_div-div-1) != -1)
1443 0 : return isl_bool_false;
1444 : } else
1445 0 : return isl_bool_false;
1446 :
1447 0 : return isl_bool_true;
1448 : }
1449 :
1450 0 : isl_bool isl_basic_set_is_div_constraint(__isl_keep isl_basic_set *bset,
1451 : isl_int *constraint, unsigned div)
1452 : {
1453 0 : return isl_basic_map_is_div_constraint(bset, constraint, div);
1454 : }
1455 :
1456 :
1457 : /* If the only constraints a div d=floor(f/m)
1458 : * appears in are its two defining constraints
1459 : *
1460 : * f - m d >=0
1461 : * -(f - (m - 1)) + m d >= 0
1462 : *
1463 : * then it can safely be removed.
1464 : */
1465 0 : static isl_bool div_is_redundant(__isl_keep isl_basic_map *bmap, int div)
1466 : {
1467 : int i;
1468 0 : unsigned pos = 1 + isl_space_dim(bmap->dim, isl_dim_all) + div;
1469 :
1470 0 : for (i = 0; i < bmap->n_eq; ++i)
1471 0 : if (!isl_int_is_zero(bmap->eq[i][pos]))
1472 0 : return isl_bool_false;
1473 :
1474 0 : for (i = 0; i < bmap->n_ineq; ++i) {
1475 : isl_bool red;
1476 :
1477 0 : if (isl_int_is_zero(bmap->ineq[i][pos]))
1478 0 : continue;
1479 0 : red = isl_basic_map_is_div_constraint(bmap, bmap->ineq[i], div);
1480 0 : if (red < 0 || !red)
1481 0 : return red;
1482 : }
1483 :
1484 0 : for (i = 0; i < bmap->n_div; ++i) {
1485 0 : if (isl_int_is_zero(bmap->div[i][0]))
1486 0 : continue;
1487 0 : if (!isl_int_is_zero(bmap->div[i][1+pos]))
1488 0 : return isl_bool_false;
1489 : }
1490 :
1491 0 : return isl_bool_true;
1492 : }
1493 :
1494 : /*
1495 : * Remove divs that don't occur in any of the constraints or other divs.
1496 : * These can arise when dropping constraints from a basic map or
1497 : * when the divs of a basic map have been temporarily aligned
1498 : * with the divs of another basic map.
1499 : */
1500 25838543167 : static __isl_give isl_basic_map *remove_redundant_divs(
1501 : __isl_take isl_basic_map *bmap)
1502 : {
1503 : int i;
1504 :
1505 25838543167 : if (!bmap)
1506 0 : return NULL;
1507 :
1508 25838543167 : for (i = bmap->n_div-1; i >= 0; --i) {
1509 : isl_bool redundant;
1510 :
1511 0 : redundant = div_is_redundant(bmap, i);
1512 0 : if (redundant < 0)
1513 0 : return isl_basic_map_free(bmap);
1514 0 : if (!redundant)
1515 0 : continue;
1516 0 : bmap = isl_basic_map_drop_div(bmap, i);
1517 : }
1518 25838543167 : return bmap;
1519 : }
1520 :
1521 : /* Mark "bmap" as final, without checking for obviously redundant
1522 : * integer divisions. This function should be used when "bmap"
1523 : * is known not to involve any such integer divisions.
1524 : */
1525 25838543167 : __isl_give isl_basic_map *isl_basic_map_mark_final(
1526 : __isl_take isl_basic_map *bmap)
1527 : {
1528 25838543167 : if (!bmap)
1529 0 : return NULL;
1530 25838543167 : ISL_F_SET(bmap, ISL_BASIC_SET_FINAL);
1531 25838543167 : return bmap;
1532 : }
1533 :
1534 : /* Mark "bmap" as final, after removing obviously redundant integer divisions.
1535 : */
1536 25838543167 : __isl_give isl_basic_map *isl_basic_map_finalize(__isl_take isl_basic_map *bmap)
1537 : {
1538 25838543167 : bmap = remove_redundant_divs(bmap);
1539 25838543167 : bmap = isl_basic_map_mark_final(bmap);
1540 25838543167 : return bmap;
1541 : }
1542 :
1543 10625867154 : struct isl_basic_set *isl_basic_set_finalize(struct isl_basic_set *bset)
1544 : {
1545 10625867154 : return bset_from_bmap(isl_basic_map_finalize(bset_to_bmap(bset)));
1546 : }
1547 :
1548 : /* Remove definition of any div that is defined in terms of the given variable.
1549 : * The div itself is not removed. Functions such as
1550 : * eliminate_divs_ineq depend on the other divs remaining in place.
1551 : */
1552 1275257 : static __isl_give isl_basic_map *remove_dependent_vars(
1553 : __isl_take isl_basic_map *bmap, int pos)
1554 : {
1555 : int i;
1556 :
1557 1275257 : if (!bmap)
1558 0 : return NULL;
1559 :
1560 1275257 : for (i = 0; i < bmap->n_div; ++i) {
1561 0 : if (isl_int_is_zero(bmap->div[i][0]))
1562 0 : continue;
1563 0 : if (isl_int_is_zero(bmap->div[i][1+1+pos]))
1564 0 : continue;
1565 0 : bmap = isl_basic_map_mark_div_unknown(bmap, i);
1566 0 : if (!bmap)
1567 0 : return NULL;
1568 : }
1569 1275257 : return bmap;
1570 : }
1571 :
1572 : /* Eliminate the specified variables from the constraints using
1573 : * Fourier-Motzkin. The variables themselves are not removed.
1574 : */
1575 752725 : __isl_give isl_basic_map *isl_basic_map_eliminate_vars(
1576 : __isl_take isl_basic_map *bmap, unsigned pos, unsigned n)
1577 : {
1578 : int d;
1579 : int i, j, k;
1580 : unsigned total;
1581 752725 : int need_gauss = 0;
1582 :
1583 752725 : if (n == 0)
1584 0 : return bmap;
1585 752725 : if (!bmap)
1586 0 : return NULL;
1587 752725 : total = isl_basic_map_total_dim(bmap);
1588 :
1589 752725 : bmap = isl_basic_map_cow(bmap);
1590 2027982 : for (d = pos + n - 1; d >= 0 && d >= pos; --d)
1591 1275257 : bmap = remove_dependent_vars(bmap, d);
1592 752725 : if (!bmap)
1593 0 : return NULL;
1594 :
1595 1505450 : for (d = pos + n - 1;
1596 752725 : d >= 0 && d >= total - bmap->n_div && d >= pos; --d)
1597 0 : isl_seq_clr(bmap->div[d-(total-bmap->n_div)], 2+total);
1598 2027982 : for (d = pos + n - 1; d >= 0 && d >= pos; --d) {
1599 : int n_lower, n_upper;
1600 1275257 : if (!bmap)
1601 0 : return NULL;
1602 2922202 : for (i = 0; i < bmap->n_eq; ++i) {
1603 1461101 : if (isl_int_is_zero(bmap->eq[i][1+d]))
1604 185844 : continue;
1605 1275257 : eliminate_var_using_equality(bmap, d, bmap->eq[i], 0, NULL);
1606 1275257 : isl_basic_map_drop_equality(bmap, i);
1607 1275257 : need_gauss = 1;
1608 1275257 : break;
1609 : }
1610 1275257 : if (i < bmap->n_eq)
1611 384265 : continue;
1612 890992 : n_lower = 0;
1613 890992 : n_upper = 0;
1614 4166777 : for (i = 0; i < bmap->n_ineq; ++i) {
1615 3275785 : if (isl_int_is_pos(bmap->ineq[i][1+d]))
1616 0 : n_lower++;
1617 3275785 : else if (isl_int_is_neg(bmap->ineq[i][1+d]))
1618 0 : n_upper++;
1619 : }
1620 890992 : bmap = isl_basic_map_extend_constraints(bmap,
1621 890992 : 0, n_lower * n_upper);
1622 890992 : if (!bmap)
1623 0 : goto error;
1624 4166777 : for (i = bmap->n_ineq - 1; i >= 0; --i) {
1625 : int last;
1626 3275785 : if (isl_int_is_zero(bmap->ineq[i][1+d]))
1627 3275785 : continue;
1628 0 : last = -1;
1629 0 : for (j = 0; j < i; ++j) {
1630 0 : if (isl_int_is_zero(bmap->ineq[j][1+d]))
1631 0 : continue;
1632 0 : last = j;
1633 0 : if (isl_int_sgn(bmap->ineq[i][1+d]) ==
1634 0 : isl_int_sgn(bmap->ineq[j][1+d]))
1635 0 : continue;
1636 0 : k = isl_basic_map_alloc_inequality(bmap);
1637 0 : if (k < 0)
1638 0 : goto error;
1639 0 : isl_seq_cpy(bmap->ineq[k], bmap->ineq[i],
1640 : 1+total);
1641 0 : isl_seq_elim(bmap->ineq[k], bmap->ineq[j],
1642 0 : 1+d, 1+total, NULL);
1643 : }
1644 0 : isl_basic_map_drop_inequality(bmap, i);
1645 0 : i = last + 1;
1646 : }
1647 890992 : if (n_lower > 0 && n_upper > 0) {
1648 0 : bmap = isl_basic_map_normalize_constraints(bmap);
1649 0 : bmap = isl_basic_map_remove_duplicate_constraints(bmap,
1650 : NULL, 0);
1651 0 : bmap = isl_basic_map_gauss(bmap, NULL);
1652 0 : bmap = isl_basic_map_remove_redundancies(bmap);
1653 0 : need_gauss = 0;
1654 0 : if (!bmap)
1655 0 : goto error;
1656 0 : if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
1657 0 : break;
1658 : }
1659 : }
1660 752725 : ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
1661 752725 : if (need_gauss)
1662 752725 : bmap = isl_basic_map_gauss(bmap, NULL);
1663 752725 : return bmap;
1664 : error:
1665 0 : isl_basic_map_free(bmap);
1666 0 : return NULL;
1667 : }
1668 :
1669 0 : struct isl_basic_set *isl_basic_set_eliminate_vars(
1670 : struct isl_basic_set *bset, unsigned pos, unsigned n)
1671 : {
1672 0 : return bset_from_bmap(isl_basic_map_eliminate_vars(bset_to_bmap(bset),
1673 : pos, n));
1674 : }
1675 :
1676 : /* Eliminate the specified n dimensions starting at first from the
1677 : * constraints, without removing the dimensions from the space.
1678 : * If the set is rational, the dimensions are eliminated using Fourier-Motzkin.
1679 : * Otherwise, they are projected out and the original space is restored.
1680 : */
1681 0 : __isl_give isl_basic_map *isl_basic_map_eliminate(
1682 : __isl_take isl_basic_map *bmap,
1683 : enum isl_dim_type type, unsigned first, unsigned n)
1684 : {
1685 : isl_space *space;
1686 :
1687 0 : if (!bmap)
1688 0 : return NULL;
1689 0 : if (n == 0)
1690 0 : return bmap;
1691 :
1692 0 : if (first + n > isl_basic_map_dim(bmap, type) || first + n < first)
1693 0 : isl_die(bmap->ctx, isl_error_invalid,
1694 : "index out of bounds", goto error);
1695 :
1696 0 : if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL)) {
1697 0 : first += isl_basic_map_offset(bmap, type) - 1;
1698 0 : bmap = isl_basic_map_eliminate_vars(bmap, first, n);
1699 0 : return isl_basic_map_finalize(bmap);
1700 : }
1701 :
1702 0 : space = isl_basic_map_get_space(bmap);
1703 0 : bmap = isl_basic_map_project_out(bmap, type, first, n);
1704 0 : bmap = isl_basic_map_insert_dims(bmap, type, first, n);
1705 0 : bmap = isl_basic_map_reset_space(bmap, space);
1706 0 : return bmap;
1707 : error:
1708 0 : isl_basic_map_free(bmap);
1709 0 : return NULL;
1710 : }
1711 :
1712 0 : __isl_give isl_basic_set *isl_basic_set_eliminate(
1713 : __isl_take isl_basic_set *bset,
1714 : enum isl_dim_type type, unsigned first, unsigned n)
1715 : {
1716 0 : return isl_basic_map_eliminate(bset, type, first, n);
1717 : }
1718 :
1719 : /* Remove all constraints from "bmap" that reference any unknown local
1720 : * variables (directly or indirectly).
1721 : *
1722 : * Dropping all constraints on a local variable will make it redundant,
1723 : * so it will get removed implicitly by
1724 : * isl_basic_map_drop_constraints_involving_dims. Some other local
1725 : * variables may also end up becoming redundant if they only appear
1726 : * in constraints together with the unknown local variable.
1727 : * Therefore, start over after calling
1728 : * isl_basic_map_drop_constraints_involving_dims.
1729 : */
1730 0 : __isl_give isl_basic_map *isl_basic_map_drop_constraint_involving_unknown_divs(
1731 : __isl_take isl_basic_map *bmap)
1732 : {
1733 : isl_bool known;
1734 : int i, n_div, o_div;
1735 :
1736 0 : known = isl_basic_map_divs_known(bmap);
1737 0 : if (known < 0)
1738 0 : return isl_basic_map_free(bmap);
1739 0 : if (known)
1740 0 : return bmap;
1741 :
1742 0 : n_div = isl_basic_map_dim(bmap, isl_dim_div);
1743 0 : o_div = isl_basic_map_offset(bmap, isl_dim_div) - 1;
1744 :
1745 0 : for (i = 0; i < n_div; ++i) {
1746 0 : known = isl_basic_map_div_is_known(bmap, i);
1747 0 : if (known < 0)
1748 0 : return isl_basic_map_free(bmap);
1749 0 : if (known)
1750 0 : continue;
1751 0 : bmap = remove_dependent_vars(bmap, o_div + i);
1752 0 : bmap = isl_basic_map_drop_constraints_involving_dims(bmap,
1753 : isl_dim_div, i, 1);
1754 0 : if (!bmap)
1755 0 : return NULL;
1756 0 : n_div = isl_basic_map_dim(bmap, isl_dim_div);
1757 0 : i = -1;
1758 : }
1759 :
1760 0 : return bmap;
1761 : }
1762 :
1763 : /* Remove all constraints from "map" that reference any unknown local
1764 : * variables (directly or indirectly).
1765 : *
1766 : * Since constraints may get dropped from the basic maps,
1767 : * they may no longer be disjoint from each other.
1768 : */
1769 0 : __isl_give isl_map *isl_map_drop_constraint_involving_unknown_divs(
1770 : __isl_take isl_map *map)
1771 : {
1772 : int i;
1773 : isl_bool known;
1774 :
1775 0 : known = isl_map_divs_known(map);
1776 0 : if (known < 0)
1777 0 : return isl_map_free(map);
1778 0 : if (known)
1779 0 : return map;
1780 :
1781 0 : map = isl_map_cow(map);
1782 0 : if (!map)
1783 0 : return NULL;
1784 :
1785 0 : for (i = 0; i < map->n; ++i) {
1786 0 : map->p[i] =
1787 0 : isl_basic_map_drop_constraint_involving_unknown_divs(
1788 0 : map->p[i]);
1789 0 : if (!map->p[i])
1790 0 : return isl_map_free(map);
1791 : }
1792 :
1793 0 : if (map->n > 1)
1794 0 : ISL_F_CLR(map, ISL_MAP_DISJOINT);
1795 :
1796 0 : return map;
1797 : }
1798 :
1799 : /* Don't assume equalities are in order, because align_divs
1800 : * may have changed the order of the divs.
1801 : */
1802 6822802 : static void compute_elimination_index(__isl_keep isl_basic_map *bmap, int *elim)
1803 : {
1804 : int d, i;
1805 : unsigned total;
1806 :
1807 6822802 : total = isl_space_dim(bmap->dim, isl_dim_all);
1808 51047165 : for (d = 0; d < total; ++d)
1809 44224363 : elim[d] = -1;
1810 19107308 : for (i = 0; i < bmap->n_eq; ++i) {
1811 90438436 : for (d = total - 1; d >= 0; --d) {
1812 45219218 : if (isl_int_is_zero(bmap->eq[i][1+d]))
1813 32934712 : continue;
1814 12284506 : elim[d] = i;
1815 12284506 : break;
1816 : }
1817 : }
1818 6822802 : }
1819 :
1820 0 : static void set_compute_elimination_index(__isl_keep isl_basic_set *bset,
1821 : int *elim)
1822 : {
1823 0 : compute_elimination_index(bset_to_bmap(bset), elim);
1824 0 : }
1825 :
1826 51711391 : static int reduced_using_equalities(isl_int *dst, isl_int *src,
1827 : __isl_keep isl_basic_map *bmap, int *elim)
1828 : {
1829 : int d;
1830 51711391 : int copied = 0;
1831 : unsigned total;
1832 :
1833 51711391 : total = isl_space_dim(bmap->dim, isl_dim_all);
1834 393221630 : for (d = total - 1; d >= 0; --d) {
1835 341510239 : if (isl_int_is_zero(src[1+d]))
1836 204316367 : continue;
1837 137193872 : if (elim[d] == -1)
1838 123152337 : continue;
1839 14041535 : if (!copied) {
1840 11620235 : isl_seq_cpy(dst, src, 1 + total);
1841 11620235 : copied = 1;
1842 : }
1843 14041535 : isl_seq_elim(dst, bmap->eq[elim[d]], 1 + d, 1 + total, NULL);
1844 : }
1845 51711391 : return copied;
1846 : }
1847 :
1848 0 : static int set_reduced_using_equalities(isl_int *dst, isl_int *src,
1849 : __isl_keep isl_basic_set *bset, int *elim)
1850 : {
1851 0 : return reduced_using_equalities(dst, src,
1852 : bset_to_bmap(bset), elim);
1853 : }
1854 :
1855 0 : static __isl_give isl_basic_set *isl_basic_set_reduce_using_equalities(
1856 : __isl_take isl_basic_set *bset, __isl_take isl_basic_set *context)
1857 : {
1858 : int i;
1859 : int *elim;
1860 :
1861 0 : if (!bset || !context)
1862 : goto error;
1863 :
1864 0 : if (context->n_eq == 0) {
1865 0 : isl_basic_set_free(context);
1866 0 : return bset;
1867 : }
1868 :
1869 0 : bset = isl_basic_set_cow(bset);
1870 0 : if (!bset)
1871 0 : goto error;
1872 :
1873 0 : elim = isl_alloc_array(bset->ctx, int, isl_basic_set_n_dim(bset));
1874 0 : if (!elim)
1875 0 : goto error;
1876 0 : set_compute_elimination_index(context, elim);
1877 0 : for (i = 0; i < bset->n_eq; ++i)
1878 0 : set_reduced_using_equalities(bset->eq[i], bset->eq[i],
1879 : context, elim);
1880 0 : for (i = 0; i < bset->n_ineq; ++i)
1881 0 : set_reduced_using_equalities(bset->ineq[i], bset->ineq[i],
1882 : context, elim);
1883 0 : isl_basic_set_free(context);
1884 0 : free(elim);
1885 0 : bset = isl_basic_set_simplify(bset);
1886 0 : bset = isl_basic_set_finalize(bset);
1887 0 : return bset;
1888 : error:
1889 0 : isl_basic_set_free(bset);
1890 0 : isl_basic_set_free(context);
1891 0 : return NULL;
1892 : }
1893 :
1894 : /* For each inequality in "ineq" that is a shifted (more relaxed)
1895 : * copy of an inequality in "context", mark the corresponding entry
1896 : * in "row" with -1.
1897 : * If an inequality only has a non-negative constant term, then
1898 : * mark it as well.
1899 : */
1900 0 : static isl_stat mark_shifted_constraints(__isl_keep isl_mat *ineq,
1901 : __isl_keep isl_basic_set *context, int *row)
1902 : {
1903 : struct isl_constraint_index ci;
1904 : int n_ineq;
1905 : unsigned total;
1906 : int k;
1907 :
1908 0 : if (!ineq || !context)
1909 0 : return isl_stat_error;
1910 0 : if (context->n_ineq == 0)
1911 0 : return isl_stat_ok;
1912 0 : if (setup_constraint_index(&ci, context) < 0)
1913 0 : return isl_stat_error;
1914 :
1915 0 : n_ineq = isl_mat_rows(ineq);
1916 0 : total = isl_mat_cols(ineq) - 1;
1917 0 : for (k = 0; k < n_ineq; ++k) {
1918 : int l;
1919 : isl_bool redundant;
1920 :
1921 0 : l = isl_seq_first_non_zero(ineq->row[k] + 1, total);
1922 0 : if (l < 0 && isl_int_is_nonneg(ineq->row[k][0])) {
1923 0 : row[k] = -1;
1924 0 : continue;
1925 : }
1926 0 : redundant = constraint_index_is_redundant(&ci, ineq->row[k]);
1927 0 : if (redundant < 0)
1928 0 : goto error;
1929 0 : if (!redundant)
1930 0 : continue;
1931 0 : row[k] = -1;
1932 : }
1933 0 : constraint_index_free(&ci);
1934 0 : return isl_stat_ok;
1935 : error:
1936 0 : constraint_index_free(&ci);
1937 0 : return isl_stat_error;
1938 : }
1939 :
1940 0 : static __isl_give isl_basic_set *remove_shifted_constraints(
1941 : __isl_take isl_basic_set *bset, __isl_keep isl_basic_set *context)
1942 : {
1943 : struct isl_constraint_index ci;
1944 : int k;
1945 :
1946 0 : if (!bset || !context)
1947 0 : return bset;
1948 :
1949 0 : if (context->n_ineq == 0)
1950 0 : return bset;
1951 0 : if (setup_constraint_index(&ci, context) < 0)
1952 0 : return bset;
1953 :
1954 0 : for (k = 0; k < bset->n_ineq; ++k) {
1955 : isl_bool redundant;
1956 :
1957 0 : redundant = constraint_index_is_redundant(&ci, bset->ineq[k]);
1958 0 : if (redundant < 0)
1959 0 : goto error;
1960 0 : if (!redundant)
1961 0 : continue;
1962 0 : bset = isl_basic_set_cow(bset);
1963 0 : if (!bset)
1964 0 : goto error;
1965 0 : isl_basic_set_drop_inequality(bset, k);
1966 0 : --k;
1967 : }
1968 0 : constraint_index_free(&ci);
1969 0 : return bset;
1970 : error:
1971 0 : constraint_index_free(&ci);
1972 0 : return bset;
1973 : }
1974 :
1975 : /* Remove constraints from "bmap" that are identical to constraints
1976 : * in "context" or that are more relaxed (greater constant term).
1977 : *
1978 : * We perform the test for shifted copies on the pure constraints
1979 : * in remove_shifted_constraints.
1980 : */
1981 0 : static __isl_give isl_basic_map *isl_basic_map_remove_shifted_constraints(
1982 : __isl_take isl_basic_map *bmap, __isl_take isl_basic_map *context)
1983 : {
1984 : isl_basic_set *bset, *bset_context;
1985 :
1986 0 : if (!bmap || !context)
1987 : goto error;
1988 :
1989 0 : if (bmap->n_ineq == 0 || context->n_ineq == 0) {
1990 0 : isl_basic_map_free(context);
1991 0 : return bmap;
1992 : }
1993 :
1994 0 : context = isl_basic_map_align_divs(context, bmap);
1995 0 : bmap = isl_basic_map_align_divs(bmap, context);
1996 :
1997 0 : bset = isl_basic_map_underlying_set(isl_basic_map_copy(bmap));
1998 0 : bset_context = isl_basic_map_underlying_set(context);
1999 0 : bset = remove_shifted_constraints(bset, bset_context);
2000 0 : isl_basic_set_free(bset_context);
2001 :
2002 0 : bmap = isl_basic_map_overlying_set(bset, bmap);
2003 :
2004 0 : return bmap;
2005 : error:
2006 0 : isl_basic_map_free(bmap);
2007 0 : isl_basic_map_free(context);
2008 0 : return NULL;
2009 : }
2010 :
2011 : /* Does the (linear part of a) constraint "c" involve any of the "len"
2012 : * "relevant" dimensions?
2013 : */
2014 0 : static int is_related(isl_int *c, int len, int *relevant)
2015 : {
2016 : int i;
2017 :
2018 0 : for (i = 0; i < len; ++i) {
2019 0 : if (!relevant[i])
2020 0 : continue;
2021 0 : if (!isl_int_is_zero(c[i]))
2022 0 : return 1;
2023 : }
2024 :
2025 0 : return 0;
2026 : }
2027 :
2028 : /* Drop constraints from "bmap" that do not involve any of
2029 : * the dimensions marked "relevant".
2030 : */
2031 0 : static __isl_give isl_basic_map *drop_unrelated_constraints(
2032 : __isl_take isl_basic_map *bmap, int *relevant)
2033 : {
2034 : int i, dim;
2035 :
2036 0 : dim = isl_basic_map_dim(bmap, isl_dim_all);
2037 0 : for (i = 0; i < dim; ++i)
2038 0 : if (!relevant[i])
2039 0 : break;
2040 0 : if (i >= dim)
2041 0 : return bmap;
2042 :
2043 0 : for (i = bmap->n_eq - 1; i >= 0; --i)
2044 0 : if (!is_related(bmap->eq[i] + 1, dim, relevant)) {
2045 0 : bmap = isl_basic_map_cow(bmap);
2046 0 : if (isl_basic_map_drop_equality(bmap, i) < 0)
2047 0 : return isl_basic_map_free(bmap);
2048 : }
2049 :
2050 0 : for (i = bmap->n_ineq - 1; i >= 0; --i)
2051 0 : if (!is_related(bmap->ineq[i] + 1, dim, relevant)) {
2052 0 : bmap = isl_basic_map_cow(bmap);
2053 0 : if (isl_basic_map_drop_inequality(bmap, i) < 0)
2054 0 : return isl_basic_map_free(bmap);
2055 : }
2056 :
2057 0 : return bmap;
2058 : }
2059 :
2060 : /* Update the groups in "group" based on the (linear part of a) constraint "c".
2061 : *
2062 : * In particular, for any variable involved in the constraint,
2063 : * find the actual group id from before and replace the group
2064 : * of the corresponding variable by the minimal group of all
2065 : * the variables involved in the constraint considered so far
2066 : * (if this minimum is smaller) or replace the minimum by this group
2067 : * (if the minimum is larger).
2068 : *
2069 : * At the end, all the variables in "c" will (indirectly) point
2070 : * to the minimal of the groups that they referred to originally.
2071 : */
2072 0 : static void update_groups(int dim, int *group, isl_int *c)
2073 : {
2074 : int j;
2075 0 : int min = dim;
2076 :
2077 0 : for (j = 0; j < dim; ++j) {
2078 0 : if (isl_int_is_zero(c[j]))
2079 0 : continue;
2080 0 : while (group[j] >= 0 && group[group[j]] != group[j])
2081 0 : group[j] = group[group[j]];
2082 0 : if (group[j] == min)
2083 0 : continue;
2084 0 : if (group[j] < min) {
2085 0 : if (min >= 0 && min < dim)
2086 0 : group[min] = group[j];
2087 0 : min = group[j];
2088 : } else
2089 0 : group[group[j]] = min;
2090 : }
2091 0 : }
2092 :
2093 : /* Allocate an array of groups of variables, one for each variable
2094 : * in "context", initialized to zero.
2095 : */
2096 0 : static int *alloc_groups(__isl_keep isl_basic_set *context)
2097 : {
2098 : isl_ctx *ctx;
2099 : int dim;
2100 :
2101 0 : dim = isl_basic_set_dim(context, isl_dim_set);
2102 0 : ctx = isl_basic_set_get_ctx(context);
2103 0 : return isl_calloc_array(ctx, int, dim);
2104 : }
2105 :
2106 : /* Drop constraints from "bmap" that only involve variables that are
2107 : * not related to any of the variables marked with a "-1" in "group".
2108 : *
2109 : * We construct groups of variables that collect variables that
2110 : * (indirectly) appear in some common constraint of "bmap".
2111 : * Each group is identified by the first variable in the group,
2112 : * except for the special group of variables that was already identified
2113 : * in the input as -1 (or are related to those variables).
2114 : * If group[i] is equal to i (or -1), then the group of i is i (or -1),
2115 : * otherwise the group of i is the group of group[i].
2116 : *
2117 : * We first initialize groups for the remaining variables.
2118 : * Then we iterate over the constraints of "bmap" and update the
2119 : * group of the variables in the constraint by the smallest group.
2120 : * Finally, we resolve indirect references to groups by running over
2121 : * the variables.
2122 : *
2123 : * After computing the groups, we drop constraints that do not involve
2124 : * any variables in the -1 group.
2125 : */
2126 0 : __isl_give isl_basic_map *isl_basic_map_drop_unrelated_constraints(
2127 : __isl_take isl_basic_map *bmap, __isl_take int *group)
2128 : {
2129 : int dim;
2130 : int i;
2131 : int last;
2132 :
2133 0 : if (!bmap)
2134 0 : return NULL;
2135 :
2136 0 : dim = isl_basic_map_dim(bmap, isl_dim_all);
2137 :
2138 0 : last = -1;
2139 0 : for (i = 0; i < dim; ++i)
2140 0 : if (group[i] >= 0)
2141 0 : last = group[i] = i;
2142 0 : if (last < 0) {
2143 0 : free(group);
2144 0 : return bmap;
2145 : }
2146 :
2147 0 : for (i = 0; i < bmap->n_eq; ++i)
2148 0 : update_groups(dim, group, bmap->eq[i] + 1);
2149 0 : for (i = 0; i < bmap->n_ineq; ++i)
2150 0 : update_groups(dim, group, bmap->ineq[i] + 1);
2151 :
2152 0 : for (i = 0; i < dim; ++i)
2153 0 : if (group[i] >= 0)
2154 0 : group[i] = group[group[i]];
2155 :
2156 0 : for (i = 0; i < dim; ++i)
2157 0 : group[i] = group[i] == -1;
2158 :
2159 0 : bmap = drop_unrelated_constraints(bmap, group);
2160 :
2161 0 : free(group);
2162 0 : return bmap;
2163 : }
2164 :
2165 : /* Drop constraints from "context" that are irrelevant for computing
2166 : * the gist of "bset".
2167 : *
2168 : * In particular, drop constraints in variables that are not related
2169 : * to any of the variables involved in the constraints of "bset"
2170 : * in the sense that there is no sequence of constraints that connects them.
2171 : *
2172 : * We first mark all variables that appear in "bset" as belonging
2173 : * to a "-1" group and then continue with group_and_drop_irrelevant_constraints.
2174 : */
2175 0 : static __isl_give isl_basic_set *drop_irrelevant_constraints(
2176 : __isl_take isl_basic_set *context, __isl_keep isl_basic_set *bset)
2177 : {
2178 : int *group;
2179 : int dim;
2180 : int i, j;
2181 :
2182 0 : if (!context || !bset)
2183 0 : return isl_basic_set_free(context);
2184 :
2185 0 : group = alloc_groups(context);
2186 :
2187 0 : if (!group)
2188 0 : return isl_basic_set_free(context);
2189 :
2190 0 : dim = isl_basic_set_dim(bset, isl_dim_set);
2191 0 : for (i = 0; i < dim; ++i) {
2192 0 : for (j = 0; j < bset->n_eq; ++j)
2193 0 : if (!isl_int_is_zero(bset->eq[j][1 + i]))
2194 : break;
2195 0 : if (j < bset->n_eq) {
2196 0 : group[i] = -1;
2197 0 : continue;
2198 : }
2199 0 : for (j = 0; j < bset->n_ineq; ++j)
2200 0 : if (!isl_int_is_zero(bset->ineq[j][1 + i]))
2201 : break;
2202 0 : if (j < bset->n_ineq)
2203 0 : group[i] = -1;
2204 : }
2205 :
2206 0 : return isl_basic_map_drop_unrelated_constraints(context, group);
2207 : }
2208 :
2209 : /* Drop constraints from "context" that are irrelevant for computing
2210 : * the gist of the inequalities "ineq".
2211 : * Inequalities in "ineq" for which the corresponding element of row
2212 : * is set to -1 have already been marked for removal and should be ignored.
2213 : *
2214 : * In particular, drop constraints in variables that are not related
2215 : * to any of the variables involved in "ineq"
2216 : * in the sense that there is no sequence of constraints that connects them.
2217 : *
2218 : * We first mark all variables that appear in "bset" as belonging
2219 : * to a "-1" group and then continue with group_and_drop_irrelevant_constraints.
2220 : */
2221 0 : static __isl_give isl_basic_set *drop_irrelevant_constraints_marked(
2222 : __isl_take isl_basic_set *context, __isl_keep isl_mat *ineq, int *row)
2223 : {
2224 : int *group;
2225 : int dim;
2226 : int i, j, n;
2227 :
2228 0 : if (!context || !ineq)
2229 0 : return isl_basic_set_free(context);
2230 :
2231 0 : group = alloc_groups(context);
2232 :
2233 0 : if (!group)
2234 0 : return isl_basic_set_free(context);
2235 :
2236 0 : dim = isl_basic_set_dim(context, isl_dim_set);
2237 0 : n = isl_mat_rows(ineq);
2238 0 : for (i = 0; i < dim; ++i) {
2239 0 : for (j = 0; j < n; ++j) {
2240 0 : if (row[j] < 0)
2241 0 : continue;
2242 0 : if (!isl_int_is_zero(ineq->row[j][1 + i]))
2243 : break;
2244 : }
2245 0 : if (j < n)
2246 0 : group[i] = -1;
2247 : }
2248 :
2249 0 : return isl_basic_map_drop_unrelated_constraints(context, group);
2250 : }
2251 :
2252 : /* Do all "n" entries of "row" contain a negative value?
2253 : */
2254 0 : static int all_neg(int *row, int n)
2255 : {
2256 : int i;
2257 :
2258 0 : for (i = 0; i < n; ++i)
2259 0 : if (row[i] >= 0)
2260 0 : return 0;
2261 :
2262 0 : return 1;
2263 : }
2264 :
2265 : /* Update the inequalities in "bset" based on the information in "row"
2266 : * and "tab".
2267 : *
2268 : * In particular, the array "row" contains either -1, meaning that
2269 : * the corresponding inequality of "bset" is redundant, or the index
2270 : * of an inequality in "tab".
2271 : *
2272 : * If the row entry is -1, then drop the inequality.
2273 : * Otherwise, if the constraint is marked redundant in the tableau,
2274 : * then drop the inequality. Similarly, if it is marked as an equality
2275 : * in the tableau, then turn the inequality into an equality and
2276 : * perform Gaussian elimination.
2277 : */
2278 0 : static __isl_give isl_basic_set *update_ineq(__isl_take isl_basic_set *bset,
2279 : __isl_keep int *row, struct isl_tab *tab)
2280 : {
2281 : int i;
2282 : unsigned n_ineq;
2283 : unsigned n_eq;
2284 0 : int found_equality = 0;
2285 :
2286 0 : if (!bset)
2287 0 : return NULL;
2288 0 : if (tab && tab->empty)
2289 0 : return isl_basic_set_set_to_empty(bset);
2290 :
2291 0 : n_ineq = bset->n_ineq;
2292 0 : for (i = n_ineq - 1; i >= 0; --i) {
2293 0 : if (row[i] < 0) {
2294 0 : if (isl_basic_set_drop_inequality(bset, i) < 0)
2295 0 : return isl_basic_set_free(bset);
2296 0 : continue;
2297 : }
2298 0 : if (!tab)
2299 0 : continue;
2300 0 : n_eq = tab->n_eq;
2301 0 : if (isl_tab_is_equality(tab, n_eq + row[i])) {
2302 0 : isl_basic_map_inequality_to_equality(bset, i);
2303 0 : found_equality = 1;
2304 0 : } else if (isl_tab_is_redundant(tab, n_eq + row[i])) {
2305 0 : if (isl_basic_set_drop_inequality(bset, i) < 0)
2306 0 : return isl_basic_set_free(bset);
2307 : }
2308 : }
2309 :
2310 0 : if (found_equality)
2311 0 : bset = isl_basic_set_gauss(bset, NULL);
2312 0 : bset = isl_basic_set_finalize(bset);
2313 0 : return bset;
2314 : }
2315 :
2316 : /* Update the inequalities in "bset" based on the information in "row"
2317 : * and "tab" and free all arguments (other than "bset").
2318 : */
2319 0 : static __isl_give isl_basic_set *update_ineq_free(
2320 : __isl_take isl_basic_set *bset, __isl_take isl_mat *ineq,
2321 : __isl_take isl_basic_set *context, __isl_take int *row,
2322 : struct isl_tab *tab)
2323 : {
2324 0 : isl_mat_free(ineq);
2325 0 : isl_basic_set_free(context);
2326 :
2327 0 : bset = update_ineq(bset, row, tab);
2328 :
2329 0 : free(row);
2330 0 : isl_tab_free(tab);
2331 0 : return bset;
2332 : }
2333 :
2334 : /* Remove all information from bset that is redundant in the context
2335 : * of context.
2336 : * "ineq" contains the (possibly transformed) inequalities of "bset",
2337 : * in the same order.
2338 : * The (explicit) equalities of "bset" are assumed to have been taken
2339 : * into account by the transformation such that only the inequalities
2340 : * are relevant.
2341 : * "context" is assumed not to be empty.
2342 : *
2343 : * "row" keeps track of the constraint index of a "bset" inequality in "tab".
2344 : * A value of -1 means that the inequality is obviously redundant and may
2345 : * not even appear in "tab".
2346 : *
2347 : * We first mark the inequalities of "bset"
2348 : * that are obviously redundant with respect to some inequality in "context".
2349 : * Then we remove those constraints from "context" that have become
2350 : * irrelevant for computing the gist of "bset".
2351 : * Note that this removal of constraints cannot be replaced by
2352 : * a factorization because factors in "bset" may still be connected
2353 : * to each other through constraints in "context".
2354 : *
2355 : * If there are any inequalities left, we construct a tableau for
2356 : * the context and then add the inequalities of "bset".
2357 : * Before adding these inequalities, we freeze all constraints such that
2358 : * they won't be considered redundant in terms of the constraints of "bset".
2359 : * Then we detect all redundant constraints (among the
2360 : * constraints that weren't frozen), first by checking for redundancy in the
2361 : * the tableau and then by checking if replacing a constraint by its negation
2362 : * would lead to an empty set. This last step is fairly expensive
2363 : * and could be optimized by more reuse of the tableau.
2364 : * Finally, we update bset according to the results.
2365 : */
2366 0 : static __isl_give isl_basic_set *uset_gist_full(__isl_take isl_basic_set *bset,
2367 : __isl_take isl_mat *ineq, __isl_take isl_basic_set *context)
2368 : {
2369 : int i, r;
2370 0 : int *row = NULL;
2371 : isl_ctx *ctx;
2372 0 : isl_basic_set *combined = NULL;
2373 0 : struct isl_tab *tab = NULL;
2374 : unsigned n_eq, context_ineq;
2375 :
2376 0 : if (!bset || !ineq || !context)
2377 : goto error;
2378 :
2379 0 : if (bset->n_ineq == 0 || isl_basic_set_plain_is_universe(context)) {
2380 0 : isl_basic_set_free(context);
2381 0 : isl_mat_free(ineq);
2382 0 : return bset;
2383 : }
2384 :
2385 0 : ctx = isl_basic_set_get_ctx(context);
2386 0 : row = isl_calloc_array(ctx, int, bset->n_ineq);
2387 0 : if (!row)
2388 0 : goto error;
2389 :
2390 0 : if (mark_shifted_constraints(ineq, context, row) < 0)
2391 0 : goto error;
2392 0 : if (all_neg(row, bset->n_ineq))
2393 0 : return update_ineq_free(bset, ineq, context, row, NULL);
2394 :
2395 0 : context = drop_irrelevant_constraints_marked(context, ineq, row);
2396 0 : if (!context)
2397 0 : goto error;
2398 0 : if (isl_basic_set_plain_is_universe(context))
2399 0 : return update_ineq_free(bset, ineq, context, row, NULL);
2400 :
2401 0 : n_eq = context->n_eq;
2402 0 : context_ineq = context->n_ineq;
2403 0 : combined = isl_basic_set_cow(isl_basic_set_copy(context));
2404 0 : combined = isl_basic_set_extend_constraints(combined, 0, bset->n_ineq);
2405 0 : tab = isl_tab_from_basic_set(combined, 0);
2406 0 : for (i = 0; i < context_ineq; ++i)
2407 0 : if (isl_tab_freeze_constraint(tab, n_eq + i) < 0)
2408 0 : goto error;
2409 0 : if (isl_tab_extend_cons(tab, bset->n_ineq) < 0)
2410 0 : goto error;
2411 0 : r = context_ineq;
2412 0 : for (i = 0; i < bset->n_ineq; ++i) {
2413 0 : if (row[i] < 0)
2414 0 : continue;
2415 0 : combined = isl_basic_set_add_ineq(combined, ineq->row[i]);
2416 0 : if (isl_tab_add_ineq(tab, ineq->row[i]) < 0)
2417 0 : goto error;
2418 0 : row[i] = r++;
2419 : }
2420 0 : if (isl_tab_detect_implicit_equalities(tab) < 0)
2421 0 : goto error;
2422 0 : if (isl_tab_detect_redundant(tab) < 0)
2423 0 : goto error;
2424 0 : for (i = bset->n_ineq - 1; i >= 0; --i) {
2425 : isl_basic_set *test;
2426 : int is_empty;
2427 :
2428 0 : if (row[i] < 0)
2429 0 : continue;
2430 0 : r = row[i];
2431 0 : if (tab->con[n_eq + r].is_redundant)
2432 0 : continue;
2433 0 : test = isl_basic_set_dup(combined);
2434 0 : if (isl_inequality_negate(test, r) < 0)
2435 0 : test = isl_basic_set_free(test);
2436 0 : test = isl_basic_set_update_from_tab(test, tab);
2437 0 : is_empty = isl_basic_set_is_empty(test);
2438 0 : isl_basic_set_free(test);
2439 0 : if (is_empty < 0)
2440 0 : goto error;
2441 0 : if (is_empty)
2442 0 : tab->con[n_eq + r].is_redundant = 1;
2443 : }
2444 0 : bset = update_ineq_free(bset, ineq, context, row, tab);
2445 0 : if (bset) {
2446 0 : ISL_F_SET(bset, ISL_BASIC_SET_NO_IMPLICIT);
2447 0 : ISL_F_SET(bset, ISL_BASIC_SET_NO_REDUNDANT);
2448 : }
2449 :
2450 0 : isl_basic_set_free(combined);
2451 0 : return bset;
2452 : error:
2453 0 : free(row);
2454 0 : isl_mat_free(ineq);
2455 0 : isl_tab_free(tab);
2456 0 : isl_basic_set_free(combined);
2457 0 : isl_basic_set_free(context);
2458 0 : isl_basic_set_free(bset);
2459 0 : return NULL;
2460 : }
2461 :
2462 : /* Extract the inequalities of "bset" as an isl_mat.
2463 : */
2464 0 : static __isl_give isl_mat *extract_ineq(__isl_keep isl_basic_set *bset)
2465 : {
2466 : unsigned total;
2467 : isl_ctx *ctx;
2468 : isl_mat *ineq;
2469 :
2470 0 : if (!bset)
2471 0 : return NULL;
2472 :
2473 0 : ctx = isl_basic_set_get_ctx(bset);
2474 0 : total = isl_basic_set_total_dim(bset);
2475 0 : ineq = isl_mat_sub_alloc6(ctx, bset->ineq, 0, bset->n_ineq,
2476 : 0, 1 + total);
2477 :
2478 0 : return ineq;
2479 : }
2480 :
2481 : /* Remove all information from "bset" that is redundant in the context
2482 : * of "context", for the case where both "bset" and "context" are
2483 : * full-dimensional.
2484 : */
2485 0 : static __isl_give isl_basic_set *uset_gist_uncompressed(
2486 : __isl_take isl_basic_set *bset, __isl_take isl_basic_set *context)
2487 : {
2488 : isl_mat *ineq;
2489 :
2490 0 : ineq = extract_ineq(bset);
2491 0 : return uset_gist_full(bset, ineq, context);
2492 : }
2493 :
2494 : /* Remove all information from "bset" that is redundant in the context
2495 : * of "context", for the case where the combined equalities of
2496 : * "bset" and "context" allow for a compression that can be obtained
2497 : * by preapplication of "T".
2498 : *
2499 : * "bset" itself is not transformed by "T". Instead, the inequalities
2500 : * are extracted from "bset" and those are transformed by "T".
2501 : * uset_gist_full then determines which of the transformed inequalities
2502 : * are redundant with respect to the transformed "context" and removes
2503 : * the corresponding inequalities from "bset".
2504 : *
2505 : * After preapplying "T" to the inequalities, any common factor is
2506 : * removed from the coefficients. If this results in a tightening
2507 : * of the constant term, then the same tightening is applied to
2508 : * the corresponding untransformed inequality in "bset".
2509 : * That is, if after plugging in T, a constraint f(x) >= 0 is of the form
2510 : *
2511 : * g f'(x) + r >= 0
2512 : *
2513 : * with 0 <= r < g, then it is equivalent to
2514 : *
2515 : * f'(x) >= 0
2516 : *
2517 : * This means that f(x) >= 0 is equivalent to f(x) - r >= 0 in the affine
2518 : * subspace compressed by T since the latter would be transformed to
2519 : *
2520 : * g f'(x) >= 0
2521 : */
2522 0 : static __isl_give isl_basic_set *uset_gist_compressed(
2523 : __isl_take isl_basic_set *bset, __isl_take isl_basic_set *context,
2524 : __isl_take isl_mat *T)
2525 : {
2526 : isl_ctx *ctx;
2527 : isl_mat *ineq;
2528 : int i, n_row, n_col;
2529 : isl_int rem;
2530 :
2531 0 : ineq = extract_ineq(bset);
2532 0 : ineq = isl_mat_product(ineq, isl_mat_copy(T));
2533 0 : context = isl_basic_set_preimage(context, T);
2534 :
2535 0 : if (!ineq || !context)
2536 : goto error;
2537 0 : if (isl_basic_set_plain_is_empty(context)) {
2538 0 : isl_mat_free(ineq);
2539 0 : isl_basic_set_free(context);
2540 0 : return isl_basic_set_set_to_empty(bset);
2541 : }
2542 :
2543 0 : ctx = isl_mat_get_ctx(ineq);
2544 0 : n_row = isl_mat_rows(ineq);
2545 0 : n_col = isl_mat_cols(ineq);
2546 0 : isl_int_init(rem);
2547 0 : for (i = 0; i < n_row; ++i) {
2548 0 : isl_seq_gcd(ineq->row[i] + 1, n_col - 1, &ctx->normalize_gcd);
2549 0 : if (isl_int_is_zero(ctx->normalize_gcd))
2550 0 : continue;
2551 0 : if (isl_int_is_one(ctx->normalize_gcd))
2552 0 : continue;
2553 0 : isl_seq_scale_down(ineq->row[i] + 1, ineq->row[i] + 1,
2554 0 : ctx->normalize_gcd, n_col - 1);
2555 0 : isl_int_fdiv_r(rem, ineq->row[i][0], ctx->normalize_gcd);
2556 0 : isl_int_fdiv_q(ineq->row[i][0],
2557 : ineq->row[i][0], ctx->normalize_gcd);
2558 0 : if (isl_int_is_zero(rem))
2559 0 : continue;
2560 0 : bset = isl_basic_set_cow(bset);
2561 0 : if (!bset)
2562 0 : break;
2563 0 : isl_int_sub(bset->ineq[i][0], bset->ineq[i][0], rem);
2564 : }
2565 0 : isl_int_clear(rem);
2566 :
2567 0 : return uset_gist_full(bset, ineq, context);
2568 : error:
2569 0 : isl_mat_free(ineq);
2570 0 : isl_basic_set_free(context);
2571 0 : isl_basic_set_free(bset);
2572 0 : return NULL;
2573 : }
2574 :
2575 : /* Project "bset" onto the variables that are involved in "template".
2576 : */
2577 0 : static __isl_give isl_basic_set *project_onto_involved(
2578 : __isl_take isl_basic_set *bset, __isl_keep isl_basic_set *template)
2579 : {
2580 : int i, n;
2581 :
2582 0 : if (!bset || !template)
2583 0 : return isl_basic_set_free(bset);
2584 :
2585 0 : n = isl_basic_set_dim(template, isl_dim_set);
2586 :
2587 0 : for (i = 0; i < n; ++i) {
2588 : isl_bool involved;
2589 :
2590 0 : involved = isl_basic_set_involves_dims(template,
2591 : isl_dim_set, i, 1);
2592 0 : if (involved < 0)
2593 0 : return isl_basic_set_free(bset);
2594 0 : if (involved)
2595 0 : continue;
2596 0 : bset = isl_basic_set_eliminate_vars(bset, i, 1);
2597 : }
2598 :
2599 0 : return bset;
2600 : }
2601 :
2602 : /* Remove all information from bset that is redundant in the context
2603 : * of context. In particular, equalities that are linear combinations
2604 : * of those in context are removed. Then the inequalities that are
2605 : * redundant in the context of the equalities and inequalities of
2606 : * context are removed.
2607 : *
2608 : * First of all, we drop those constraints from "context"
2609 : * that are irrelevant for computing the gist of "bset".
2610 : * Alternatively, we could factorize the intersection of "context" and "bset".
2611 : *
2612 : * We first compute the intersection of the integer affine hulls
2613 : * of "bset" and "context",
2614 : * compute the gist inside this intersection and then reduce
2615 : * the constraints with respect to the equalities of the context
2616 : * that only involve variables already involved in the input.
2617 : *
2618 : * If two constraints are mutually redundant, then uset_gist_full
2619 : * will remove the second of those constraints. We therefore first
2620 : * sort the constraints so that constraints not involving existentially
2621 : * quantified variables are given precedence over those that do.
2622 : * We have to perform this sorting before the variable compression,
2623 : * because that may effect the order of the variables.
2624 : */
2625 0 : static __isl_give isl_basic_set *uset_gist(__isl_take isl_basic_set *bset,
2626 : __isl_take isl_basic_set *context)
2627 : {
2628 : isl_mat *eq;
2629 : isl_mat *T;
2630 : isl_basic_set *aff;
2631 : isl_basic_set *aff_context;
2632 : unsigned total;
2633 :
2634 0 : if (!bset || !context)
2635 : goto error;
2636 :
2637 0 : context = drop_irrelevant_constraints(context, bset);
2638 :
2639 0 : bset = isl_basic_set_detect_equalities(bset);
2640 0 : aff = isl_basic_set_copy(bset);
2641 0 : aff = isl_basic_set_plain_affine_hull(aff);
2642 0 : context = isl_basic_set_detect_equalities(context);
2643 0 : aff_context = isl_basic_set_copy(context);
2644 0 : aff_context = isl_basic_set_plain_affine_hull(aff_context);
2645 0 : aff = isl_basic_set_intersect(aff, aff_context);
2646 0 : if (!aff)
2647 0 : goto error;
2648 0 : if (isl_basic_set_plain_is_empty(aff)) {
2649 0 : isl_basic_set_free(bset);
2650 0 : isl_basic_set_free(context);
2651 0 : return aff;
2652 : }
2653 0 : bset = isl_basic_set_sort_constraints(bset);
2654 0 : if (aff->n_eq == 0) {
2655 0 : isl_basic_set_free(aff);
2656 0 : return uset_gist_uncompressed(bset, context);
2657 : }
2658 0 : total = isl_basic_set_total_dim(bset);
2659 0 : eq = isl_mat_sub_alloc6(bset->ctx, aff->eq, 0, aff->n_eq, 0, 1 + total);
2660 0 : eq = isl_mat_cow(eq);
2661 0 : T = isl_mat_variable_compression(eq, NULL);
2662 0 : isl_basic_set_free(aff);
2663 0 : if (T && T->n_col == 0) {
2664 0 : isl_mat_free(T);
2665 0 : isl_basic_set_free(context);
2666 0 : return isl_basic_set_set_to_empty(bset);
2667 : }
2668 :
2669 0 : aff_context = isl_basic_set_affine_hull(isl_basic_set_copy(context));
2670 0 : aff_context = project_onto_involved(aff_context, bset);
2671 :
2672 0 : bset = uset_gist_compressed(bset, context, T);
2673 0 : bset = isl_basic_set_reduce_using_equalities(bset, aff_context);
2674 :
2675 0 : if (bset) {
2676 0 : ISL_F_SET(bset, ISL_BASIC_SET_NO_IMPLICIT);
2677 0 : ISL_F_SET(bset, ISL_BASIC_SET_NO_REDUNDANT);
2678 : }
2679 :
2680 0 : return bset;
2681 : error:
2682 0 : isl_basic_set_free(bset);
2683 0 : isl_basic_set_free(context);
2684 0 : return NULL;
2685 : }
2686 :
2687 : /* Return the number of equality constraints in "bmap" that involve
2688 : * local variables. This function assumes that Gaussian elimination
2689 : * has been applied to the equality constraints.
2690 : */
2691 0 : static int n_div_eq(__isl_keep isl_basic_map *bmap)
2692 : {
2693 : int i;
2694 : int total, n_div;
2695 :
2696 0 : if (!bmap)
2697 0 : return -1;
2698 :
2699 0 : if (bmap->n_eq == 0)
2700 0 : return 0;
2701 :
2702 0 : total = isl_basic_map_dim(bmap, isl_dim_all);
2703 0 : n_div = isl_basic_map_dim(bmap, isl_dim_div);
2704 0 : total -= n_div;
2705 :
2706 0 : for (i = 0; i < bmap->n_eq; ++i)
2707 0 : if (isl_seq_first_non_zero(bmap->eq[i] + 1 + total,
2708 : n_div) == -1)
2709 0 : return i;
2710 :
2711 0 : return bmap->n_eq;
2712 : }
2713 :
2714 : /* Construct a basic map in "space" defined by the equality constraints in "eq".
2715 : * The constraints are assumed not to involve any local variables.
2716 : */
2717 0 : static __isl_give isl_basic_map *basic_map_from_equalities(
2718 : __isl_take isl_space *space, __isl_take isl_mat *eq)
2719 : {
2720 : int i, k;
2721 0 : isl_basic_map *bmap = NULL;
2722 :
2723 0 : if (!space || !eq)
2724 : goto error;
2725 :
2726 0 : if (1 + isl_space_dim(space, isl_dim_all) != eq->n_col)
2727 0 : isl_die(isl_space_get_ctx(space), isl_error_internal,
2728 : "unexpected number of columns", goto error);
2729 :
2730 0 : bmap = isl_basic_map_alloc_space(isl_space_copy(space),
2731 : 0, eq->n_row, 0);
2732 0 : for (i = 0; i < eq->n_row; ++i) {
2733 0 : k = isl_basic_map_alloc_equality(bmap);
2734 0 : if (k < 0)
2735 0 : goto error;
2736 0 : isl_seq_cpy(bmap->eq[k], eq->row[i], eq->n_col);
2737 : }
2738 :
2739 0 : isl_space_free(space);
2740 0 : isl_mat_free(eq);
2741 0 : return bmap;
2742 : error:
2743 0 : isl_space_free(space);
2744 0 : isl_mat_free(eq);
2745 0 : isl_basic_map_free(bmap);
2746 0 : return NULL;
2747 : }
2748 :
2749 : /* Construct and return a variable compression based on the equality
2750 : * constraints in "bmap1" and "bmap2" that do not involve the local variables.
2751 : * "n1" is the number of (initial) equality constraints in "bmap1"
2752 : * that do involve local variables.
2753 : * "n2" is the number of (initial) equality constraints in "bmap2"
2754 : * that do involve local variables.
2755 : * "total" is the total number of other variables.
2756 : * This function assumes that Gaussian elimination
2757 : * has been applied to the equality constraints in both "bmap1" and "bmap2"
2758 : * such that the equality constraints not involving local variables
2759 : * are those that start at "n1" or "n2".
2760 : *
2761 : * If either of "bmap1" and "bmap2" does not have such equality constraints,
2762 : * then simply compute the compression based on the equality constraints
2763 : * in the other basic map.
2764 : * Otherwise, combine the equality constraints from both into a new
2765 : * basic map such that Gaussian elimination can be applied to this combination
2766 : * and then construct a variable compression from the resulting
2767 : * equality constraints.
2768 : */
2769 0 : static __isl_give isl_mat *combined_variable_compression(
2770 : __isl_keep isl_basic_map *bmap1, int n1,
2771 : __isl_keep isl_basic_map *bmap2, int n2, int total)
2772 : {
2773 : isl_ctx *ctx;
2774 : isl_mat *E1, *E2, *V;
2775 : isl_basic_map *bmap;
2776 :
2777 0 : ctx = isl_basic_map_get_ctx(bmap1);
2778 0 : if (bmap1->n_eq == n1) {
2779 0 : E2 = isl_mat_sub_alloc6(ctx, bmap2->eq,
2780 0 : n2, bmap2->n_eq - n2, 0, 1 + total);
2781 0 : return isl_mat_variable_compression(E2, NULL);
2782 : }
2783 0 : if (bmap2->n_eq == n2) {
2784 0 : E1 = isl_mat_sub_alloc6(ctx, bmap1->eq,
2785 0 : n1, bmap1->n_eq - n1, 0, 1 + total);
2786 0 : return isl_mat_variable_compression(E1, NULL);
2787 : }
2788 0 : E1 = isl_mat_sub_alloc6(ctx, bmap1->eq,
2789 0 : n1, bmap1->n_eq - n1, 0, 1 + total);
2790 0 : E2 = isl_mat_sub_alloc6(ctx, bmap2->eq,
2791 0 : n2, bmap2->n_eq - n2, 0, 1 + total);
2792 0 : E1 = isl_mat_concat(E1, E2);
2793 0 : bmap = basic_map_from_equalities(isl_basic_map_get_space(bmap1), E1);
2794 0 : bmap = isl_basic_map_gauss(bmap, NULL);
2795 0 : if (!bmap)
2796 0 : return NULL;
2797 0 : E1 = isl_mat_sub_alloc6(ctx, bmap->eq, 0, bmap->n_eq, 0, 1 + total);
2798 0 : V = isl_mat_variable_compression(E1, NULL);
2799 0 : isl_basic_map_free(bmap);
2800 :
2801 0 : return V;
2802 : }
2803 :
2804 : /* Extract the stride constraints from "bmap", compressed
2805 : * with respect to both the stride constraints in "context" and
2806 : * the remaining equality constraints in both "bmap" and "context".
2807 : * "bmap_n_eq" is the number of (initial) stride constraints in "bmap".
2808 : * "context_n_eq" is the number of (initial) stride constraints in "context".
2809 : *
2810 : * Let x be all variables in "bmap" (and "context") other than the local
2811 : * variables. First compute a variable compression
2812 : *
2813 : * x = V x'
2814 : *
2815 : * based on the non-stride equality constraints in "bmap" and "context".
2816 : * Consider the stride constraints of "context",
2817 : *
2818 : * A(x) + B(y) = 0
2819 : *
2820 : * with y the local variables and plug in the variable compression,
2821 : * resulting in
2822 : *
2823 : * A(V x') + B(y) = 0
2824 : *
2825 : * Use these constraints to compute a parameter compression on x'
2826 : *
2827 : * x' = T x''
2828 : *
2829 : * Now consider the stride constraints of "bmap"
2830 : *
2831 : * C(x) + D(y) = 0
2832 : *
2833 : * and plug in x = V*T x''.
2834 : * That is, return A = [C*V*T D].
2835 : */
2836 0 : static __isl_give isl_mat *extract_compressed_stride_constraints(
2837 : __isl_keep isl_basic_map *bmap, int bmap_n_eq,
2838 : __isl_keep isl_basic_map *context, int context_n_eq)
2839 : {
2840 : int total, n_div;
2841 : isl_ctx *ctx;
2842 : isl_mat *A, *B, *T, *V;
2843 :
2844 0 : total = isl_basic_map_dim(context, isl_dim_all);
2845 0 : n_div = isl_basic_map_dim(context, isl_dim_div);
2846 0 : total -= n_div;
2847 :
2848 0 : ctx = isl_basic_map_get_ctx(bmap);
2849 :
2850 0 : V = combined_variable_compression(bmap, bmap_n_eq,
2851 : context, context_n_eq, total);
2852 :
2853 0 : A = isl_mat_sub_alloc6(ctx, context->eq, 0, context_n_eq, 0, 1 + total);
2854 0 : B = isl_mat_sub_alloc6(ctx, context->eq,
2855 0 : 0, context_n_eq, 1 + total, n_div);
2856 0 : A = isl_mat_product(A, isl_mat_copy(V));
2857 0 : T = isl_mat_parameter_compression_ext(A, B);
2858 0 : T = isl_mat_product(V, T);
2859 :
2860 0 : n_div = isl_basic_map_dim(bmap, isl_dim_div);
2861 0 : T = isl_mat_diagonal(T, isl_mat_identity(ctx, n_div));
2862 :
2863 0 : A = isl_mat_sub_alloc6(ctx, bmap->eq,
2864 0 : 0, bmap_n_eq, 0, 1 + total + n_div);
2865 0 : A = isl_mat_product(A, T);
2866 :
2867 0 : return A;
2868 : }
2869 :
2870 : /* Remove the prime factors from *g that have an exponent that
2871 : * is strictly smaller than the exponent in "c".
2872 : * All exponents in *g are known to be smaller than or equal
2873 : * to those in "c".
2874 : *
2875 : * That is, if *g is equal to
2876 : *
2877 : * p_1^{e_1} p_2^{e_2} ... p_n^{e_n}
2878 : *
2879 : * and "c" is equal to
2880 : *
2881 : * p_1^{f_1} p_2^{f_2} ... p_n^{f_n}
2882 : *
2883 : * then update *g to
2884 : *
2885 : * p_1^{e_1 * (e_1 = f_1)} p_2^{e_2 * (e_2 = f_2)} ...
2886 : * p_n^{e_n * (e_n = f_n)}
2887 : *
2888 : * If e_i = f_i, then c / *g does not have any p_i factors and therefore
2889 : * neither does the gcd of *g and c / *g.
2890 : * If e_i < f_i, then the gcd of *g and c / *g has a positive
2891 : * power min(e_i, s_i) of p_i with s_i = f_i - e_i among its factors.
2892 : * Dividing *g by this gcd therefore strictly reduces the exponent
2893 : * of the prime factors that need to be removed, while leaving the
2894 : * other prime factors untouched.
2895 : * Repeating this process until gcd(*g, c / *g) = 1 therefore
2896 : * removes all undesired factors, without removing any others.
2897 : */
2898 0 : static void remove_incomplete_powers(isl_int *g, isl_int c)
2899 : {
2900 : isl_int t;
2901 :
2902 0 : isl_int_init(t);
2903 : for (;;) {
2904 0 : isl_int_divexact(t, c, *g);
2905 0 : isl_int_gcd(t, t, *g);
2906 0 : if (isl_int_is_one(t))
2907 0 : break;
2908 0 : isl_int_divexact(*g, *g, t);
2909 0 : }
2910 0 : isl_int_clear(t);
2911 0 : }
2912 :
2913 : /* Reduce the "n" stride constraints in "bmap" based on a copy "A"
2914 : * of the same stride constraints in a compressed space that exploits
2915 : * all equalities in the context and the other equalities in "bmap".
2916 : *
2917 : * If the stride constraints of "bmap" are of the form
2918 : *
2919 : * C(x) + D(y) = 0
2920 : *
2921 : * then A is of the form
2922 : *
2923 : * B(x') + D(y) = 0
2924 : *
2925 : * If any of these constraints involves only a single local variable y,
2926 : * then the constraint appears as
2927 : *
2928 : * f(x) + m y_i = 0
2929 : *
2930 : * in "bmap" and as
2931 : *
2932 : * h(x') + m y_i = 0
2933 : *
2934 : * in "A".
2935 : *
2936 : * Let g be the gcd of m and the coefficients of h.
2937 : * Then, in particular, g is a divisor of the coefficients of h and
2938 : *
2939 : * f(x) = h(x')
2940 : *
2941 : * is known to be a multiple of g.
2942 : * If some prime factor in m appears with the same exponent in g,
2943 : * then it can be removed from m because f(x) is already known
2944 : * to be a multiple of g and therefore in particular of this power
2945 : * of the prime factors.
2946 : * Prime factors that appear with a smaller exponent in g cannot
2947 : * be removed from m.
2948 : * Let g' be the divisor of g containing all prime factors that
2949 : * appear with the same exponent in m and g, then
2950 : *
2951 : * f(x) + m y_i = 0
2952 : *
2953 : * can be replaced by
2954 : *
2955 : * f(x) + m/g' y_i' = 0
2956 : *
2957 : * Note that (if g' != 1) this changes the explicit representation
2958 : * of y_i to that of y_i', so the integer division at position i
2959 : * is marked unknown and later recomputed by a call to
2960 : * isl_basic_map_gauss.
2961 : */
2962 0 : static __isl_give isl_basic_map *reduce_stride_constraints(
2963 : __isl_take isl_basic_map *bmap, int n, __isl_keep isl_mat *A)
2964 : {
2965 : int i;
2966 : int total, n_div;
2967 0 : int any = 0;
2968 : isl_int gcd;
2969 :
2970 0 : if (!bmap || !A)
2971 0 : return isl_basic_map_free(bmap);
2972 :
2973 0 : total = isl_basic_map_dim(bmap, isl_dim_all);
2974 0 : n_div = isl_basic_map_dim(bmap, isl_dim_div);
2975 0 : total -= n_div;
2976 :
2977 0 : isl_int_init(gcd);
2978 0 : for (i = 0; i < n; ++i) {
2979 : int div;
2980 :
2981 0 : div = isl_seq_first_non_zero(bmap->eq[i] + 1 + total, n_div);
2982 0 : if (div < 0)
2983 0 : isl_die(isl_basic_map_get_ctx(bmap), isl_error_internal,
2984 : "equality constraints modified unexpectedly",
2985 : goto error);
2986 0 : if (isl_seq_first_non_zero(bmap->eq[i] + 1 + total + div + 1,
2987 0 : n_div - div - 1) != -1)
2988 0 : continue;
2989 0 : if (isl_mat_row_gcd(A, i, &gcd) < 0)
2990 0 : goto error;
2991 0 : if (isl_int_is_one(gcd))
2992 0 : continue;
2993 0 : remove_incomplete_powers(&gcd, bmap->eq[i][1 + total + div]);
2994 0 : if (isl_int_is_one(gcd))
2995 0 : continue;
2996 0 : isl_int_divexact(bmap->eq[i][1 + total + div],
2997 : bmap->eq[i][1 + total + div], gcd);
2998 0 : bmap = isl_basic_map_mark_div_unknown(bmap, div);
2999 0 : if (!bmap)
3000 0 : goto error;
3001 0 : any = 1;
3002 : }
3003 0 : isl_int_clear(gcd);
3004 :
3005 0 : if (any)
3006 0 : bmap = isl_basic_map_gauss(bmap, NULL);
3007 :
3008 0 : return bmap;
3009 : error:
3010 0 : isl_int_clear(gcd);
3011 0 : isl_basic_map_free(bmap);
3012 0 : return NULL;
3013 : }
3014 :
3015 : /* Simplify the stride constraints in "bmap" based on
3016 : * the remaining equality constraints in "bmap" and all equality
3017 : * constraints in "context".
3018 : * Only do this if both "bmap" and "context" have stride constraints.
3019 : *
3020 : * First extract a copy of the stride constraints in "bmap" in a compressed
3021 : * space exploiting all the other equality constraints and then
3022 : * use this compressed copy to simplify the original stride constraints.
3023 : */
3024 0 : static __isl_give isl_basic_map *gist_strides(__isl_take isl_basic_map *bmap,
3025 : __isl_keep isl_basic_map *context)
3026 : {
3027 : int bmap_n_eq, context_n_eq;
3028 : isl_mat *A;
3029 :
3030 0 : if (!bmap || !context)
3031 0 : return isl_basic_map_free(bmap);
3032 :
3033 0 : bmap_n_eq = n_div_eq(bmap);
3034 0 : context_n_eq = n_div_eq(context);
3035 :
3036 0 : if (bmap_n_eq < 0 || context_n_eq < 0)
3037 0 : return isl_basic_map_free(bmap);
3038 0 : if (bmap_n_eq == 0 || context_n_eq == 0)
3039 0 : return bmap;
3040 :
3041 0 : A = extract_compressed_stride_constraints(bmap, bmap_n_eq,
3042 : context, context_n_eq);
3043 0 : bmap = reduce_stride_constraints(bmap, bmap_n_eq, A);
3044 :
3045 0 : isl_mat_free(A);
3046 :
3047 0 : return bmap;
3048 : }
3049 :
3050 : /* Return a basic map that has the same intersection with "context" as "bmap"
3051 : * and that is as "simple" as possible.
3052 : *
3053 : * The core computation is performed on the pure constraints.
3054 : * When we add back the meaning of the integer divisions, we need
3055 : * to (re)introduce the div constraints. If we happen to have
3056 : * discovered that some of these integer divisions are equal to
3057 : * some affine combination of other variables, then these div
3058 : * constraints may end up getting simplified in terms of the equalities,
3059 : * resulting in extra inequalities on the other variables that
3060 : * may have been removed already or that may not even have been
3061 : * part of the input. We try and remove those constraints of
3062 : * this form that are most obviously redundant with respect to
3063 : * the context. We also remove those div constraints that are
3064 : * redundant with respect to the other constraints in the result.
3065 : *
3066 : * The stride constraints among the equality constraints in "bmap" are
3067 : * also simplified with respecting to the other equality constraints
3068 : * in "bmap" and with respect to all equality constraints in "context".
3069 : */
3070 0 : __isl_give isl_basic_map *isl_basic_map_gist(__isl_take isl_basic_map *bmap,
3071 : __isl_take isl_basic_map *context)
3072 : {
3073 : isl_basic_set *bset, *eq;
3074 : isl_basic_map *eq_bmap;
3075 : unsigned total, n_div, extra, n_eq, n_ineq;
3076 :
3077 0 : if (!bmap || !context)
3078 : goto error;
3079 :
3080 0 : if (isl_basic_map_plain_is_universe(bmap)) {
3081 0 : isl_basic_map_free(context);
3082 0 : return bmap;
3083 : }
3084 0 : if (isl_basic_map_plain_is_empty(context)) {
3085 0 : isl_space *space = isl_basic_map_get_space(bmap);
3086 0 : isl_basic_map_free(bmap);
3087 0 : isl_basic_map_free(context);
3088 0 : return isl_basic_map_universe(space);
3089 : }
3090 0 : if (isl_basic_map_plain_is_empty(bmap)) {
3091 0 : isl_basic_map_free(context);
3092 0 : return bmap;
3093 : }
3094 :
3095 0 : bmap = isl_basic_map_remove_redundancies(bmap);
3096 0 : context = isl_basic_map_remove_redundancies(context);
3097 0 : context = isl_basic_map_align_divs(context, bmap);
3098 0 : if (!context)
3099 0 : goto error;
3100 :
3101 0 : n_div = isl_basic_map_dim(context, isl_dim_div);
3102 0 : total = isl_basic_map_dim(bmap, isl_dim_all);
3103 0 : extra = n_div - isl_basic_map_dim(bmap, isl_dim_div);
3104 :
3105 0 : bset = isl_basic_map_underlying_set(isl_basic_map_copy(bmap));
3106 0 : bset = isl_basic_set_add_dims(bset, isl_dim_set, extra);
3107 0 : bset = uset_gist(bset,
3108 : isl_basic_map_underlying_set(isl_basic_map_copy(context)));
3109 0 : bset = isl_basic_set_project_out(bset, isl_dim_set, total, extra);
3110 :
3111 0 : if (!bset || bset->n_eq == 0 || n_div == 0 ||
3112 0 : isl_basic_set_plain_is_empty(bset)) {
3113 0 : isl_basic_map_free(context);
3114 0 : return isl_basic_map_overlying_set(bset, bmap);
3115 : }
3116 :
3117 0 : n_eq = bset->n_eq;
3118 0 : n_ineq = bset->n_ineq;
3119 0 : eq = isl_basic_set_copy(bset);
3120 0 : eq = isl_basic_set_cow(eq);
3121 0 : if (isl_basic_set_free_inequality(eq, n_ineq) < 0)
3122 0 : eq = isl_basic_set_free(eq);
3123 0 : if (isl_basic_set_free_equality(bset, n_eq) < 0)
3124 0 : bset = isl_basic_set_free(bset);
3125 :
3126 0 : eq_bmap = isl_basic_map_overlying_set(eq, isl_basic_map_copy(bmap));
3127 0 : eq_bmap = gist_strides(eq_bmap, context);
3128 0 : eq_bmap = isl_basic_map_remove_shifted_constraints(eq_bmap, context);
3129 0 : bmap = isl_basic_map_overlying_set(bset, bmap);
3130 0 : bmap = isl_basic_map_intersect(bmap, eq_bmap);
3131 0 : bmap = isl_basic_map_remove_redundancies(bmap);
3132 :
3133 0 : return bmap;
3134 : error:
3135 0 : isl_basic_map_free(bmap);
3136 0 : isl_basic_map_free(context);
3137 0 : return NULL;
3138 : }
3139 :
3140 : /*
3141 : * Assumes context has no implicit divs.
3142 : */
3143 0 : __isl_give isl_map *isl_map_gist_basic_map(__isl_take isl_map *map,
3144 : __isl_take isl_basic_map *context)
3145 : {
3146 : int i;
3147 :
3148 0 : if (!map || !context)
3149 : goto error;
3150 :
3151 0 : if (isl_basic_map_plain_is_empty(context)) {
3152 0 : isl_space *space = isl_map_get_space(map);
3153 0 : isl_map_free(map);
3154 0 : isl_basic_map_free(context);
3155 0 : return isl_map_universe(space);
3156 : }
3157 :
3158 0 : context = isl_basic_map_remove_redundancies(context);
3159 0 : map = isl_map_cow(map);
3160 0 : if (!map || !context)
3161 : goto error;
3162 0 : isl_assert(map->ctx, isl_space_is_equal(map->dim, context->dim), goto error);
3163 0 : map = isl_map_compute_divs(map);
3164 0 : if (!map)
3165 0 : goto error;
3166 0 : for (i = map->n - 1; i >= 0; --i) {
3167 0 : map->p[i] = isl_basic_map_gist(map->p[i],
3168 : isl_basic_map_copy(context));
3169 0 : if (!map->p[i])
3170 0 : goto error;
3171 0 : if (isl_basic_map_plain_is_empty(map->p[i])) {
3172 0 : isl_basic_map_free(map->p[i]);
3173 0 : if (i != map->n - 1)
3174 0 : map->p[i] = map->p[map->n - 1];
3175 0 : map->n--;
3176 : }
3177 : }
3178 0 : isl_basic_map_free(context);
3179 0 : ISL_F_CLR(map, ISL_MAP_NORMALIZED);
3180 0 : return map;
3181 : error:
3182 0 : isl_map_free(map);
3183 0 : isl_basic_map_free(context);
3184 0 : return NULL;
3185 : }
3186 :
3187 : /* Drop all inequalities from "bmap" that also appear in "context".
3188 : * "context" is assumed to have only known local variables and
3189 : * the initial local variables of "bmap" are assumed to be the same
3190 : * as those of "context".
3191 : * The constraints of both "bmap" and "context" are assumed
3192 : * to have been sorted using isl_basic_map_sort_constraints.
3193 : *
3194 : * Run through the inequality constraints of "bmap" and "context"
3195 : * in sorted order.
3196 : * If a constraint of "bmap" involves variables not in "context",
3197 : * then it cannot appear in "context".
3198 : * If a matching constraint is found, it is removed from "bmap".
3199 : */
3200 0 : static __isl_give isl_basic_map *drop_inequalities(
3201 : __isl_take isl_basic_map *bmap, __isl_keep isl_basic_map *context)
3202 : {
3203 : int i1, i2;
3204 : unsigned total, extra;
3205 :
3206 0 : if (!bmap || !context)
3207 0 : return isl_basic_map_free(bmap);
3208 :
3209 0 : total = isl_basic_map_total_dim(context);
3210 0 : extra = isl_basic_map_total_dim(bmap) - total;
3211 :
3212 0 : i1 = bmap->n_ineq - 1;
3213 0 : i2 = context->n_ineq - 1;
3214 0 : while (bmap && i1 >= 0 && i2 >= 0) {
3215 : int cmp;
3216 :
3217 0 : if (isl_seq_first_non_zero(bmap->ineq[i1] + 1 + total,
3218 : extra) != -1) {
3219 0 : --i1;
3220 0 : continue;
3221 : }
3222 0 : cmp = isl_basic_map_constraint_cmp(context, bmap->ineq[i1],
3223 0 : context->ineq[i2]);
3224 0 : if (cmp < 0) {
3225 0 : --i2;
3226 0 : continue;
3227 : }
3228 0 : if (cmp > 0) {
3229 0 : --i1;
3230 0 : continue;
3231 : }
3232 0 : if (isl_int_eq(bmap->ineq[i1][0], context->ineq[i2][0])) {
3233 0 : bmap = isl_basic_map_cow(bmap);
3234 0 : if (isl_basic_map_drop_inequality(bmap, i1) < 0)
3235 0 : bmap = isl_basic_map_free(bmap);
3236 : }
3237 0 : --i1;
3238 0 : --i2;
3239 : }
3240 :
3241 0 : return bmap;
3242 : }
3243 :
3244 : /* Drop all equalities from "bmap" that also appear in "context".
3245 : * "context" is assumed to have only known local variables and
3246 : * the initial local variables of "bmap" are assumed to be the same
3247 : * as those of "context".
3248 : *
3249 : * Run through the equality constraints of "bmap" and "context"
3250 : * in sorted order.
3251 : * If a constraint of "bmap" involves variables not in "context",
3252 : * then it cannot appear in "context".
3253 : * If a matching constraint is found, it is removed from "bmap".
3254 : */
3255 0 : static __isl_give isl_basic_map *drop_equalities(
3256 : __isl_take isl_basic_map *bmap, __isl_keep isl_basic_map *context)
3257 : {
3258 : int i1, i2;
3259 : unsigned total, extra;
3260 :
3261 0 : if (!bmap || !context)
3262 0 : return isl_basic_map_free(bmap);
3263 :
3264 0 : total = isl_basic_map_total_dim(context);
3265 0 : extra = isl_basic_map_total_dim(bmap) - total;
3266 :
3267 0 : i1 = bmap->n_eq - 1;
3268 0 : i2 = context->n_eq - 1;
3269 :
3270 0 : while (bmap && i1 >= 0 && i2 >= 0) {
3271 : int last1, last2;
3272 :
3273 0 : if (isl_seq_first_non_zero(bmap->eq[i1] + 1 + total,
3274 : extra) != -1)
3275 0 : break;
3276 0 : last1 = isl_seq_last_non_zero(bmap->eq[i1] + 1, total);
3277 0 : last2 = isl_seq_last_non_zero(context->eq[i2] + 1, total);
3278 0 : if (last1 > last2) {
3279 0 : --i2;
3280 0 : continue;
3281 : }
3282 0 : if (last1 < last2) {
3283 0 : --i1;
3284 0 : continue;
3285 : }
3286 0 : if (isl_seq_eq(bmap->eq[i1], context->eq[i2], 1 + total)) {
3287 0 : bmap = isl_basic_map_cow(bmap);
3288 0 : if (isl_basic_map_drop_equality(bmap, i1) < 0)
3289 0 : bmap = isl_basic_map_free(bmap);
3290 : }
3291 0 : --i1;
3292 0 : --i2;
3293 : }
3294 :
3295 0 : return bmap;
3296 : }
3297 :
3298 : /* Remove the constraints in "context" from "bmap".
3299 : * "context" is assumed to have explicit representations
3300 : * for all local variables.
3301 : *
3302 : * First align the divs of "bmap" to those of "context" and
3303 : * sort the constraints. Then drop all constraints from "bmap"
3304 : * that appear in "context".
3305 : */
3306 0 : __isl_give isl_basic_map *isl_basic_map_plain_gist(
3307 : __isl_take isl_basic_map *bmap, __isl_take isl_basic_map *context)
3308 : {
3309 : isl_bool done, known;
3310 :
3311 0 : done = isl_basic_map_plain_is_universe(context);
3312 0 : if (done == isl_bool_false)
3313 0 : done = isl_basic_map_plain_is_universe(bmap);
3314 0 : if (done == isl_bool_false)
3315 0 : done = isl_basic_map_plain_is_empty(context);
3316 0 : if (done == isl_bool_false)
3317 0 : done = isl_basic_map_plain_is_empty(bmap);
3318 0 : if (done < 0)
3319 0 : goto error;
3320 0 : if (done) {
3321 0 : isl_basic_map_free(context);
3322 0 : return bmap;
3323 : }
3324 0 : known = isl_basic_map_divs_known(context);
3325 0 : if (known < 0)
3326 0 : goto error;
3327 0 : if (!known)
3328 0 : isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
3329 : "context has unknown divs", goto error);
3330 :
3331 0 : bmap = isl_basic_map_align_divs(bmap, context);
3332 0 : bmap = isl_basic_map_gauss(bmap, NULL);
3333 0 : bmap = isl_basic_map_sort_constraints(bmap);
3334 0 : context = isl_basic_map_sort_constraints(context);
3335 :
3336 0 : bmap = drop_inequalities(bmap, context);
3337 0 : bmap = drop_equalities(bmap, context);
3338 :
3339 0 : isl_basic_map_free(context);
3340 0 : bmap = isl_basic_map_finalize(bmap);
3341 0 : return bmap;
3342 : error:
3343 0 : isl_basic_map_free(bmap);
3344 0 : isl_basic_map_free(context);
3345 0 : return NULL;
3346 : }
3347 :
3348 : /* Replace "map" by the disjunct at position "pos" and free "context".
3349 : */
3350 0 : static __isl_give isl_map *replace_by_disjunct(__isl_take isl_map *map,
3351 : int pos, __isl_take isl_basic_map *context)
3352 : {
3353 : isl_basic_map *bmap;
3354 :
3355 0 : bmap = isl_basic_map_copy(map->p[pos]);
3356 0 : isl_map_free(map);
3357 0 : isl_basic_map_free(context);
3358 0 : return isl_map_from_basic_map(bmap);
3359 : }
3360 :
3361 : /* Remove the constraints in "context" from "map".
3362 : * If any of the disjuncts in the result turns out to be the universe,
3363 : * then return this universe.
3364 : * "context" is assumed to have explicit representations
3365 : * for all local variables.
3366 : */
3367 0 : __isl_give isl_map *isl_map_plain_gist_basic_map(__isl_take isl_map *map,
3368 : __isl_take isl_basic_map *context)
3369 : {
3370 : int i;
3371 : isl_bool univ, known;
3372 :
3373 0 : univ = isl_basic_map_plain_is_universe(context);
3374 0 : if (univ < 0)
3375 0 : goto error;
3376 0 : if (univ) {
3377 0 : isl_basic_map_free(context);
3378 0 : return map;
3379 : }
3380 0 : known = isl_basic_map_divs_known(context);
3381 0 : if (known < 0)
3382 0 : goto error;
3383 0 : if (!known)
3384 0 : isl_die(isl_map_get_ctx(map), isl_error_invalid,
3385 : "context has unknown divs", goto error);
3386 :
3387 0 : map = isl_map_cow(map);
3388 0 : if (!map)
3389 0 : goto error;
3390 0 : for (i = 0; i < map->n; ++i) {
3391 0 : map->p[i] = isl_basic_map_plain_gist(map->p[i],
3392 : isl_basic_map_copy(context));
3393 0 : univ = isl_basic_map_plain_is_universe(map->p[i]);
3394 0 : if (univ < 0)
3395 0 : goto error;
3396 0 : if (univ && map->n > 1)
3397 0 : return replace_by_disjunct(map, i, context);
3398 : }
3399 :
3400 0 : isl_basic_map_free(context);
3401 0 : ISL_F_CLR(map, ISL_MAP_NORMALIZED);
3402 0 : if (map->n > 1)
3403 0 : ISL_F_CLR(map, ISL_MAP_DISJOINT);
3404 0 : return map;
3405 : error:
3406 0 : isl_map_free(map);
3407 0 : isl_basic_map_free(context);
3408 0 : return NULL;
3409 : }
3410 :
3411 : /* Remove the constraints in "context" from "set".
3412 : * If any of the disjuncts in the result turns out to be the universe,
3413 : * then return this universe.
3414 : * "context" is assumed to have explicit representations
3415 : * for all local variables.
3416 : */
3417 0 : __isl_give isl_set *isl_set_plain_gist_basic_set(__isl_take isl_set *set,
3418 : __isl_take isl_basic_set *context)
3419 : {
3420 0 : return set_from_map(isl_map_plain_gist_basic_map(set_to_map(set),
3421 : bset_to_bmap(context)));
3422 : }
3423 :
3424 : /* Remove the constraints in "context" from "map".
3425 : * If any of the disjuncts in the result turns out to be the universe,
3426 : * then return this universe.
3427 : * "context" is assumed to consist of a single disjunct and
3428 : * to have explicit representations for all local variables.
3429 : */
3430 0 : __isl_give isl_map *isl_map_plain_gist(__isl_take isl_map *map,
3431 : __isl_take isl_map *context)
3432 : {
3433 : isl_basic_map *hull;
3434 :
3435 0 : hull = isl_map_unshifted_simple_hull(context);
3436 0 : return isl_map_plain_gist_basic_map(map, hull);
3437 : }
3438 :
3439 : /* Replace "map" by a universe map in the same space and free "drop".
3440 : */
3441 0 : static __isl_give isl_map *replace_by_universe(__isl_take isl_map *map,
3442 : __isl_take isl_map *drop)
3443 : {
3444 : isl_map *res;
3445 :
3446 0 : res = isl_map_universe(isl_map_get_space(map));
3447 0 : isl_map_free(map);
3448 0 : isl_map_free(drop);
3449 0 : return res;
3450 : }
3451 :
3452 : /* Return a map that has the same intersection with "context" as "map"
3453 : * and that is as "simple" as possible.
3454 : *
3455 : * If "map" is already the universe, then we cannot make it any simpler.
3456 : * Similarly, if "context" is the universe, then we cannot exploit it
3457 : * to simplify "map"
3458 : * If "map" and "context" are identical to each other, then we can
3459 : * return the corresponding universe.
3460 : *
3461 : * If either "map" or "context" consists of multiple disjuncts,
3462 : * then check if "context" happens to be a subset of "map",
3463 : * in which case all constraints can be removed.
3464 : * In case of multiple disjuncts, the standard procedure
3465 : * may not be able to detect that all constraints can be removed.
3466 : *
3467 : * If none of these cases apply, we have to work a bit harder.
3468 : * During this computation, we make use of a single disjunct context,
3469 : * so if the original context consists of more than one disjunct
3470 : * then we need to approximate the context by a single disjunct set.
3471 : * Simply taking the simple hull may drop constraints that are
3472 : * only implicitly available in each disjunct. We therefore also
3473 : * look for constraints among those defining "map" that are valid
3474 : * for the context. These can then be used to simplify away
3475 : * the corresponding constraints in "map".
3476 : */
3477 0 : static __isl_give isl_map *map_gist(__isl_take isl_map *map,
3478 : __isl_take isl_map *context)
3479 : {
3480 : int equal;
3481 : int is_universe;
3482 : int single_disjunct_map, single_disjunct_context;
3483 : isl_bool subset;
3484 : isl_basic_map *hull;
3485 :
3486 0 : is_universe = isl_map_plain_is_universe(map);
3487 0 : if (is_universe >= 0 && !is_universe)
3488 0 : is_universe = isl_map_plain_is_universe(context);
3489 0 : if (is_universe < 0)
3490 0 : goto error;
3491 0 : if (is_universe) {
3492 0 : isl_map_free(context);
3493 0 : return map;
3494 : }
3495 :
3496 0 : equal = isl_map_plain_is_equal(map, context);
3497 0 : if (equal < 0)
3498 0 : goto error;
3499 0 : if (equal)
3500 0 : return replace_by_universe(map, context);
3501 :
3502 0 : single_disjunct_map = isl_map_n_basic_map(map) == 1;
3503 0 : single_disjunct_context = isl_map_n_basic_map(context) == 1;
3504 0 : if (!single_disjunct_map || !single_disjunct_context) {
3505 0 : subset = isl_map_is_subset(context, map);
3506 0 : if (subset < 0)
3507 0 : goto error;
3508 0 : if (subset)
3509 0 : return replace_by_universe(map, context);
3510 : }
3511 :
3512 0 : context = isl_map_compute_divs(context);
3513 0 : if (!context)
3514 0 : goto error;
3515 0 : if (single_disjunct_context) {
3516 0 : hull = isl_map_simple_hull(context);
3517 : } else {
3518 : isl_ctx *ctx;
3519 : isl_map_list *list;
3520 :
3521 0 : ctx = isl_map_get_ctx(map);
3522 0 : list = isl_map_list_alloc(ctx, 2);
3523 0 : list = isl_map_list_add(list, isl_map_copy(context));
3524 0 : list = isl_map_list_add(list, isl_map_copy(map));
3525 0 : hull = isl_map_unshifted_simple_hull_from_map_list(context,
3526 : list);
3527 : }
3528 0 : return isl_map_gist_basic_map(map, hull);
3529 : error:
3530 0 : isl_map_free(map);
3531 0 : isl_map_free(context);
3532 0 : return NULL;
3533 : }
3534 :
3535 0 : __isl_give isl_map *isl_map_gist(__isl_take isl_map *map,
3536 : __isl_take isl_map *context)
3537 : {
3538 0 : return isl_map_align_params_map_map_and(map, context, &map_gist);
3539 : }
3540 :
3541 0 : struct isl_basic_set *isl_basic_set_gist(struct isl_basic_set *bset,
3542 : struct isl_basic_set *context)
3543 : {
3544 0 : return bset_from_bmap(isl_basic_map_gist(bset_to_bmap(bset),
3545 : bset_to_bmap(context)));
3546 : }
3547 :
3548 0 : __isl_give isl_set *isl_set_gist_basic_set(__isl_take isl_set *set,
3549 : __isl_take isl_basic_set *context)
3550 : {
3551 0 : return set_from_map(isl_map_gist_basic_map(set_to_map(set),
3552 : bset_to_bmap(context)));
3553 : }
3554 :
3555 0 : __isl_give isl_set *isl_set_gist_params_basic_set(__isl_take isl_set *set,
3556 : __isl_take isl_basic_set *context)
3557 : {
3558 0 : isl_space *space = isl_set_get_space(set);
3559 0 : isl_basic_set *dom_context = isl_basic_set_universe(space);
3560 0 : dom_context = isl_basic_set_intersect_params(dom_context, context);
3561 0 : return isl_set_gist_basic_set(set, dom_context);
3562 : }
3563 :
3564 0 : __isl_give isl_set *isl_set_gist(__isl_take isl_set *set,
3565 : __isl_take isl_set *context)
3566 : {
3567 0 : return set_from_map(isl_map_gist(set_to_map(set), set_to_map(context)));
3568 : }
3569 :
3570 : /* Compute the gist of "bmap" with respect to the constraints "context"
3571 : * on the domain.
3572 : */
3573 0 : __isl_give isl_basic_map *isl_basic_map_gist_domain(
3574 : __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *context)
3575 : {
3576 0 : isl_space *space = isl_basic_map_get_space(bmap);
3577 0 : isl_basic_map *bmap_context = isl_basic_map_universe(space);
3578 :
3579 0 : bmap_context = isl_basic_map_intersect_domain(bmap_context, context);
3580 0 : return isl_basic_map_gist(bmap, bmap_context);
3581 : }
3582 :
3583 0 : __isl_give isl_map *isl_map_gist_domain(__isl_take isl_map *map,
3584 : __isl_take isl_set *context)
3585 : {
3586 0 : isl_map *map_context = isl_map_universe(isl_map_get_space(map));
3587 0 : map_context = isl_map_intersect_domain(map_context, context);
3588 0 : return isl_map_gist(map, map_context);
3589 : }
3590 :
3591 0 : __isl_give isl_map *isl_map_gist_range(__isl_take isl_map *map,
3592 : __isl_take isl_set *context)
3593 : {
3594 0 : isl_map *map_context = isl_map_universe(isl_map_get_space(map));
3595 0 : map_context = isl_map_intersect_range(map_context, context);
3596 0 : return isl_map_gist(map, map_context);
3597 : }
3598 :
3599 0 : __isl_give isl_map *isl_map_gist_params(__isl_take isl_map *map,
3600 : __isl_take isl_set *context)
3601 : {
3602 0 : isl_map *map_context = isl_map_universe(isl_map_get_space(map));
3603 0 : map_context = isl_map_intersect_params(map_context, context);
3604 0 : return isl_map_gist(map, map_context);
3605 : }
3606 :
3607 0 : __isl_give isl_set *isl_set_gist_params(__isl_take isl_set *set,
3608 : __isl_take isl_set *context)
3609 : {
3610 0 : return isl_map_gist_params(set, context);
3611 : }
3612 :
3613 : /* Quick check to see if two basic maps are disjoint.
3614 : * In particular, we reduce the equalities and inequalities of
3615 : * one basic map in the context of the equalities of the other
3616 : * basic map and check if we get a contradiction.
3617 : */
3618 7182202 : isl_bool isl_basic_map_plain_is_disjoint(__isl_keep isl_basic_map *bmap1,
3619 : __isl_keep isl_basic_map *bmap2)
3620 : {
3621 7182202 : struct isl_vec *v = NULL;
3622 7182202 : int *elim = NULL;
3623 : unsigned total;
3624 : int i;
3625 :
3626 7182202 : if (!bmap1 || !bmap2)
3627 0 : return isl_bool_error;
3628 7182202 : isl_assert(bmap1->ctx, isl_space_is_equal(bmap1->dim, bmap2->dim),
3629 : return isl_bool_error);
3630 7182202 : if (bmap1->n_div || bmap2->n_div)
3631 0 : return isl_bool_false;
3632 7182202 : if (!bmap1->n_eq && !bmap2->n_eq)
3633 3502132 : return isl_bool_false;
3634 :
3635 3680070 : total = isl_space_dim(bmap1->dim, isl_dim_all);
3636 3680070 : if (total == 0)
3637 0 : return isl_bool_false;
3638 3680070 : v = isl_vec_alloc(bmap1->ctx, 1 + total);
3639 3680070 : if (!v)
3640 0 : goto error;
3641 3680070 : elim = isl_alloc_array(bmap1->ctx, int, total);
3642 3680070 : if (!elim)
3643 0 : goto error;
3644 3680070 : compute_elimination_index(bmap1, elim);
3645 12086956 : for (i = 0; i < bmap2->n_eq; ++i) {
3646 : int reduced;
3647 8714484 : reduced = reduced_using_equalities(v->block.data, bmap2->eq[i],
3648 : bmap1, elim);
3649 9245802 : if (reduced && !isl_int_is_zero(v->block.data[0]) &&
3650 531318 : isl_seq_first_non_zero(v->block.data + 1, total) == -1)
3651 307598 : goto disjoint;
3652 : }
3653 32151434 : for (i = 0; i < bmap2->n_ineq; ++i) {
3654 : int reduced;
3655 29008702 : reduced = reduced_using_equalities(v->block.data,
3656 29008702 : bmap2->ineq[i], bmap1, elim);
3657 30278354 : if (reduced && isl_int_is_neg(v->block.data[0]) &&
3658 1269652 : isl_seq_first_non_zero(v->block.data + 1, total) == -1)
3659 229740 : goto disjoint;
3660 : }
3661 3142732 : compute_elimination_index(bmap2, elim);
3662 15599315 : for (i = 0; i < bmap1->n_ineq; ++i) {
3663 : int reduced;
3664 13988205 : reduced = reduced_using_equalities(v->block.data,
3665 13988205 : bmap1->ineq[i], bmap2, elim);
3666 16599944 : if (reduced && isl_int_is_neg(v->block.data[0]) &&
3667 2611739 : isl_seq_first_non_zero(v->block.data + 1, total) == -1)
3668 1531622 : goto disjoint;
3669 : }
3670 1611110 : isl_vec_free(v);
3671 1611110 : free(elim);
3672 1611110 : return isl_bool_false;
3673 : disjoint:
3674 2068960 : isl_vec_free(v);
3675 2068960 : free(elim);
3676 2068960 : return isl_bool_true;
3677 : error:
3678 0 : isl_vec_free(v);
3679 0 : free(elim);
3680 0 : return isl_bool_error;
3681 : }
3682 :
3683 0 : int isl_basic_set_plain_is_disjoint(__isl_keep isl_basic_set *bset1,
3684 : __isl_keep isl_basic_set *bset2)
3685 : {
3686 0 : return isl_basic_map_plain_is_disjoint(bset_to_bmap(bset1),
3687 : bset_to_bmap(bset2));
3688 : }
3689 :
3690 : /* Does "test" hold for all pairs of basic maps in "map1" and "map2"?
3691 : */
3692 328482 : static isl_bool all_pairs(__isl_keep isl_map *map1, __isl_keep isl_map *map2,
3693 : isl_bool (*test)(__isl_keep isl_basic_map *bmap1,
3694 : __isl_keep isl_basic_map *bmap2))
3695 : {
3696 : int i, j;
3697 :
3698 328482 : if (!map1 || !map2)
3699 0 : return isl_bool_error;
3700 :
3701 2974126 : for (i = 0; i < map1->n; ++i) {
3702 9827846 : for (j = 0; j < map2->n; ++j) {
3703 7182202 : isl_bool d = test(map1->p[i], map2->p[j]);
3704 7182202 : if (d != isl_bool_true)
3705 240268 : return d;
3706 : }
3707 : }
3708 :
3709 88214 : return isl_bool_true;
3710 : }
3711 :
3712 : /* Are "map1" and "map2" obviously disjoint, based on information
3713 : * that can be derived without looking at the individual basic maps?
3714 : *
3715 : * In particular, if one of them is empty or if they live in different spaces
3716 : * (ignoring parameters), then they are clearly disjoint.
3717 : */
3718 3023681 : static isl_bool isl_map_plain_is_disjoint_global(__isl_keep isl_map *map1,
3719 : __isl_keep isl_map *map2)
3720 : {
3721 : isl_bool disjoint;
3722 : isl_bool match;
3723 :
3724 3023681 : if (!map1 || !map2)
3725 0 : return isl_bool_error;
3726 :
3727 3023681 : disjoint = isl_map_plain_is_empty(map1);
3728 3023681 : if (disjoint < 0 || disjoint)
3729 87582 : return disjoint;
3730 :
3731 2936099 : disjoint = isl_map_plain_is_empty(map2);
3732 2936099 : if (disjoint < 0 || disjoint)
3733 656571 : return disjoint;
3734 :
3735 2279528 : match = isl_space_tuple_is_equal(map1->dim, isl_dim_in,
3736 : map2->dim, isl_dim_in);
3737 2279528 : if (match < 0 || !match)
3738 0 : return match < 0 ? isl_bool_error : isl_bool_true;
3739 :
3740 2279528 : match = isl_space_tuple_is_equal(map1->dim, isl_dim_out,
3741 : map2->dim, isl_dim_out);
3742 2279528 : if (match < 0 || !match)
3743 0 : return match < 0 ? isl_bool_error : isl_bool_true;
3744 :
3745 2279528 : return isl_bool_false;
3746 : }
3747 :
3748 : /* Are "map1" and "map2" obviously disjoint?
3749 : *
3750 : * If one of them is empty or if they live in different spaces (ignoring
3751 : * parameters), then they are clearly disjoint.
3752 : * This is checked by isl_map_plain_is_disjoint_global.
3753 : *
3754 : * If they have different parameters, then we skip any further tests.
3755 : *
3756 : * If they are obviously equal, but not obviously empty, then we will
3757 : * not be able to detect if they are disjoint.
3758 : *
3759 : * Otherwise we check if each basic map in "map1" is obviously disjoint
3760 : * from each basic map in "map2".
3761 : */
3762 0 : isl_bool isl_map_plain_is_disjoint(__isl_keep isl_map *map1,
3763 : __isl_keep isl_map *map2)
3764 : {
3765 : isl_bool disjoint;
3766 : isl_bool intersect;
3767 : isl_bool match;
3768 :
3769 0 : disjoint = isl_map_plain_is_disjoint_global(map1, map2);
3770 0 : if (disjoint < 0 || disjoint)
3771 0 : return disjoint;
3772 :
3773 0 : match = isl_map_has_equal_params(map1, map2);
3774 0 : if (match < 0 || !match)
3775 0 : return match < 0 ? isl_bool_error : isl_bool_false;
3776 :
3777 0 : intersect = isl_map_plain_is_equal(map1, map2);
3778 0 : if (intersect < 0 || intersect)
3779 0 : return intersect < 0 ? isl_bool_error : isl_bool_false;
3780 :
3781 0 : return all_pairs(map1, map2, &isl_basic_map_plain_is_disjoint);
3782 : }
3783 :
3784 : /* Are "map1" and "map2" disjoint?
3785 : * The parameters are assumed to have been aligned.
3786 : *
3787 : * In particular, check whether all pairs of basic maps are disjoint.
3788 : */
3789 328482 : static isl_bool isl_map_is_disjoint_aligned(__isl_keep isl_map *map1,
3790 : __isl_keep isl_map *map2)
3791 : {
3792 328482 : return all_pairs(map1, map2, &isl_basic_map_is_disjoint);
3793 : }
3794 :
3795 : /* Are "map1" and "map2" disjoint?
3796 : *
3797 : * They are disjoint if they are "obviously disjoint" or if one of them
3798 : * is empty. Otherwise, they are not disjoint if one of them is universal.
3799 : * If the two inputs are (obviously) equal and not empty, then they are
3800 : * not disjoint.
3801 : * If none of these cases apply, then check if all pairs of basic maps
3802 : * are disjoint after aligning the parameters.
3803 : */
3804 3023681 : isl_bool isl_map_is_disjoint(__isl_keep isl_map *map1, __isl_keep isl_map *map2)
3805 : {
3806 : isl_bool disjoint;
3807 : isl_bool intersect;
3808 :
3809 3023681 : disjoint = isl_map_plain_is_disjoint_global(map1, map2);
3810 3023681 : if (disjoint < 0 || disjoint)
3811 744153 : return disjoint;
3812 :
3813 2279528 : disjoint = isl_map_is_empty(map1);
3814 2279528 : if (disjoint < 0 || disjoint)
3815 0 : return disjoint;
3816 :
3817 2279528 : disjoint = isl_map_is_empty(map2);
3818 2279528 : if (disjoint < 0 || disjoint)
3819 0 : return disjoint;
3820 :
3821 2279528 : intersect = isl_map_plain_is_universe(map1);
3822 2279528 : if (intersect < 0 || intersect)
3823 1917279 : return intersect < 0 ? isl_bool_error : isl_bool_false;
3824 :
3825 362249 : intersect = isl_map_plain_is_universe(map2);
3826 362249 : if (intersect < 0 || intersect)
3827 33767 : return intersect < 0 ? isl_bool_error : isl_bool_false;
3828 :
3829 328482 : intersect = isl_map_plain_is_equal(map1, map2);
3830 328482 : if (intersect < 0 || intersect)
3831 0 : return isl_bool_not(intersect);
3832 :
3833 328482 : return isl_map_align_params_map_map_and_test(map1, map2,
3834 : &isl_map_is_disjoint_aligned);
3835 : }
3836 :
3837 : /* Are "bmap1" and "bmap2" disjoint?
3838 : *
3839 : * They are disjoint if they are "obviously disjoint" or if one of them
3840 : * is empty. Otherwise, they are not disjoint if one of them is universal.
3841 : * If none of these cases apply, we compute the intersection and see if
3842 : * the result is empty.
3843 : */
3844 7182202 : isl_bool isl_basic_map_is_disjoint(__isl_keep isl_basic_map *bmap1,
3845 : __isl_keep isl_basic_map *bmap2)
3846 : {
3847 : isl_bool disjoint;
3848 : isl_bool intersect;
3849 : isl_basic_map *test;
3850 :
3851 7182202 : disjoint = isl_basic_map_plain_is_disjoint(bmap1, bmap2);
3852 7182202 : if (disjoint < 0 || disjoint)
3853 2068960 : return disjoint;
3854 :
3855 5113242 : disjoint = isl_basic_map_is_empty(bmap1);
3856 5113242 : if (disjoint < 0 || disjoint)
3857 0 : return disjoint;
3858 :
3859 5113242 : disjoint = isl_basic_map_is_empty(bmap2);
3860 5113242 : if (disjoint < 0 || disjoint)
3861 0 : return disjoint;
3862 :
3863 5113242 : intersect = isl_basic_map_plain_is_universe(bmap1);
3864 5113242 : if (intersect < 0 || intersect)
3865 0 : return intersect < 0 ? isl_bool_error : isl_bool_false;
3866 :
3867 5113242 : intersect = isl_basic_map_plain_is_universe(bmap2);
3868 5113242 : if (intersect < 0 || intersect)
3869 0 : return intersect < 0 ? isl_bool_error : isl_bool_false;
3870 :
3871 5113242 : test = isl_basic_map_intersect(isl_basic_map_copy(bmap1),
3872 : isl_basic_map_copy(bmap2));
3873 5113242 : disjoint = isl_basic_map_is_empty(test);
3874 5113242 : isl_basic_map_free(test);
3875 :
3876 5113242 : return disjoint;
3877 : }
3878 :
3879 : /* Are "bset1" and "bset2" disjoint?
3880 : */
3881 0 : isl_bool isl_basic_set_is_disjoint(__isl_keep isl_basic_set *bset1,
3882 : __isl_keep isl_basic_set *bset2)
3883 : {
3884 0 : return isl_basic_map_is_disjoint(bset1, bset2);
3885 : }
3886 :
3887 0 : isl_bool isl_set_plain_is_disjoint(__isl_keep isl_set *set1,
3888 : __isl_keep isl_set *set2)
3889 : {
3890 0 : return isl_map_plain_is_disjoint(set_to_map(set1), set_to_map(set2));
3891 : }
3892 :
3893 : /* Are "set1" and "set2" disjoint?
3894 : */
3895 0 : isl_bool isl_set_is_disjoint(__isl_keep isl_set *set1, __isl_keep isl_set *set2)
3896 : {
3897 0 : return isl_map_is_disjoint(set1, set2);
3898 : }
3899 :
3900 : /* Is "v" equal to 0, 1 or -1?
3901 : */
3902 0 : static int is_zero_or_one(isl_int v)
3903 : {
3904 0 : return isl_int_is_zero(v) || isl_int_is_one(v) || isl_int_is_negone(v);
3905 : }
3906 :
3907 : /* Check if we can combine a given div with lower bound l and upper
3908 : * bound u with some other div and if so return that other div.
3909 : * Otherwise return -1.
3910 : *
3911 : * We first check that
3912 : * - the bounds are opposites of each other (except for the constant
3913 : * term)
3914 : * - the bounds do not reference any other div
3915 : * - no div is defined in terms of this div
3916 : *
3917 : * Let m be the size of the range allowed on the div by the bounds.
3918 : * That is, the bounds are of the form
3919 : *
3920 : * e <= a <= e + m - 1
3921 : *
3922 : * with e some expression in the other variables.
3923 : * We look for another div b such that no third div is defined in terms
3924 : * of this second div b and such that in any constraint that contains
3925 : * a (except for the given lower and upper bound), also contains b
3926 : * with a coefficient that is m times that of b.
3927 : * That is, all constraints (except for the lower and upper bound)
3928 : * are of the form
3929 : *
3930 : * e + f (a + m b) >= 0
3931 : *
3932 : * Furthermore, in the constraints that only contain b, the coefficient
3933 : * of b should be equal to 1 or -1.
3934 : * If so, we return b so that "a + m b" can be replaced by
3935 : * a single div "c = a + m b".
3936 : */
3937 0 : static int div_find_coalesce(struct isl_basic_map *bmap, int *pairs,
3938 : unsigned div, unsigned l, unsigned u)
3939 : {
3940 : int i, j;
3941 : unsigned dim;
3942 0 : int coalesce = -1;
3943 :
3944 0 : if (bmap->n_div <= 1)
3945 0 : return -1;
3946 0 : dim = isl_space_dim(bmap->dim, isl_dim_all);
3947 0 : if (isl_seq_first_non_zero(bmap->ineq[l] + 1 + dim, div) != -1)
3948 0 : return -1;
3949 0 : if (isl_seq_first_non_zero(bmap->ineq[l] + 1 + dim + div + 1,
3950 0 : bmap->n_div - div - 1) != -1)
3951 0 : return -1;
3952 0 : if (!isl_seq_is_neg(bmap->ineq[l] + 1, bmap->ineq[u] + 1,
3953 0 : dim + bmap->n_div))
3954 0 : return -1;
3955 :
3956 0 : for (i = 0; i < bmap->n_div; ++i) {
3957 0 : if (isl_int_is_zero(bmap->div[i][0]))
3958 0 : continue;
3959 0 : if (!isl_int_is_zero(bmap->div[i][1 + 1 + dim + div]))
3960 0 : return -1;
3961 : }
3962 :
3963 0 : isl_int_add(bmap->ineq[l][0], bmap->ineq[l][0], bmap->ineq[u][0]);
3964 0 : if (isl_int_is_neg(bmap->ineq[l][0])) {
3965 0 : isl_int_sub(bmap->ineq[l][0],
3966 : bmap->ineq[l][0], bmap->ineq[u][0]);
3967 0 : bmap = isl_basic_map_copy(bmap);
3968 0 : bmap = isl_basic_map_set_to_empty(bmap);
3969 0 : isl_basic_map_free(bmap);
3970 0 : return -1;
3971 : }
3972 0 : isl_int_add_ui(bmap->ineq[l][0], bmap->ineq[l][0], 1);
3973 0 : for (i = 0; i < bmap->n_div; ++i) {
3974 0 : if (i == div)
3975 0 : continue;
3976 0 : if (!pairs[i])
3977 0 : continue;
3978 0 : for (j = 0; j < bmap->n_div; ++j) {
3979 0 : if (isl_int_is_zero(bmap->div[j][0]))
3980 0 : continue;
3981 0 : if (!isl_int_is_zero(bmap->div[j][1 + 1 + dim + i]))
3982 : break;
3983 : }
3984 0 : if (j < bmap->n_div)
3985 0 : continue;
3986 0 : for (j = 0; j < bmap->n_ineq; ++j) {
3987 : int valid;
3988 0 : if (j == l || j == u)
3989 0 : continue;
3990 0 : if (isl_int_is_zero(bmap->ineq[j][1 + dim + div])) {
3991 0 : if (is_zero_or_one(bmap->ineq[j][1 + dim + i]))
3992 0 : continue;
3993 0 : break;
3994 : }
3995 0 : if (isl_int_is_zero(bmap->ineq[j][1 + dim + i]))
3996 0 : break;
3997 0 : isl_int_mul(bmap->ineq[j][1 + dim + div],
3998 : bmap->ineq[j][1 + dim + div],
3999 : bmap->ineq[l][0]);
4000 0 : valid = isl_int_eq(bmap->ineq[j][1 + dim + div],
4001 : bmap->ineq[j][1 + dim + i]);
4002 0 : isl_int_divexact(bmap->ineq[j][1 + dim + div],
4003 : bmap->ineq[j][1 + dim + div],
4004 : bmap->ineq[l][0]);
4005 0 : if (!valid)
4006 0 : break;
4007 : }
4008 0 : if (j < bmap->n_ineq)
4009 0 : continue;
4010 0 : coalesce = i;
4011 0 : break;
4012 : }
4013 0 : isl_int_sub_ui(bmap->ineq[l][0], bmap->ineq[l][0], 1);
4014 0 : isl_int_sub(bmap->ineq[l][0], bmap->ineq[l][0], bmap->ineq[u][0]);
4015 0 : return coalesce;
4016 : }
4017 :
4018 : /* Internal data structure used during the construction and/or evaluation of
4019 : * an inequality that ensures that a pair of bounds always allows
4020 : * for an integer value.
4021 : *
4022 : * "tab" is the tableau in which the inequality is evaluated. It may
4023 : * be NULL until it is actually needed.
4024 : * "v" contains the inequality coefficients.
4025 : * "g", "fl" and "fu" are temporary scalars used during the construction and
4026 : * evaluation.
4027 : */
4028 : struct test_ineq_data {
4029 : struct isl_tab *tab;
4030 : isl_vec *v;
4031 : isl_int g;
4032 : isl_int fl;
4033 : isl_int fu;
4034 : };
4035 :
4036 : /* Free all the memory allocated by the fields of "data".
4037 : */
4038 0 : static void test_ineq_data_clear(struct test_ineq_data *data)
4039 : {
4040 0 : isl_tab_free(data->tab);
4041 0 : isl_vec_free(data->v);
4042 0 : isl_int_clear(data->g);
4043 0 : isl_int_clear(data->fl);
4044 0 : isl_int_clear(data->fu);
4045 0 : }
4046 :
4047 : /* Is the inequality stored in data->v satisfied by "bmap"?
4048 : * That is, does it only attain non-negative values?
4049 : * data->tab is a tableau corresponding to "bmap".
4050 : */
4051 0 : static isl_bool test_ineq_is_satisfied(__isl_keep isl_basic_map *bmap,
4052 : struct test_ineq_data *data)
4053 : {
4054 : isl_ctx *ctx;
4055 : enum isl_lp_result res;
4056 :
4057 0 : ctx = isl_basic_map_get_ctx(bmap);
4058 0 : if (!data->tab)
4059 0 : data->tab = isl_tab_from_basic_map(bmap, 0);
4060 0 : res = isl_tab_min(data->tab, data->v->el, ctx->one, &data->g, NULL, 0);
4061 0 : if (res == isl_lp_error)
4062 0 : return isl_bool_error;
4063 0 : return res == isl_lp_ok && isl_int_is_nonneg(data->g);
4064 : }
4065 :
4066 : /* Given a lower and an upper bound on div i, do they always allow
4067 : * for an integer value of the given div?
4068 : * Determine this property by constructing an inequality
4069 : * such that the property is guaranteed when the inequality is nonnegative.
4070 : * The lower bound is inequality l, while the upper bound is inequality u.
4071 : * The constructed inequality is stored in data->v.
4072 : *
4073 : * Let the upper bound be
4074 : *
4075 : * -n_u a + e_u >= 0
4076 : *
4077 : * and the lower bound
4078 : *
4079 : * n_l a + e_l >= 0
4080 : *
4081 : * Let n_u = f_u g and n_l = f_l g, with g = gcd(n_u, n_l).
4082 : * We have
4083 : *
4084 : * - f_u e_l <= f_u f_l g a <= f_l e_u
4085 : *
4086 : * Since all variables are integer valued, this is equivalent to
4087 : *
4088 : * - f_u e_l - (f_u - 1) <= f_u f_l g a <= f_l e_u + (f_l - 1)
4089 : *
4090 : * If this interval is at least f_u f_l g, then it contains at least
4091 : * one integer value for a.
4092 : * That is, the test constraint is
4093 : *
4094 : * f_l e_u + f_u e_l + f_l - 1 + f_u - 1 + 1 >= f_u f_l g
4095 : *
4096 : * or
4097 : *
4098 : * f_l e_u + f_u e_l + f_l - 1 + f_u - 1 + 1 - f_u f_l g >= 0
4099 : *
4100 : * If the coefficients of f_l e_u + f_u e_l have a common divisor g',
4101 : * then the constraint can be scaled down by a factor g',
4102 : * with the constant term replaced by
4103 : * floor((f_l e_{u,0} + f_u e_{l,0} + f_l - 1 + f_u - 1 + 1 - f_u f_l g)/g').
4104 : * Note that the result of applying Fourier-Motzkin to this pair
4105 : * of constraints is
4106 : *
4107 : * f_l e_u + f_u e_l >= 0
4108 : *
4109 : * If the constant term of the scaled down version of this constraint,
4110 : * i.e., floor((f_l e_{u,0} + f_u e_{l,0})/g') is equal to the constant
4111 : * term of the scaled down test constraint, then the test constraint
4112 : * is known to hold and no explicit evaluation is required.
4113 : * This is essentially the Omega test.
4114 : *
4115 : * If the test constraint consists of only a constant term, then
4116 : * it is sufficient to look at the sign of this constant term.
4117 : */
4118 0 : static isl_bool int_between_bounds(__isl_keep isl_basic_map *bmap, int i,
4119 : int l, int u, struct test_ineq_data *data)
4120 : {
4121 : unsigned offset, n_div;
4122 0 : offset = isl_basic_map_offset(bmap, isl_dim_div);
4123 0 : n_div = isl_basic_map_dim(bmap, isl_dim_div);
4124 :
4125 0 : isl_int_gcd(data->g,
4126 : bmap->ineq[l][offset + i], bmap->ineq[u][offset + i]);
4127 0 : isl_int_divexact(data->fl, bmap->ineq[l][offset + i], data->g);
4128 0 : isl_int_divexact(data->fu, bmap->ineq[u][offset + i], data->g);
4129 0 : isl_int_neg(data->fu, data->fu);
4130 0 : isl_seq_combine(data->v->el, data->fl, bmap->ineq[u],
4131 0 : data->fu, bmap->ineq[l], offset + n_div);
4132 0 : isl_int_mul(data->g, data->g, data->fl);
4133 0 : isl_int_mul(data->g, data->g, data->fu);
4134 0 : isl_int_sub(data->g, data->g, data->fl);
4135 0 : isl_int_sub(data->g, data->g, data->fu);
4136 0 : isl_int_add_ui(data->g, data->g, 1);
4137 0 : isl_int_sub(data->fl, data->v->el[0], data->g);
4138 :
4139 0 : isl_seq_gcd(data->v->el + 1, offset - 1 + n_div, &data->g);
4140 0 : if (isl_int_is_zero(data->g))
4141 0 : return isl_int_is_nonneg(data->fl);
4142 0 : if (isl_int_is_one(data->g)) {
4143 0 : isl_int_set(data->v->el[0], data->fl);
4144 0 : return test_ineq_is_satisfied(bmap, data);
4145 : }
4146 0 : isl_int_fdiv_q(data->fl, data->fl, data->g);
4147 0 : isl_int_fdiv_q(data->v->el[0], data->v->el[0], data->g);
4148 0 : if (isl_int_eq(data->fl, data->v->el[0]))
4149 0 : return isl_bool_true;
4150 0 : isl_int_set(data->v->el[0], data->fl);
4151 0 : isl_seq_scale_down(data->v->el + 1, data->v->el + 1, data->g,
4152 0 : offset - 1 + n_div);
4153 :
4154 0 : return test_ineq_is_satisfied(bmap, data);
4155 : }
4156 :
4157 : /* Remove more kinds of divs that are not strictly needed.
4158 : * In particular, if all pairs of lower and upper bounds on a div
4159 : * are such that they allow at least one integer value of the div,
4160 : * then we can eliminate the div using Fourier-Motzkin without
4161 : * introducing any spurious solutions.
4162 : *
4163 : * If at least one of the two constraints has a unit coefficient for the div,
4164 : * then the presence of such a value is guaranteed so there is no need to check.
4165 : * In particular, the value attained by the bound with unit coefficient
4166 : * can serve as this intermediate value.
4167 : */
4168 0 : static __isl_give isl_basic_map *drop_more_redundant_divs(
4169 : __isl_take isl_basic_map *bmap, __isl_take int *pairs, int n)
4170 : {
4171 : isl_ctx *ctx;
4172 0 : struct test_ineq_data data = { NULL, NULL };
4173 : unsigned off, n_div;
4174 0 : int remove = -1;
4175 :
4176 0 : isl_int_init(data.g);
4177 0 : isl_int_init(data.fl);
4178 0 : isl_int_init(data.fu);
4179 :
4180 0 : if (!bmap)
4181 0 : goto error;
4182 :
4183 0 : ctx = isl_basic_map_get_ctx(bmap);
4184 0 : off = isl_basic_map_offset(bmap, isl_dim_div);
4185 0 : n_div = isl_basic_map_dim(bmap, isl_dim_div);
4186 0 : data.v = isl_vec_alloc(ctx, off + n_div);
4187 0 : if (!data.v)
4188 0 : goto error;
4189 :
4190 0 : while (n > 0) {
4191 : int i, l, u;
4192 0 : int best = -1;
4193 : isl_bool has_int;
4194 :
4195 0 : for (i = 0; i < n_div; ++i) {
4196 0 : if (!pairs[i])
4197 0 : continue;
4198 0 : if (best >= 0 && pairs[best] <= pairs[i])
4199 0 : continue;
4200 0 : best = i;
4201 : }
4202 :
4203 0 : i = best;
4204 0 : for (l = 0; l < bmap->n_ineq; ++l) {
4205 0 : if (!isl_int_is_pos(bmap->ineq[l][off + i]))
4206 0 : continue;
4207 0 : if (isl_int_is_one(bmap->ineq[l][off + i]))
4208 0 : continue;
4209 0 : for (u = 0; u < bmap->n_ineq; ++u) {
4210 0 : if (!isl_int_is_neg(bmap->ineq[u][off + i]))
4211 0 : continue;
4212 0 : if (isl_int_is_negone(bmap->ineq[u][off + i]))
4213 0 : continue;
4214 0 : has_int = int_between_bounds(bmap, i, l, u,
4215 : &data);
4216 0 : if (has_int < 0)
4217 0 : goto error;
4218 0 : if (data.tab && data.tab->empty)
4219 0 : break;
4220 0 : if (!has_int)
4221 0 : break;
4222 : }
4223 0 : if (u < bmap->n_ineq)
4224 0 : break;
4225 : }
4226 0 : if (data.tab && data.tab->empty) {
4227 0 : bmap = isl_basic_map_set_to_empty(bmap);
4228 0 : break;
4229 : }
4230 0 : if (l == bmap->n_ineq) {
4231 0 : remove = i;
4232 0 : break;
4233 : }
4234 0 : pairs[i] = 0;
4235 0 : --n;
4236 : }
4237 :
4238 0 : test_ineq_data_clear(&data);
4239 :
4240 0 : free(pairs);
4241 :
4242 0 : if (remove < 0)
4243 0 : return bmap;
4244 :
4245 0 : bmap = isl_basic_map_remove_dims(bmap, isl_dim_div, remove, 1);
4246 0 : return isl_basic_map_drop_redundant_divs(bmap);
4247 : error:
4248 0 : free(pairs);
4249 0 : isl_basic_map_free(bmap);
4250 0 : test_ineq_data_clear(&data);
4251 0 : return NULL;
4252 : }
4253 :
4254 : /* Given a pair of divs div1 and div2 such that, except for the lower bound l
4255 : * and the upper bound u, div1 always occurs together with div2 in the form
4256 : * (div1 + m div2), where m is the constant range on the variable div1
4257 : * allowed by l and u, replace the pair div1 and div2 by a single
4258 : * div that is equal to div1 + m div2.
4259 : *
4260 : * The new div will appear in the location that contains div2.
4261 : * We need to modify all constraints that contain
4262 : * div2 = (div - div1) / m
4263 : * The coefficient of div2 is known to be equal to 1 or -1.
4264 : * (If a constraint does not contain div2, it will also not contain div1.)
4265 : * If the constraint also contains div1, then we know they appear
4266 : * as f (div1 + m div2) and we can simply replace (div1 + m div2) by div,
4267 : * i.e., the coefficient of div is f.
4268 : *
4269 : * Otherwise, we first need to introduce div1 into the constraint.
4270 : * Let l be
4271 : *
4272 : * div1 + f >=0
4273 : *
4274 : * and u
4275 : *
4276 : * -div1 + f' >= 0
4277 : *
4278 : * A lower bound on div2
4279 : *
4280 : * div2 + t >= 0
4281 : *
4282 : * can be replaced by
4283 : *
4284 : * m div2 + div1 + m t + f >= 0
4285 : *
4286 : * An upper bound
4287 : *
4288 : * -div2 + t >= 0
4289 : *
4290 : * can be replaced by
4291 : *
4292 : * -(m div2 + div1) + m t + f' >= 0
4293 : *
4294 : * These constraint are those that we would obtain from eliminating
4295 : * div1 using Fourier-Motzkin.
4296 : *
4297 : * After all constraints have been modified, we drop the lower and upper
4298 : * bound and then drop div1.
4299 : * Since the new div is only placed in the same location that used
4300 : * to store div2, but otherwise has a different meaning, any possible
4301 : * explicit representation of the original div2 is removed.
4302 : */
4303 0 : static __isl_give isl_basic_map *coalesce_divs(__isl_take isl_basic_map *bmap,
4304 : unsigned div1, unsigned div2, unsigned l, unsigned u)
4305 : {
4306 : isl_ctx *ctx;
4307 : isl_int m;
4308 : unsigned dim, total;
4309 : int i;
4310 :
4311 0 : ctx = isl_basic_map_get_ctx(bmap);
4312 :
4313 0 : dim = isl_space_dim(bmap->dim, isl_dim_all);
4314 0 : total = 1 + dim + bmap->n_div;
4315 :
4316 0 : isl_int_init(m);
4317 0 : isl_int_add(m, bmap->ineq[l][0], bmap->ineq[u][0]);
4318 0 : isl_int_add_ui(m, m, 1);
4319 :
4320 0 : for (i = 0; i < bmap->n_ineq; ++i) {
4321 0 : if (i == l || i == u)
4322 0 : continue;
4323 0 : if (isl_int_is_zero(bmap->ineq[i][1 + dim + div2]))
4324 0 : continue;
4325 0 : if (isl_int_is_zero(bmap->ineq[i][1 + dim + div1])) {
4326 0 : if (isl_int_is_pos(bmap->ineq[i][1 + dim + div2]))
4327 0 : isl_seq_combine(bmap->ineq[i], m, bmap->ineq[i],
4328 0 : ctx->one, bmap->ineq[l], total);
4329 : else
4330 0 : isl_seq_combine(bmap->ineq[i], m, bmap->ineq[i],
4331 0 : ctx->one, bmap->ineq[u], total);
4332 : }
4333 0 : isl_int_set(bmap->ineq[i][1 + dim + div2],
4334 : bmap->ineq[i][1 + dim + div1]);
4335 0 : isl_int_set_si(bmap->ineq[i][1 + dim + div1], 0);
4336 : }
4337 :
4338 0 : isl_int_clear(m);
4339 0 : if (l > u) {
4340 0 : isl_basic_map_drop_inequality(bmap, l);
4341 0 : isl_basic_map_drop_inequality(bmap, u);
4342 : } else {
4343 0 : isl_basic_map_drop_inequality(bmap, u);
4344 0 : isl_basic_map_drop_inequality(bmap, l);
4345 : }
4346 0 : bmap = isl_basic_map_mark_div_unknown(bmap, div2);
4347 0 : bmap = isl_basic_map_drop_div(bmap, div1);
4348 0 : return bmap;
4349 : }
4350 :
4351 : /* First check if we can coalesce any pair of divs and
4352 : * then continue with dropping more redundant divs.
4353 : *
4354 : * We loop over all pairs of lower and upper bounds on a div
4355 : * with coefficient 1 and -1, respectively, check if there
4356 : * is any other div "c" with which we can coalesce the div
4357 : * and if so, perform the coalescing.
4358 : */
4359 0 : static __isl_give isl_basic_map *coalesce_or_drop_more_redundant_divs(
4360 : __isl_take isl_basic_map *bmap, int *pairs, int n)
4361 : {
4362 : int i, l, u;
4363 : unsigned dim;
4364 :
4365 0 : dim = isl_space_dim(bmap->dim, isl_dim_all);
4366 :
4367 0 : for (i = 0; i < bmap->n_div; ++i) {
4368 0 : if (!pairs[i])
4369 0 : continue;
4370 0 : for (l = 0; l < bmap->n_ineq; ++l) {
4371 0 : if (!isl_int_is_one(bmap->ineq[l][1 + dim + i]))
4372 0 : continue;
4373 0 : for (u = 0; u < bmap->n_ineq; ++u) {
4374 : int c;
4375 :
4376 0 : if (!isl_int_is_negone(bmap->ineq[u][1+dim+i]))
4377 0 : continue;
4378 0 : c = div_find_coalesce(bmap, pairs, i, l, u);
4379 0 : if (c < 0)
4380 0 : continue;
4381 0 : free(pairs);
4382 0 : bmap = coalesce_divs(bmap, i, c, l, u);
4383 0 : return isl_basic_map_drop_redundant_divs(bmap);
4384 : }
4385 : }
4386 : }
4387 :
4388 0 : if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY)) {
4389 0 : free(pairs);
4390 0 : return bmap;
4391 : }
4392 :
4393 0 : return drop_more_redundant_divs(bmap, pairs, n);
4394 : }
4395 :
4396 : /* Are the "n" coefficients starting at "first" of inequality constraints
4397 : * "i" and "j" of "bmap" equal to each other?
4398 : */
4399 0 : static int is_parallel_part(__isl_keep isl_basic_map *bmap, int i, int j,
4400 : int first, int n)
4401 : {
4402 0 : return isl_seq_eq(bmap->ineq[i] + first, bmap->ineq[j] + first, n);
4403 : }
4404 :
4405 : /* Are the "n" coefficients starting at "first" of inequality constraints
4406 : * "i" and "j" of "bmap" opposite to each other?
4407 : */
4408 0 : static int is_opposite_part(__isl_keep isl_basic_map *bmap, int i, int j,
4409 : int first, int n)
4410 : {
4411 0 : return isl_seq_is_neg(bmap->ineq[i] + first, bmap->ineq[j] + first, n);
4412 : }
4413 :
4414 : /* Are inequality constraints "i" and "j" of "bmap" opposite to each other,
4415 : * apart from the constant term?
4416 : */
4417 0 : static isl_bool is_opposite(__isl_keep isl_basic_map *bmap, int i, int j)
4418 : {
4419 : unsigned total;
4420 :
4421 0 : total = isl_basic_map_dim(bmap, isl_dim_all);
4422 0 : return is_opposite_part(bmap, i, j, 1, total);
4423 : }
4424 :
4425 : /* Are inequality constraints "i" and "j" of "bmap" equal to each other,
4426 : * apart from the constant term and the coefficient at position "pos"?
4427 : */
4428 0 : static int is_parallel_except(__isl_keep isl_basic_map *bmap, int i, int j,
4429 : int pos)
4430 : {
4431 : unsigned total;
4432 :
4433 0 : total = isl_basic_map_dim(bmap, isl_dim_all);
4434 0 : return is_parallel_part(bmap, i, j, 1, pos - 1) &&
4435 0 : is_parallel_part(bmap, i, j, pos + 1, total - pos);
4436 : }
4437 :
4438 : /* Are inequality constraints "i" and "j" of "bmap" opposite to each other,
4439 : * apart from the constant term and the coefficient at position "pos"?
4440 : */
4441 0 : static int is_opposite_except(__isl_keep isl_basic_map *bmap, int i, int j,
4442 : int pos)
4443 : {
4444 : unsigned total;
4445 :
4446 0 : total = isl_basic_map_dim(bmap, isl_dim_all);
4447 0 : return is_opposite_part(bmap, i, j, 1, pos - 1) &&
4448 0 : is_opposite_part(bmap, i, j, pos + 1, total - pos);
4449 : }
4450 :
4451 : /* Restart isl_basic_map_drop_redundant_divs after "bmap" has
4452 : * been modified, simplying it if "simplify" is set.
4453 : * Free the temporary data structure "pairs" that was associated
4454 : * to the old version of "bmap".
4455 : */
4456 0 : static __isl_give isl_basic_map *drop_redundant_divs_again(
4457 : __isl_take isl_basic_map *bmap, __isl_take int *pairs, int simplify)
4458 : {
4459 0 : if (simplify)
4460 0 : bmap = isl_basic_map_simplify(bmap);
4461 0 : free(pairs);
4462 0 : return isl_basic_map_drop_redundant_divs(bmap);
4463 : }
4464 :
4465 : /* Is "div" the single unknown existentially quantified variable
4466 : * in inequality constraint "ineq" of "bmap"?
4467 : * "div" is known to have a non-zero coefficient in "ineq".
4468 : */
4469 0 : static isl_bool single_unknown(__isl_keep isl_basic_map *bmap, int ineq,
4470 : int div)
4471 : {
4472 : int i;
4473 : unsigned n_div, o_div;
4474 : isl_bool known;
4475 :
4476 0 : known = isl_basic_map_div_is_known(bmap, div);
4477 0 : if (known < 0 || known)
4478 0 : return isl_bool_not(known);
4479 0 : n_div = isl_basic_map_dim(bmap, isl_dim_div);
4480 0 : if (n_div == 1)
4481 0 : return isl_bool_true;
4482 0 : o_div = isl_basic_map_offset(bmap, isl_dim_div);
4483 0 : for (i = 0; i < n_div; ++i) {
4484 : isl_bool known;
4485 :
4486 0 : if (i == div)
4487 0 : continue;
4488 0 : if (isl_int_is_zero(bmap->ineq[ineq][o_div + i]))
4489 0 : continue;
4490 0 : known = isl_basic_map_div_is_known(bmap, i);
4491 0 : if (known < 0 || !known)
4492 0 : return known;
4493 : }
4494 :
4495 0 : return isl_bool_true;
4496 : }
4497 :
4498 : /* Does integer division "div" have coefficient 1 in inequality constraint
4499 : * "ineq" of "map"?
4500 : */
4501 0 : static isl_bool has_coef_one(__isl_keep isl_basic_map *bmap, int div, int ineq)
4502 : {
4503 : unsigned o_div;
4504 :
4505 0 : o_div = isl_basic_map_offset(bmap, isl_dim_div);
4506 0 : if (isl_int_is_one(bmap->ineq[ineq][o_div + div]))
4507 0 : return isl_bool_true;
4508 :
4509 0 : return isl_bool_false;
4510 : }
4511 :
4512 : /* Turn inequality constraint "ineq" of "bmap" into an equality and
4513 : * then try and drop redundant divs again,
4514 : * freeing the temporary data structure "pairs" that was associated
4515 : * to the old version of "bmap".
4516 : */
4517 0 : static __isl_give isl_basic_map *set_eq_and_try_again(
4518 : __isl_take isl_basic_map *bmap, int ineq, __isl_take int *pairs)
4519 : {
4520 0 : bmap = isl_basic_map_cow(bmap);
4521 0 : isl_basic_map_inequality_to_equality(bmap, ineq);
4522 0 : return drop_redundant_divs_again(bmap, pairs, 1);
4523 : }
4524 :
4525 : /* Drop the integer division at position "div", along with the two
4526 : * inequality constraints "ineq1" and "ineq2" in which it appears
4527 : * from "bmap" and then try and drop redundant divs again,
4528 : * freeing the temporary data structure "pairs" that was associated
4529 : * to the old version of "bmap".
4530 : */
4531 0 : static __isl_give isl_basic_map *drop_div_and_try_again(
4532 : __isl_take isl_basic_map *bmap, int div, int ineq1, int ineq2,
4533 : __isl_take int *pairs)
4534 : {
4535 0 : if (ineq1 > ineq2) {
4536 0 : isl_basic_map_drop_inequality(bmap, ineq1);
4537 0 : isl_basic_map_drop_inequality(bmap, ineq2);
4538 : } else {
4539 0 : isl_basic_map_drop_inequality(bmap, ineq2);
4540 0 : isl_basic_map_drop_inequality(bmap, ineq1);
4541 : }
4542 0 : bmap = isl_basic_map_drop_div(bmap, div);
4543 0 : return drop_redundant_divs_again(bmap, pairs, 0);
4544 : }
4545 :
4546 : /* Given two inequality constraints
4547 : *
4548 : * f(x) + n d + c >= 0, (ineq)
4549 : *
4550 : * with d the variable at position "pos", and
4551 : *
4552 : * f(x) + c0 >= 0, (lower)
4553 : *
4554 : * compute the maximal value of the lower bound ceil((-f(x) - c)/n)
4555 : * determined by the first constraint.
4556 : * That is, store
4557 : *
4558 : * ceil((c0 - c)/n)
4559 : *
4560 : * in *l.
4561 : */
4562 0 : static void lower_bound_from_parallel(__isl_keep isl_basic_map *bmap,
4563 : int ineq, int lower, int pos, isl_int *l)
4564 : {
4565 0 : isl_int_neg(*l, bmap->ineq[ineq][0]);
4566 0 : isl_int_add(*l, *l, bmap->ineq[lower][0]);
4567 0 : isl_int_cdiv_q(*l, *l, bmap->ineq[ineq][pos]);
4568 0 : }
4569 :
4570 : /* Given two inequality constraints
4571 : *
4572 : * f(x) + n d + c >= 0, (ineq)
4573 : *
4574 : * with d the variable at position "pos", and
4575 : *
4576 : * -f(x) - c0 >= 0, (upper)
4577 : *
4578 : * compute the minimal value of the lower bound ceil((-f(x) - c)/n)
4579 : * determined by the first constraint.
4580 : * That is, store
4581 : *
4582 : * ceil((-c1 - c)/n)
4583 : *
4584 : * in *u.
4585 : */
4586 0 : static void lower_bound_from_opposite(__isl_keep isl_basic_map *bmap,
4587 : int ineq, int upper, int pos, isl_int *u)
4588 : {
4589 0 : isl_int_neg(*u, bmap->ineq[ineq][0]);
4590 0 : isl_int_sub(*u, *u, bmap->ineq[upper][0]);
4591 0 : isl_int_cdiv_q(*u, *u, bmap->ineq[ineq][pos]);
4592 0 : }
4593 :
4594 : /* Given a lower bound constraint "ineq" on "div" in "bmap",
4595 : * does the corresponding lower bound have a fixed value in "bmap"?
4596 : *
4597 : * In particular, "ineq" is of the form
4598 : *
4599 : * f(x) + n d + c >= 0
4600 : *
4601 : * with n > 0, c the constant term and
4602 : * d the existentially quantified variable "div".
4603 : * That is, the lower bound is
4604 : *
4605 : * ceil((-f(x) - c)/n)
4606 : *
4607 : * Look for a pair of constraints
4608 : *
4609 : * f(x) + c0 >= 0
4610 : * -f(x) + c1 >= 0
4611 : *
4612 : * i.e., -c1 <= -f(x) <= c0, that fix ceil((-f(x) - c)/n) to a constant value.
4613 : * That is, check that
4614 : *
4615 : * ceil((-c1 - c)/n) = ceil((c0 - c)/n)
4616 : *
4617 : * If so, return the index of inequality f(x) + c0 >= 0.
4618 : * Otherwise, return -1.
4619 : */
4620 0 : static int lower_bound_is_cst(__isl_keep isl_basic_map *bmap, int div, int ineq)
4621 : {
4622 : int i;
4623 0 : int lower = -1, upper = -1;
4624 : unsigned o_div;
4625 : isl_int l, u;
4626 : int equal;
4627 :
4628 0 : o_div = isl_basic_map_offset(bmap, isl_dim_div);
4629 0 : for (i = 0; i < bmap->n_ineq && (lower < 0 || upper < 0); ++i) {
4630 0 : if (i == ineq)
4631 0 : continue;
4632 0 : if (!isl_int_is_zero(bmap->ineq[i][o_div + div]))
4633 0 : continue;
4634 0 : if (lower < 0 &&
4635 0 : is_parallel_except(bmap, ineq, i, o_div + div)) {
4636 0 : lower = i;
4637 0 : continue;
4638 : }
4639 0 : if (upper < 0 &&
4640 0 : is_opposite_except(bmap, ineq, i, o_div + div)) {
4641 0 : upper = i;
4642 : }
4643 : }
4644 :
4645 0 : if (lower < 0 || upper < 0)
4646 0 : return -1;
4647 :
4648 0 : isl_int_init(l);
4649 0 : isl_int_init(u);
4650 :
4651 0 : lower_bound_from_parallel(bmap, ineq, lower, o_div + div, &l);
4652 0 : lower_bound_from_opposite(bmap, ineq, upper, o_div + div, &u);
4653 :
4654 0 : equal = isl_int_eq(l, u);
4655 :
4656 0 : isl_int_clear(l);
4657 0 : isl_int_clear(u);
4658 :
4659 0 : return equal ? lower : -1;
4660 : }
4661 :
4662 : /* Given a lower bound constraint "ineq" on the existentially quantified
4663 : * variable "div", such that the corresponding lower bound has
4664 : * a fixed value in "bmap", assign this fixed value to the variable and
4665 : * then try and drop redundant divs again,
4666 : * freeing the temporary data structure "pairs" that was associated
4667 : * to the old version of "bmap".
4668 : * "lower" determines the constant value for the lower bound.
4669 : *
4670 : * In particular, "ineq" is of the form
4671 : *
4672 : * f(x) + n d + c >= 0,
4673 : *
4674 : * while "lower" is of the form
4675 : *
4676 : * f(x) + c0 >= 0
4677 : *
4678 : * The lower bound is ceil((-f(x) - c)/n) and its constant value
4679 : * is ceil((c0 - c)/n).
4680 : */
4681 0 : static __isl_give isl_basic_map *fix_cst_lower(__isl_take isl_basic_map *bmap,
4682 : int div, int ineq, int lower, int *pairs)
4683 : {
4684 : isl_int c;
4685 : unsigned o_div;
4686 :
4687 0 : isl_int_init(c);
4688 :
4689 0 : o_div = isl_basic_map_offset(bmap, isl_dim_div);
4690 0 : lower_bound_from_parallel(bmap, ineq, lower, o_div + div, &c);
4691 0 : bmap = isl_basic_map_fix(bmap, isl_dim_div, div, c);
4692 0 : free(pairs);
4693 :
4694 0 : isl_int_clear(c);
4695 :
4696 0 : return isl_basic_map_drop_redundant_divs(bmap);
4697 : }
4698 :
4699 : /* Remove divs that are not strictly needed based on the inequality
4700 : * constraints.
4701 : * In particular, if a div only occurs positively (or negatively)
4702 : * in constraints, then it can simply be dropped.
4703 : * Also, if a div occurs in only two constraints and if moreover
4704 : * those two constraints are opposite to each other, except for the constant
4705 : * term and if the sum of the constant terms is such that for any value
4706 : * of the other values, there is always at least one integer value of the
4707 : * div, i.e., if one plus this sum is greater than or equal to
4708 : * the (absolute value) of the coefficient of the div in the constraints,
4709 : * then we can also simply drop the div.
4710 : *
4711 : * If an existentially quantified variable does not have an explicit
4712 : * representation, appears in only a single lower bound that does not
4713 : * involve any other such existentially quantified variables and appears
4714 : * in this lower bound with coefficient 1,
4715 : * then fix the variable to the value of the lower bound. That is,
4716 : * turn the inequality into an equality.
4717 : * If for any value of the other variables, there is any value
4718 : * for the existentially quantified variable satisfying the constraints,
4719 : * then this lower bound also satisfies the constraints.
4720 : * It is therefore safe to pick this lower bound.
4721 : *
4722 : * The same reasoning holds even if the coefficient is not one.
4723 : * However, fixing the variable to the value of the lower bound may
4724 : * in general introduce an extra integer division, in which case
4725 : * it may be better to pick another value.
4726 : * If this integer division has a known constant value, then plugging
4727 : * in this constant value removes the existentially quantified variable
4728 : * completely. In particular, if the lower bound is of the form
4729 : * ceil((-f(x) - c)/n) and there are two constraints, f(x) + c0 >= 0 and
4730 : * -f(x) + c1 >= 0 such that ceil((-c1 - c)/n) = ceil((c0 - c)/n),
4731 : * then the existentially quantified variable can be assigned this
4732 : * shared value.
4733 : *
4734 : * We skip divs that appear in equalities or in the definition of other divs.
4735 : * Divs that appear in the definition of other divs usually occur in at least
4736 : * 4 constraints, but the constraints may have been simplified.
4737 : *
4738 : * If any divs are left after these simple checks then we move on
4739 : * to more complicated cases in drop_more_redundant_divs.
4740 : */
4741 0 : static __isl_give isl_basic_map *isl_basic_map_drop_redundant_divs_ineq(
4742 : __isl_take isl_basic_map *bmap)
4743 : {
4744 : int i, j;
4745 : unsigned off;
4746 0 : int *pairs = NULL;
4747 0 : int n = 0;
4748 :
4749 0 : if (!bmap)
4750 0 : goto error;
4751 0 : if (bmap->n_div == 0)
4752 0 : return bmap;
4753 :
4754 0 : off = isl_space_dim(bmap->dim, isl_dim_all);
4755 0 : pairs = isl_calloc_array(bmap->ctx, int, bmap->n_div);
4756 0 : if (!pairs)
4757 0 : goto error;
4758 :
4759 0 : for (i = 0; i < bmap->n_div; ++i) {
4760 : int pos, neg;
4761 : int last_pos, last_neg;
4762 : int redundant;
4763 : int defined;
4764 : isl_bool opp, set_div;
4765 :
4766 0 : defined = !isl_int_is_zero(bmap->div[i][0]);
4767 0 : for (j = i; j < bmap->n_div; ++j)
4768 0 : if (!isl_int_is_zero(bmap->div[j][1 + 1 + off + i]))
4769 : break;
4770 0 : if (j < bmap->n_div)
4771 0 : continue;
4772 0 : for (j = 0; j < bmap->n_eq; ++j)
4773 0 : if (!isl_int_is_zero(bmap->eq[j][1 + off + i]))
4774 : break;
4775 0 : if (j < bmap->n_eq)
4776 0 : continue;
4777 0 : ++n;
4778 0 : pos = neg = 0;
4779 0 : for (j = 0; j < bmap->n_ineq; ++j) {
4780 0 : if (isl_int_is_pos(bmap->ineq[j][1 + off + i])) {
4781 0 : last_pos = j;
4782 0 : ++pos;
4783 : }
4784 0 : if (isl_int_is_neg(bmap->ineq[j][1 + off + i])) {
4785 0 : last_neg = j;
4786 0 : ++neg;
4787 : }
4788 : }
4789 0 : pairs[i] = pos * neg;
4790 0 : if (pairs[i] == 0) {
4791 0 : for (j = bmap->n_ineq - 1; j >= 0; --j)
4792 0 : if (!isl_int_is_zero(bmap->ineq[j][1+off+i]))
4793 0 : isl_basic_map_drop_inequality(bmap, j);
4794 0 : bmap = isl_basic_map_drop_div(bmap, i);
4795 0 : return drop_redundant_divs_again(bmap, pairs, 0);
4796 : }
4797 0 : if (pairs[i] != 1)
4798 0 : opp = isl_bool_false;
4799 : else
4800 0 : opp = is_opposite(bmap, last_pos, last_neg);
4801 0 : if (opp < 0)
4802 0 : goto error;
4803 0 : if (!opp) {
4804 : int lower;
4805 : isl_bool single, one;
4806 :
4807 0 : if (pos != 1)
4808 0 : continue;
4809 0 : single = single_unknown(bmap, last_pos, i);
4810 0 : if (single < 0)
4811 0 : goto error;
4812 0 : if (!single)
4813 0 : continue;
4814 0 : one = has_coef_one(bmap, i, last_pos);
4815 0 : if (one < 0)
4816 0 : goto error;
4817 0 : if (one)
4818 0 : return set_eq_and_try_again(bmap, last_pos,
4819 : pairs);
4820 0 : lower = lower_bound_is_cst(bmap, i, last_pos);
4821 0 : if (lower >= 0)
4822 0 : return fix_cst_lower(bmap, i, last_pos, lower,
4823 : pairs);
4824 0 : continue;
4825 : }
4826 :
4827 0 : isl_int_add(bmap->ineq[last_pos][0],
4828 : bmap->ineq[last_pos][0], bmap->ineq[last_neg][0]);
4829 0 : isl_int_add_ui(bmap->ineq[last_pos][0],
4830 : bmap->ineq[last_pos][0], 1);
4831 0 : redundant = isl_int_ge(bmap->ineq[last_pos][0],
4832 : bmap->ineq[last_pos][1+off+i]);
4833 0 : isl_int_sub_ui(bmap->ineq[last_pos][0],
4834 : bmap->ineq[last_pos][0], 1);
4835 0 : isl_int_sub(bmap->ineq[last_pos][0],
4836 : bmap->ineq[last_pos][0], bmap->ineq[last_neg][0]);
4837 0 : if (redundant)
4838 0 : return drop_div_and_try_again(bmap, i,
4839 : last_pos, last_neg, pairs);
4840 0 : if (defined)
4841 0 : set_div = isl_bool_false;
4842 : else
4843 0 : set_div = ok_to_set_div_from_bound(bmap, i, last_pos);
4844 0 : if (set_div < 0)
4845 0 : return isl_basic_map_free(bmap);
4846 0 : if (set_div) {
4847 0 : bmap = set_div_from_lower_bound(bmap, i, last_pos);
4848 0 : return drop_redundant_divs_again(bmap, pairs, 1);
4849 : }
4850 0 : pairs[i] = 0;
4851 0 : --n;
4852 : }
4853 :
4854 0 : if (n > 0)
4855 0 : return coalesce_or_drop_more_redundant_divs(bmap, pairs, n);
4856 :
4857 0 : free(pairs);
4858 0 : return bmap;
4859 : error:
4860 0 : free(pairs);
4861 0 : isl_basic_map_free(bmap);
4862 0 : return NULL;
4863 : }
4864 :
4865 : /* Consider the coefficients at "c" as a row vector and replace
4866 : * them with their product with "T". "T" is assumed to be a square matrix.
4867 : */
4868 0 : static isl_stat preimage(isl_int *c, __isl_keep isl_mat *T)
4869 : {
4870 : int n;
4871 : isl_ctx *ctx;
4872 : isl_vec *v;
4873 :
4874 0 : if (!T)
4875 0 : return isl_stat_error;
4876 0 : n = isl_mat_rows(T);
4877 0 : if (isl_seq_first_non_zero(c, n) == -1)
4878 0 : return isl_stat_ok;
4879 0 : ctx = isl_mat_get_ctx(T);
4880 0 : v = isl_vec_alloc(ctx, n);
4881 0 : if (!v)
4882 0 : return isl_stat_error;
4883 0 : isl_seq_swp_or_cpy(v->el, c, n);
4884 0 : v = isl_vec_mat_product(v, isl_mat_copy(T));
4885 0 : if (!v)
4886 0 : return isl_stat_error;
4887 0 : isl_seq_swp_or_cpy(c, v->el, n);
4888 0 : isl_vec_free(v);
4889 :
4890 0 : return isl_stat_ok;
4891 : }
4892 :
4893 : /* Plug in T for the variables in "bmap" starting at "pos".
4894 : * T is a linear unimodular matrix, i.e., without constant term.
4895 : */
4896 0 : static __isl_give isl_basic_map *isl_basic_map_preimage_vars(
4897 : __isl_take isl_basic_map *bmap, unsigned pos, __isl_take isl_mat *T)
4898 : {
4899 : int i;
4900 : unsigned n, total;
4901 :
4902 0 : bmap = isl_basic_map_cow(bmap);
4903 0 : if (!bmap || !T)
4904 : goto error;
4905 :
4906 0 : n = isl_mat_cols(T);
4907 0 : if (n != isl_mat_rows(T))
4908 0 : isl_die(isl_mat_get_ctx(T), isl_error_invalid,
4909 : "expecting square matrix", goto error);
4910 :
4911 0 : total = isl_basic_map_dim(bmap, isl_dim_all);
4912 0 : if (pos + n > total || pos + n < pos)
4913 0 : isl_die(isl_mat_get_ctx(T), isl_error_invalid,
4914 : "invalid range", goto error);
4915 :
4916 0 : for (i = 0; i < bmap->n_eq; ++i)
4917 0 : if (preimage(bmap->eq[i] + 1 + pos, T) < 0)
4918 0 : goto error;
4919 0 : for (i = 0; i < bmap->n_ineq; ++i)
4920 0 : if (preimage(bmap->ineq[i] + 1 + pos, T) < 0)
4921 0 : goto error;
4922 0 : for (i = 0; i < bmap->n_div; ++i) {
4923 0 : if (isl_basic_map_div_is_marked_unknown(bmap, i))
4924 0 : continue;
4925 0 : if (preimage(bmap->div[i] + 1 + 1 + pos, T) < 0)
4926 0 : goto error;
4927 : }
4928 :
4929 0 : isl_mat_free(T);
4930 0 : return bmap;
4931 : error:
4932 0 : isl_basic_map_free(bmap);
4933 0 : isl_mat_free(T);
4934 0 : return NULL;
4935 : }
4936 :
4937 : /* Remove divs that are not strictly needed.
4938 : *
4939 : * First look for an equality constraint involving two or more
4940 : * existentially quantified variables without an explicit
4941 : * representation. Replace the combination that appears
4942 : * in the equality constraint by a single existentially quantified
4943 : * variable such that the equality can be used to derive
4944 : * an explicit representation for the variable.
4945 : * If there are no more such equality constraints, then continue
4946 : * with isl_basic_map_drop_redundant_divs_ineq.
4947 : *
4948 : * In particular, if the equality constraint is of the form
4949 : *
4950 : * f(x) + \sum_i c_i a_i = 0
4951 : *
4952 : * with a_i existentially quantified variable without explicit
4953 : * representation, then apply a transformation on the existentially
4954 : * quantified variables to turn the constraint into
4955 : *
4956 : * f(x) + g a_1' = 0
4957 : *
4958 : * with g the gcd of the c_i.
4959 : * In order to easily identify which existentially quantified variables
4960 : * have a complete explicit representation, i.e., without being defined
4961 : * in terms of other existentially quantified variables without
4962 : * an explicit representation, the existentially quantified variables
4963 : * are first sorted.
4964 : *
4965 : * The variable transformation is computed by extending the row
4966 : * [c_1/g ... c_n/g] to a unimodular matrix, obtaining the transformation
4967 : *
4968 : * [a_1'] [c_1/g ... c_n/g] [ a_1 ]
4969 : * [a_2'] [ a_2 ]
4970 : * ... = U ....
4971 : * [a_n'] [ a_n ]
4972 : *
4973 : * with [c_1/g ... c_n/g] representing the first row of U.
4974 : * The inverse of U is then plugged into the original constraints.
4975 : * The call to isl_basic_map_simplify makes sure the explicit
4976 : * representation for a_1' is extracted from the equality constraint.
4977 : */
4978 0 : __isl_give isl_basic_map *isl_basic_map_drop_redundant_divs(
4979 : __isl_take isl_basic_map *bmap)
4980 : {
4981 : int first;
4982 : int i;
4983 : unsigned o_div, n_div;
4984 : int l;
4985 : isl_ctx *ctx;
4986 : isl_mat *T;
4987 :
4988 0 : if (!bmap)
4989 0 : return NULL;
4990 0 : if (isl_basic_map_divs_known(bmap))
4991 0 : return isl_basic_map_drop_redundant_divs_ineq(bmap);
4992 0 : if (bmap->n_eq == 0)
4993 0 : return isl_basic_map_drop_redundant_divs_ineq(bmap);
4994 0 : bmap = isl_basic_map_sort_divs(bmap);
4995 0 : if (!bmap)
4996 0 : return NULL;
4997 :
4998 0 : first = isl_basic_map_first_unknown_div(bmap);
4999 0 : if (first < 0)
5000 0 : return isl_basic_map_free(bmap);
5001 :
5002 0 : o_div = isl_basic_map_offset(bmap, isl_dim_div);
5003 0 : n_div = isl_basic_map_dim(bmap, isl_dim_div);
5004 :
5005 0 : for (i = 0; i < bmap->n_eq; ++i) {
5006 0 : l = isl_seq_first_non_zero(bmap->eq[i] + o_div + first,
5007 : n_div - (first));
5008 0 : if (l < 0)
5009 0 : continue;
5010 0 : l += first;
5011 0 : if (isl_seq_first_non_zero(bmap->eq[i] + o_div + l + 1,
5012 0 : n_div - (l + 1)) == -1)
5013 0 : continue;
5014 0 : break;
5015 : }
5016 0 : if (i >= bmap->n_eq)
5017 0 : return isl_basic_map_drop_redundant_divs_ineq(bmap);
5018 :
5019 0 : ctx = isl_basic_map_get_ctx(bmap);
5020 0 : T = isl_mat_alloc(ctx, n_div - l, n_div - l);
5021 0 : if (!T)
5022 0 : return isl_basic_map_free(bmap);
5023 0 : isl_seq_cpy(T->row[0], bmap->eq[i] + o_div + l, n_div - l);
5024 0 : T = isl_mat_normalize_row(T, 0);
5025 0 : T = isl_mat_unimodular_complete(T, 1);
5026 0 : T = isl_mat_right_inverse(T);
5027 :
5028 0 : for (i = l; i < n_div; ++i)
5029 0 : bmap = isl_basic_map_mark_div_unknown(bmap, i);
5030 0 : bmap = isl_basic_map_preimage_vars(bmap, o_div - 1 + l, T);
5031 0 : bmap = isl_basic_map_simplify(bmap);
5032 :
5033 0 : return isl_basic_map_drop_redundant_divs(bmap);
5034 : }
5035 :
5036 : /* Does "bmap" satisfy any equality that involves more than 2 variables
5037 : * and/or has coefficients different from -1 and 1?
5038 : */
5039 1129738 : static int has_multiple_var_equality(__isl_keep isl_basic_map *bmap)
5040 : {
5041 : int i;
5042 : unsigned total;
5043 :
5044 1129738 : total = isl_basic_map_dim(bmap, isl_dim_all);
5045 :
5046 6238794 : for (i = 0; i < bmap->n_eq; ++i) {
5047 : int j, k;
5048 :
5049 5128768 : j = isl_seq_first_non_zero(bmap->eq[i] + 1, total);
5050 5128768 : if (j < 0)
5051 0 : continue;
5052 5149001 : if (!isl_int_is_one(bmap->eq[i][1 + j]) &&
5053 20233 : !isl_int_is_negone(bmap->eq[i][1 + j]))
5054 15549 : return 1;
5055 :
5056 5113219 : j += 1;
5057 5113219 : k = isl_seq_first_non_zero(bmap->eq[i] + 1 + j, total - j);
5058 5113219 : if (k < 0)
5059 5104383 : continue;
5060 8836 : j += k;
5061 12426 : if (!isl_int_is_one(bmap->eq[i][1 + j]) &&
5062 3590 : !isl_int_is_negone(bmap->eq[i][1 + j]))
5063 2904 : return 1;
5064 :
5065 5932 : j += 1;
5066 5932 : k = isl_seq_first_non_zero(bmap->eq[i] + 1 + j, total - j);
5067 5932 : if (k >= 0)
5068 1259 : return 1;
5069 : }
5070 :
5071 1110026 : return 0;
5072 : }
5073 :
5074 : /* Remove any common factor g from the constraint coefficients in "v".
5075 : * The constant term is stored in the first position and is replaced
5076 : * by floor(c/g). If any common factor is removed and if this results
5077 : * in a tightening of the constraint, then set *tightened.
5078 : */
5079 95064 : static __isl_give isl_vec *normalize_constraint(__isl_take isl_vec *v,
5080 : int *tightened)
5081 : {
5082 : isl_ctx *ctx;
5083 :
5084 95064 : if (!v)
5085 0 : return NULL;
5086 95064 : ctx = isl_vec_get_ctx(v);
5087 95064 : isl_seq_gcd(v->el + 1, v->size - 1, &ctx->normalize_gcd);
5088 95064 : if (isl_int_is_zero(ctx->normalize_gcd))
5089 0 : return v;
5090 95064 : if (isl_int_is_one(ctx->normalize_gcd))
5091 54592 : return v;
5092 40472 : v = isl_vec_cow(v);
5093 40472 : if (!v)
5094 0 : return NULL;
5095 40472 : if (tightened && !isl_int_is_divisible_by(v->el[0], ctx->normalize_gcd))
5096 8766 : *tightened = 1;
5097 40472 : isl_int_fdiv_q(v->el[0], v->el[0], ctx->normalize_gcd);
5098 40472 : isl_seq_scale_down(v->el + 1, v->el + 1, ctx->normalize_gcd,
5099 40472 : v->size - 1);
5100 40472 : return v;
5101 : }
5102 :
5103 : /* If "bmap" is an integer set that satisfies any equality involving
5104 : * more than 2 variables and/or has coefficients different from -1 and 1,
5105 : * then use variable compression to reduce the coefficients by removing
5106 : * any (hidden) common factor.
5107 : * In particular, apply the variable compression to each constraint,
5108 : * factor out any common factor in the non-constant coefficients and
5109 : * then apply the inverse of the compression.
5110 : * At the end, we mark the basic map as having reduced constants.
5111 : * If this flag is still set on the next invocation of this function,
5112 : * then we skip the computation.
5113 : *
5114 : * Removing a common factor may result in a tightening of some of
5115 : * the constraints. If this happens, then we may end up with two
5116 : * opposite inequalities that can be replaced by an equality.
5117 : * We therefore call isl_basic_map_detect_inequality_pairs,
5118 : * which checks for such pairs of inequalities as well as eliminate_divs_eq
5119 : * and isl_basic_map_gauss if such a pair was found.
5120 : *
5121 : * Tightening may also result in some other constraints becoming
5122 : * (rationally) redundant with respect to the tightened constraint
5123 : * (in combination with other constraints). The basic map may
5124 : * therefore no longer be assumed to have no redundant constraints.
5125 : *
5126 : * Note that this function may leave the result in an inconsistent state.
5127 : * In particular, the constraints may not be gaussed.
5128 : * Unfortunately, isl_map_coalesce actually depends on this inconsistent state
5129 : * for some of the test cases to pass successfully.
5130 : * Any potential modification of the representation is therefore only
5131 : * performed on a single copy of the basic map.
5132 : */
5133 2898920016 : __isl_give isl_basic_map *isl_basic_map_reduce_coefficients(
5134 : __isl_take isl_basic_map *bmap)
5135 : {
5136 : unsigned total;
5137 : isl_ctx *ctx;
5138 : isl_vec *v;
5139 : isl_mat *eq, *T, *T2;
5140 : int i;
5141 : int tightened;
5142 :
5143 2898920016 : if (!bmap)
5144 0 : return NULL;
5145 2898920016 : if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_REDUCED_COEFFICIENTS))
5146 209 : return bmap;
5147 2898919807 : if (isl_basic_map_is_rational(bmap))
5148 2897679204 : return bmap;
5149 1240603 : if (bmap->n_eq == 0)
5150 110865 : return bmap;
5151 1129738 : if (!has_multiple_var_equality(bmap))
5152 1110026 : return bmap;
5153 :
5154 19712 : total = isl_basic_map_dim(bmap, isl_dim_all);
5155 19712 : ctx = isl_basic_map_get_ctx(bmap);
5156 19712 : v = isl_vec_alloc(ctx, 1 + total);
5157 19712 : if (!v)
5158 0 : return isl_basic_map_free(bmap);
5159 :
5160 19712 : eq = isl_mat_sub_alloc6(ctx, bmap->eq, 0, bmap->n_eq, 0, 1 + total);
5161 19712 : T = isl_mat_variable_compression(eq, &T2);
5162 19712 : if (!T || !T2)
5163 : goto error;
5164 19712 : if (T->n_col == 0) {
5165 0 : isl_mat_free(T);
5166 0 : isl_mat_free(T2);
5167 0 : isl_vec_free(v);
5168 0 : return isl_basic_map_set_to_empty(bmap);
5169 : }
5170 :
5171 19712 : bmap = isl_basic_map_cow(bmap);
5172 19712 : if (!bmap)
5173 0 : goto error;
5174 :
5175 19712 : tightened = 0;
5176 114776 : for (i = 0; i < bmap->n_ineq; ++i) {
5177 95064 : isl_seq_cpy(v->el, bmap->ineq[i], 1 + total);
5178 95064 : v = isl_vec_mat_product(v, isl_mat_copy(T));
5179 95064 : v = normalize_constraint(v, &tightened);
5180 95064 : v = isl_vec_mat_product(v, isl_mat_copy(T2));
5181 95064 : if (!v)
5182 0 : goto error;
5183 95064 : isl_seq_cpy(bmap->ineq[i], v->el, 1 + total);
5184 : }
5185 :
5186 19712 : isl_mat_free(T);
5187 19712 : isl_mat_free(T2);
5188 19712 : isl_vec_free(v);
5189 :
5190 19712 : ISL_F_SET(bmap, ISL_BASIC_MAP_REDUCED_COEFFICIENTS);
5191 :
5192 19712 : if (tightened) {
5193 7523 : int progress = 0;
5194 :
5195 7523 : ISL_F_CLR(bmap, ISL_BASIC_MAP_NO_REDUNDANT);
5196 7523 : bmap = isl_basic_map_detect_inequality_pairs(bmap, &progress);
5197 7523 : if (progress) {
5198 218 : bmap = eliminate_divs_eq(bmap, &progress);
5199 218 : bmap = isl_basic_map_gauss(bmap, NULL);
5200 : }
5201 : }
5202 :
5203 19712 : return bmap;
5204 : error:
5205 0 : isl_mat_free(T);
5206 0 : isl_mat_free(T2);
5207 0 : isl_vec_free(v);
5208 0 : return isl_basic_map_free(bmap);
5209 : }
5210 :
5211 : /* Shift the integer division at position "div" of "bmap"
5212 : * by "shift" times the variable at position "pos".
5213 : * "pos" is as determined by isl_basic_map_offset, i.e., pos == 0
5214 : * corresponds to the constant term.
5215 : *
5216 : * That is, if the integer division has the form
5217 : *
5218 : * floor(f(x)/d)
5219 : *
5220 : * then replace it by
5221 : *
5222 : * floor((f(x) + shift * d * x_pos)/d) - shift * x_pos
5223 : */
5224 0 : __isl_give isl_basic_map *isl_basic_map_shift_div(
5225 : __isl_take isl_basic_map *bmap, int div, int pos, isl_int shift)
5226 : {
5227 : int i;
5228 : unsigned total;
5229 :
5230 0 : if (isl_int_is_zero(shift))
5231 0 : return bmap;
5232 0 : if (!bmap)
5233 0 : return NULL;
5234 :
5235 0 : total = isl_basic_map_dim(bmap, isl_dim_all);
5236 0 : total -= isl_basic_map_dim(bmap, isl_dim_div);
5237 :
5238 0 : isl_int_addmul(bmap->div[div][1 + pos], shift, bmap->div[div][0]);
5239 :
5240 0 : for (i = 0; i < bmap->n_eq; ++i) {
5241 0 : if (isl_int_is_zero(bmap->eq[i][1 + total + div]))
5242 0 : continue;
5243 0 : isl_int_submul(bmap->eq[i][pos],
5244 : shift, bmap->eq[i][1 + total + div]);
5245 : }
5246 0 : for (i = 0; i < bmap->n_ineq; ++i) {
5247 0 : if (isl_int_is_zero(bmap->ineq[i][1 + total + div]))
5248 0 : continue;
5249 0 : isl_int_submul(bmap->ineq[i][pos],
5250 : shift, bmap->ineq[i][1 + total + div]);
5251 : }
5252 0 : for (i = 0; i < bmap->n_div; ++i) {
5253 0 : if (isl_int_is_zero(bmap->div[i][0]))
5254 0 : continue;
5255 0 : if (isl_int_is_zero(bmap->div[i][1 + 1 + total + div]))
5256 0 : continue;
5257 0 : isl_int_submul(bmap->div[i][1 + pos],
5258 : shift, bmap->div[i][1 + 1 + total + div]);
5259 : }
5260 :
5261 0 : return bmap;
5262 : }
|
