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