1/* 2 * Copyright 2008-2009 Katholieke Universiteit Leuven 3 * Copyright 2012-2013 Ecole Normale Superieure 4 * Copyright 2014-2015 INRIA Rocquencourt 5 * Copyright 2016 Sven Verdoolaege 6 * 7 * Use of this software is governed by the MIT license 8 * 9 * Written by Sven Verdoolaege, K.U.Leuven, Departement 10 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium 11 * and Ecole Normale Superieure, 45 rue d���Ulm, 75230 Paris, France 12 * and Inria Paris - Rocquencourt, Domaine de Voluceau - Rocquencourt, 13 * B.P. 105 - 78153 Le Chesnay, France 14 */ 15 16#include <isl_ctx_private.h> 17#include <isl_map_private.h> 18#include "isl_equalities.h" 19#include <isl/map.h> 20#include <isl_seq.h> 21#include "isl_tab.h" 22#include <isl_space_private.h> 23#include <isl_mat_private.h> 24#include <isl_vec_private.h> 25 26#include <bset_to_bmap.c> 27#include <bset_from_bmap.c> 28#include <set_to_map.c> 29#include <set_from_map.c> 30 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-open 116 * 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 == 0 118 * corresponds to the constant term. 119 * 120 * That is, the integer division is of the form 121 * 122 * floor((... + (c * d + r) * x_pos + ...)/d) 123 * 124 * with -d < 2 * r <= d. 125 * Replace it by 126 * 127 * floor((... + r * x_pos + ...)/d) + c * x_pos 128 * 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 and 158 * 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) of 178 * integer division "div", if needed. 179 * In particular, make sure all coefficients lie in 180 * 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) of 208 * the known integer divisions, if needed 209 * In particular, make sure all coefficients lie in 210 * 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 form 236 * 237 * floor((a + m f(x))/(m d)) 238 * 239 * then replace it by 240 * 241 * floor((floor(a/m) + f(x))/d) 242 * 243 * The difference {a/m}/d in the argument satisfies 0 <= {a/m}/d < 1/d 244 * 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 form 273 * 274 * floor((a + m f(x))/(m d)) 275 * 276 * and replace it by 277 * 278 * floor((floor(a/m) + f(x))/d) 279 * 280 * The difference {a/m}/d in the argument satisfies 0 <= {a/m}/d < 1/d 281 * 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 k 346 * if the equality contains any divs. 347 * If keep_divs is set, then the divs have been ordered 348 * and we can keep the definition as long as the result 349 * 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 } else 358 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 not 384 * 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 equalities 414 */ 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 inequalities 456 */ 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" involve 494 * 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 constraint 523 * at position "eq" in "bmap" is the local variable at position "div". 524 * It can therefore be used to extract an explicit representation 525 * for that variable. 526 * Do so unless the local variable already has an explicit representation or 527 * the explicit representation would involve any other local variables 528 * that in turn do not have an explicit representation. 529 * An equality constraint involving local variables without an explicit 530 * representation can be used in isl_basic_map_drop_redundant_divs 531 * to separate out an independent local variable. Introducing 532 * 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 form 537 * 538 * f(x) + n e >= 0 539 * 540 * with n a positive number. The explicit representation derived from 541 * this constraint is 542 * 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 equality 578 * constraints of "bmap". 579 * That is, put them into row-echelon form, starting from the last column 580 * backward and use them to eliminate the corresponding coefficients 581 * from all constraints. 582 * 583 * If "progress" is not NULL, then it gets set if the elimination 584 * results in any changes. 585 * The elimination process may result in some equality constraints 586 * getting interchanged or removed. 587 * If "swap" or "drop" are not NULL, then they get called when 588 * two equality constraints get interchanged or 589 * 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 equality 592 * 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, some 688 * of which are NULL. The inequalities are hashed on the coefficients 689 * except the constant term. 690 * "size" is the number of elements in the array and is always a power of two 691 * "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 of 734 * 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 where 737 * 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 of 751 * 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 very 753 * same inequality. 754 * If there is no such inequality, then return the position where 755 * 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 respect 788 * to the constraints in "ci"? 789 * 790 * Look for an inequality in "ci" with the same coefficients and then 791 * check if the contant term of "ineq" is greater than or equal 792 * to the constant term of that inequality. If so, "ineq" is clearly 793 * redundant. 794 * 795 * Note that hash_index_ineq ignores a stored constraint if it has 796 * the same address as the passed inequality. It is ok to pass 797 * the address of a local variable here since it will never be 798 * 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 make 812 * sure we do it from last div to first div, in order not to 813 * change the position of the other divs that still need to 814 * 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 equalities 919 * of the form 920 * 921 * a x = m * e_i 922 * 923 * and we want to replace the set of e_i by a minimal set and 924 * such that the new e_i have a canonical representation in terms 925 * of the vector x. 926 * If any of the equalities involves more than one divs, then 927 * we currently simply bail out. 928 * 929 * Let us first additionally assume that all equalities involve 930 * a div. The equalities then express modulo constraints on the 931 * remaining variables and we can use "parameter compression" 932 * to find a minimal set of constraints. The result is a transformation 933 * 934 * x = T(x') = x_0 + G x' 935 * 936 * with G a lower-triangular matrix with all elements below the diagonal 937 * 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 affine 939 * T matrix. 940 * The rows i of G with a 1 on the diagonal do not impose any modulo 941 * constraint and simply express x_i = x'_i. 942 * For each of the remaining rows i, we introduce a div and a corresponding 943 * equality. In particular 944 * 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 the 948 * corresponding div (if g_kk != 1). 949 * 950 * If there are any equalities not involving any div, then we 951 * first apply a variable compression on the variables x: 952 * 953 * x = C x'' x'' = C_2 x 954 * 955 * and perform the above parameter compression on A C instead of on A. 956 * The resulting compression is then of the form 957 * 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 else 1066 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 else 1083 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 else 1091 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, we 1133 * do not allow such a definition if the resulting expression would refer to 1134 * any other undefined divs or if any known div is defined in 1135 * 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 inequality 1171 * (disregarding the div that it would define) is in an earlier position 1172 * 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 them 1197 * to obtain an expression for one of the hitherto unknown divs or 1198 * 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 one 1201 * of the divs, then this pair can be used to define the div. 1202 * To avoid the introduction of circular definitions of divs, we 1203 * do not use the pair if the resulting expression would refer to 1204 * any other undefined divs or if any known div is defined in 1205 * 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 else 1231 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 these 1286 * changes since the contents of the hash 1287 * 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 } else 1295 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 integral 1326 * (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 elimination 1329 * results in any changes. 1330 * 1331 * That is, replace 1332 * 1333 * floor(e/m) + f >= 0 1334 * 1335 * by 1336 * 1337 * e + m f >= 0 1338 * 1339 * and 1340 * 1341 * -floor(e/m) + f >= 0 1342 * 1343 * by 1344 * 1345 * -e + m f + m - 1 >= 0 1346 * 1347 * The first conversion is valid because floor(e/m) >= -f is equivalent 1348 * to e/m >= -f because -f is an integral expression. 1349 * The second conversion follows from the fact that 1350 * 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 eliminated 1355 * due to being redundant with respect to the constraint that is 1356 * being modified by this function. The modified constraint may 1357 * no longer imply this div constraint, so we add it back to make 1358 * 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 else 1391 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 with 1411 * 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 elimination 1414 * results in any changes. 1415 * 1416 * We skip integral divs, i.e., those with denominator 1, as we would 1417 * risk eliminating the div from the div constraints. We do not need 1418 * to handle those divs here anyway since the div constraints will turn 1419 * out to form an equality and this equality can then be used to eliminate 1420 * 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 every 1453 * 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 with 1461 * a (positive or negative) unit coefficient. 1462 * If "progress" is not NULL, then it gets set if the elimination 1463 * 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 selects 1472 * integer divisions that only appear with 1473 * a (positive or negative) unit coefficient 1474 * (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 with 1506 * a (positive or negative) unit coefficient, 1507 * but only if they do not appear in any other constraints 1508 * (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 } else 1586 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 constraints 1593 * 1594 * f - m d >=0 1595 * -(f - (m - 1)) + m d >= 0 1596 * 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 or 1635 * when the divs of a basic map have been temporarily aligned 1636 * 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 redundant 1664 * 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 as 1693 * 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 using 1716 * 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 the 1821 * 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 local 1860 * variables (directly or indirectly). 1861 * 1862 * Dropping all constraints on a local variable will make it redundant, 1863 * so it will get removed implicitly by 1864 * isl_basic_map_drop_constraints_involving_dims. Some other local 1865 * variables may also end up becoming redundant if they only appear 1866 * in constraints together with the unknown local variable. 1867 * Therefore, start over after calling 1868 * 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 local 1907 * 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 local 1920 * 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_divs 1956 * 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 entry 2051 * in "row" with -1. 2052 * If an inequality only has a non-negative constant term, then 2053 * 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 constraints 2134 * in "context" or that are more relaxed (greater constant term). 2135 * 2136 * We perform the test for shifted copies on the pure constraints 2137 * 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 of 2188 * 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 group 2226 * of the corresponding variable by the minimal group of all 2227 * the variables involved in the constraint considered so far 2228 * (if this minimum is smaller) or replace the minimum by this group 2229 * (if the minimum is larger). 2230 * 2231 * At the end, all the variables in "c" will (indirectly) point 2232 * 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 } else 2251 group[group[j]] = min; 2252 } 2253} 2254 2255/* Allocate an array of groups of variables, one for each variable 2256 * 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 are 2271 * not related to any of the variables marked with a "-1" in "group". 2272 * 2273 * We construct groups of variables that collect variables that 2274 * (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 identified 2277 * 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 the 2283 * group of the variables in the constraint by the smallest group. 2284 * Finally, we resolve indirect references to groups by running over 2285 * the variables. 2286 * 2287 * After computing the groups, we drop constraints that do not involve 2288 * 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 computing 2329 * the gist of "bset". 2330 * 2331 * In particular, drop constraints in variables that are not related 2332 * 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 belonging 2336 * 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 computing 2373 * the gist of the inequalities "ineq". 2374 * Inequalities in "ineq" for which the corresponding element of row 2375 * 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 related 2378 * 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 belonging 2382 * 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 that 2433 * the corresponding inequality of "bset" is redundant, or the index 2434 * 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 equality 2439 * in the tableau, then turn the inequality into an equality and 2440 * 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 context 2499 * 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 taken 2503 * into account by the transformation such that only the inequalities 2504 * 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 may 2509 * 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 become 2514 * irrelevant for computing the gist of "bset". 2515 * Note that this removal of constraints cannot be replaced by 2516 * a factorization because factors in "bset" may still be connected 2517 * to each other through constraints in "context". 2518 * 2519 * If there are any inequalities left, we construct a tableau for 2520 * the context and then add the inequalities of "bset". 2521 * Before adding these inequalities, we freeze all constraints such that 2522 * they won't be considered redundant in terms of the constraints of "bset". 2523 * Then we detect all redundant constraints (among the 2524 * constraints that weren't frozen), first by checking for redundancy in the 2525 * the tableau and then by checking if replacing a constraint by its negation 2526 * would lead to an empty set. This last step is fairly expensive 2527 * 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 context 2645 * of "context", for the case where both "bset" and "context" are 2646 * 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 context 2670 * of "context", for the case where the combined equalities of 2671 * "bset" and "context" allow for a compression that can be obtained 2672 * by preapplication of "T". 2673 * If the compression of "context" is empty, meaning that "bset" and 2674 * "context" do not intersect, then return the empty set. 2675 * 2676 * "bset" itself is not transformed by "T". Instead, the inequalities 2677 * are extracted from "bset" and those are transformed by "T". 2678 * uset_gist_full then determines which of the transformed inequalities 2679 * are redundant with respect to the transformed "context" and removes 2680 * the corresponding inequalities from "bset". 2681 * 2682 * After preapplying "T" to the inequalities, any common factor is 2683 * removed from the coefficients. If this results in a tightening 2684 * of the constant term, then the same tightening is applied to 2685 * the corresponding untransformed inequality in "bset". 2686 * That is, if after plugging in T, a constraint f(x) >= 0 is of the form 2687 * 2688 * g f'(x) + r >= 0 2689 * 2690 * with 0 <= r < g, then it is equivalent to 2691 * 2692 * f'(x) >= 0 2693 * 2694 * This means that f(x) >= 0 is equivalent to f(x) - r >= 0 in the affine 2695 * subspace compressed by T since the latter would be transformed to 2696 * 2697 * g f'(x) >= 0 2698 */ 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 context 2783 * of context. In particular, equalities that are linear combinations 2784 * of those in context are removed. Then the inequalities that are 2785 * redundant in the context of the equalities and inequalities of 2786 * 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 hulls 2793 * of "bset" and "context", 2794 * compute the gist inside this intersection and then reduce 2795 * the constraints with respect to the equalities of the context 2796 * 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_full 2801 * will remove the second of those constraints. We therefore first 2802 * sort the constraints so that constraints not involving existentially 2803 * 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 involve 2870 * local variables. This function assumes that Gaussian elimination 2871 * 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 equality 2936 * 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 elimination 2943 * has been applied to the equality constraints in both "bmap1" and "bmap2" 2944 * such that the equality constraints not involving local variables 2945 * 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 constraints 2949 * in the other basic map. 2950 * Otherwise, combine the equality constraints from both into a new 2951 * basic map such that Gaussian elimination can be applied to this combination 2952 * and then construct a variable compression from the resulting 2953 * 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", compressed 2991 * with respect to both the stride constraints in "context" and 2992 * 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 local 2997 * variables. First compute a variable compression 2998 * 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) = 0 3005 * 3006 * with y the local variables and plug in the variable compression, 3007 * resulting in 3008 * 3009 * A(V x') + B(y) = 0 3010 * 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) = 0 3018 * 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 else 3052 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 that 3062 * is strictly smaller than the exponent in "c". 3063 * All exponents in *g are known to be smaller than or equal 3064 * to those in "c". 3065 * 3066 * That is, if *g is equal to 3067 * 3068 * p_1^{e_1} p_2^{e_2} ... p_n^{e_n} 3069 * 3070 * and "c" is equal to 3071 * 3072 * p_1^{f_1} p_2^{f_2} ... p_n^{f_n} 3073 * 3074 * then update *g to 3075 * 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 therefore 3080 * 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 positive 3082 * 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 exponent 3084 * of the prime factors that need to be removed, while leaving the 3085 * other prime factors untouched. 3086 * Repeating this process until gcd(*g, c / *g) = 1 therefore 3087 * 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 exploits 3106 * all equalities in the context and the other equalities in "bmap". 3107 * 3108 * If the stride constraints of "bmap" are of the form 3109 * 3110 * C(x) + D(y) = 0 3111 * 3112 * then A is of the form 3113 * 3114 * B(x') + D(y) = 0 3115 * 3116 * If any of these constraints involves only a single local variable y, 3117 * then the constraint appears as 3118 * 3119 * f(x) + m y_i = 0 3120 * 3121 * in "bmap" and as 3122 * 3123 * h(x') + m y_i = 0 3124 * 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 and 3129 * 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 known 3135 * to be a multiple of g and therefore in particular of this power 3136 * of the prime factors. 3137 * Prime factors that appear with a smaller exponent in g cannot 3138 * be removed from m. 3139 * Let g' be the divisor of g containing all prime factors that 3140 * appear with the same exponent in m and g, then 3141 * 3142 * f(x) + m y_i = 0 3143 * 3144 * can be replaced by 3145 * 3146 * f(x) + m/g' y_i' = 0 3147 * 3148 * Note that (if g' != 1) this changes the explicit representation 3149 * of y_i to that of y_i', so the integer division at position i 3150 * is marked unknown and later recomputed by a call to 3151 * 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 on 3206 * the remaining equality constraints in "bmap" and all equality 3207 * 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 compressed 3211 * space exploiting all the other equality constraints and then 3212 * 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 need 3245 * to (re)introduce the div constraints. If we happen to have 3246 * discovered that some of these integer divisions are equal to 3247 * some affine combination of other variables, then these div 3248 * constraints may end up getting simplified in terms of the equalities, 3249 * resulting in extra inequalities on the other variables that 3250 * may have been removed already or that may not even have been 3251 * part of the input. We try and remove those constraints of 3252 * this form that are most obviously redundant with respect to 3253 * the context. We also remove those div constraints that are 3254 * redundant with respect to the other constraints in the result. 3255 * 3256 * The stride constraints among the equality constraints in "bmap" are 3257 * also simplified with respecting to the other equality constraints 3258 * 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 and 3379 * the initial local variables of "bmap" are assumed to be the same 3380 * as those of "context". 3381 * The constraints of both "bmap" and "context" are assumed 3382 * 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 and 3438 * the initial local variables of "bmap" are assumed to be the same 3439 * 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 representations 3494 * for all local variables. 3495 * 3496 * First align the divs of "bmap" to those of "context" and 3497 * 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 representations 3560 * 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 representations 3610 * 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 and 3623 * 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 it 3652 * to simplify "map" 3653 * If "map" and "context" are identical to each other, then we can 3654 * 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 procedure 3660 * 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 disjunct 3665 * then we need to approximate the context by a single disjunct set. 3666 * Simply taking the simple hull may drop constraints that are 3667 * only implicitly available in each disjunct. We therefore also 3668 * look for constraints among those defining "map" that are valid 3669 * for the context. These can then be used to simplify away 3670 * 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 of 3807 * one basic map in the context of the equalities of the other 3808 * 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 information 3905 * 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 spaces 3908 * (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 (ignoring 3941 * 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 will 3947 * not be able to detect if they are disjoint. 3948 * 3949 * Otherwise we check if each basic map in "map1" is obviously disjoint 3950 * 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 them 3988 * 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 are 3990 * not disjoint. 3991 * If none of these cases apply, then check if all pairs of basic maps 3992 * 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 them 4030 * 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 if 4032 * 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 constraints 4098 * "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 upper 4120 * 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 that 4125 * - the bounds are opposites of each other (except for the constant 4126 * term) 4127 * - the bounds do not reference any other div 4128 * - no div is defined in terms of this div 4129 * 4130 * Let m be the size of the range allowed on the div by the bounds. 4131 * That is, the bounds are of the form 4132 * 4133 * e <= a <= e + m - 1 4134 * 4135 * with e some expression in the other variables. 4136 * We look for another div b such that no third div is defined in terms 4137 * of this second div b and such that in any constraint that contains 4138 * a (except for the given lower and upper bound), also contains b 4139 * 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 form 4142 * 4143 * e + f (a + m b) >= 0 4144 * 4145 * Furthermore, in the constraints that only contain b, the coefficient 4146 * of b should be equal to 1 or -1. 4147 * If so, we return b so that "a + m b" can be replaced by 4148 * 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 of 4238 * an inequality that ensures that a pair of bounds always allows 4239 * for an integer value. 4240 * 4241 * "tab" is the tableau in which the inequality is evaluated. It may 4242 * 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 and 4245 * 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 allow 4286 * for an integer value of the given div? 4287 * Determine this property by constructing an inequality 4288 * 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 be 4293 * 4294 * -n_u a + e_u >= 0 4295 * 4296 * and the lower bound 4297 * 4298 * n_l a + e_l >= 0 4299 * 4300 * Let n_u = f_u g and n_l = f_l g, with g = gcd(n_u, n_l). 4301 * We have 4302 * 4303 * - f_u e_l <= f_u f_l g a <= f_l e_u 4304 * 4305 * Since all variables are integer valued, this is equivalent to 4306 * 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 least 4310 * one integer value for a. 4311 * That is, the test constraint is 4312 * 4313 * f_l e_u + f_u e_l + f_l - 1 + f_u - 1 + 1 >= f_u f_l g 4314 * 4315 * or 4316 * 4317 * f_l e_u + f_u e_l + f_l - 1 + f_u - 1 + 1 - f_u f_l g >= 0 4318 * 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 by 4322 * 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 pair 4324 * of constraints is 4325 * 4326 * f_l e_u + f_u e_l >= 0 4327 * 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 constant 4330 * term of the scaled down test constraint, then the test constraint 4331 * 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, then 4335 * 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 div 4382 * are such that they allow at least one integer value of the div, 4383 * then we can eliminate the div using Fourier-Motzkin without 4384 * 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 coefficient 4389 * 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 l 4479 * and the upper bound u, div1 always occurs together with div2 in the form 4480 * (div1 + m div2), where m is the constant range on the variable div1 4481 * allowed by l and u, replace the pair div1 and div2 by a single 4482 * 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 contain 4486 * div2 = (div - div1) / m 4487 * 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 appear 4490 * 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 be 4495 * 4496 * div1 + f >=0 4497 * 4498 * and u 4499 * 4500 * -div1 + f' >= 0 4501 * 4502 * A lower bound on div2 4503 * 4504 * div2 + t >= 0 4505 * 4506 * can be replaced by 4507 * 4508 * m div2 + div1 + m t + f >= 0 4509 * 4510 * An upper bound 4511 * 4512 * -div2 + t >= 0 4513 * 4514 * can be replaced by 4515 * 4516 * -(m div2 + div1) + m t + f' >= 0 4517 * 4518 * These constraint are those that we would obtain from eliminating 4519 * div1 using Fourier-Motzkin. 4520 * 4521 * After all constraints have been modified, we drop the lower and upper 4522 * bound and then drop div1. 4523 * Since the new div is only placed in the same location that used 4524 * to store div2, but otherwise has a different meaning, any possible 4525 * 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 else 4557 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 and 4579 * then continue with dropping more redundant divs. 4580 * 4581 * We loop over all pairs of lower and upper bounds on a div 4582 * with coefficient 1 and -1, respectively, check if there 4583 * is any other div "c" with which we can coalesce the div 4584 * 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 constraints 4634 * "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" has 4673 * been modified, simplying it if "simplify" is set. 4674 * Free the temporary data structure "pairs" that was associated 4675 * 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 variable 4687 * 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 constraint 4723 * "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 and 4737 * then try and drop redundant divs again, 4738 * freeing the temporary data structure "pairs" that was associated 4739 * 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 two 4750 * inequality constraints "ineq1" and "ineq2" in which it appears 4751 * from "bmap" and then try and drop redundant divs again, 4752 * freeing the temporary data structure "pairs" that was associated 4753 * 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 constraints 4771 * 4772 * f(x) + n d + c >= 0, (ineq) 4773 * 4774 * with d the variable at position "pos", and 4775 * 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, store 4781 * 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 constraints 4795 * 4796 * f(x) + n d + c >= 0, (ineq) 4797 * 4798 * with d the variable at position "pos", and 4799 * 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, store 4805 * 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 form 4822 * 4823 * f(x) + n d + c >= 0 4824 * 4825 * with n > 0, c the constant term and 4826 * d the existentially quantified variable "div". 4827 * That is, the lower bound is 4828 * 4829 * ceil((-f(x) - c)/n) 4830 * 4831 * Look for a pair of constraints 4832 * 4833 * f(x) + c0 >= 0 4834 * -f(x) + c1 >= 0 4835 * 4836 * i.e., -c1 <= -f(x) <= c0, that fix ceil((-f(x) - c)/n) to a constant value. 4837 * That is, check that 4838 * 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 quantified 4897 * variable "div", such that the corresponding lower bound has 4898 * a fixed value in "bmap", assign this fixed value to the variable and 4899 * then try and drop redundant divs again, 4900 * freeing the temporary data structure "pairs" that was associated 4901 * to the old version of "bmap". 4902 * "lower" determines the constant value for the lower bound. 4903 * 4904 * In particular, "ineq" is of the form 4905 * 4906 * f(x) + n d + c >= 0, 4907 * 4908 * while "lower" is of the form 4909 * 4910 * f(x) + c0 >= 0 4911 * 4912 * The lower bound is ceil((-f(x) - c)/n) and its constant value 4913 * 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 later 4936 * 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 inequality 4964 * 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 moreover 4968 * those two constraints are opposite to each other, except for the constant 4969 * term and if the sum of the constant terms is such that for any value 4970 * of the other values, there is always at least one integer value of the 4971 * div, i.e., if one plus this sum is greater than or equal to 4972 * 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 explicit 4976 * representation, appears in only a single lower bound that does not 4977 * involve any other such existentially quantified variables and appears 4978 * 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 value 4982 * 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 may 4988 * in general introduce an extra integer division, in which case 4989 * it may be better to pick another value. 4990 * If this integer division has a known constant value, then plugging 4991 * in this constant value removes the existentially quantified variable 4992 * completely. In particular, if the lower bound is of the form 4993 * ceil((-f(x) - c)/n) and there are two constraints, f(x) + c0 >= 0 and 4994 * -f(x) + c1 >= 0 such that ceil((-c1 - c)/n) = ceil((c0 - c)/n), 4995 * then the existentially quantified variable can be assigned this 4996 * 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 least 5000 * 4 constraints, but the constraints may have been simplified. 5001 * 5002 * If any divs are left after these simple checks then we move on 5003 * 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 else 5070 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 else 5115 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 replace 5138 * 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 more 5211 * existentially quantified variables without an explicit 5212 * representation. Replace the combination that appears 5213 * in the equality constraint by a single existentially quantified 5214 * variable such that the equality can be used to derive 5215 * an explicit representation for the variable. 5216 * If there are no more such equality constraints, then continue 5217 * with isl_basic_map_drop_redundant_divs_ineq. 5218 * 5219 * In particular, if the equality constraint is of the form 5220 * 5221 * f(x) + \sum_i c_i a_i = 0 5222 * 5223 * with a_i existentially quantified variable without explicit 5224 * representation, then apply a transformation on the existentially 5225 * quantified variables to turn the constraint into 5226 * 5227 * f(x) + g a_1' = 0 5228 * 5229 * with g the gcd of the c_i. 5230 * In order to easily identify which existentially quantified variables 5231 * have a complete explicit representation, i.e., without being defined 5232 * in terms of other existentially quantified variables without 5233 * an explicit representation, the existentially quantified variables 5234 * are first sorted. 5235 * 5236 * The variable transformation is computed by extending the row 5237 * [c_1/g ... c_n/g] to a unimodular matrix, obtaining the transformation 5238 * 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 explicit 5247 * 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 variables 5311 * 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 replaced 5352 * by floor(c/g). If any common factor is removed and if this results 5353 * 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 involving 5380 * more than 2 variables and/or has coefficients different from -1 and 1, 5381 * then use variable compression to reduce the coefficients by removing 5382 * 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 and 5385 * 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 of 5391 * the constraints. If this happens, then we may end up with two 5392 * 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_eq 5395 * and isl_basic_map_gauss if such a pair was found. 5396 * 5397 * Tightening may also result in some other constraints becoming 5398 * (rationally) redundant with respect to the tightened constraint 5399 * (in combination with other constraints). The basic map may 5400 * 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 state 5405 * for some of the test cases to pass successfully. 5406 * Any potential modification of the representation is therefore only 5407 * 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 == 0 5496 * corresponds to the constant term. 5497 * 5498 * That is, if the integer division has the form 5499 * 5500 * floor(f(x)/d) 5501 * 5502 * then replace it by 5503 * 5504 * floor((f(x) + shift * d * x_pos)/d) - shift * x_pos 5505 */ 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